# Doubly Robust Calibration of Prediction Sets under Covariate Shift

Yachong Yang Email and address<EMAIL_ADDRESS>Academic Research
Building, 265 S 37th St, Philadalphia, PA, US. Department of Statistics,
University of Pennsylvania  Arun Kumar Kuchibhotla Email and address:
<EMAIL_ADDRESS>Baker Hall, 4909 Frew St, Pittsburgh, PA, US. Department of
Statistics & Data Science, Carnegie Mellon University  Eric Tchetgen Tchetgen
Email and address<EMAIL_ADDRESS>Academic Research Building, 265 S
37th St, Philadalphia, PA, US. Department of Statistics, University of
Pennsylvania

###### Abstract

Conformal prediction has received tremendous attention in recent years and has
offered new solutions to problems in missing data and causal inference; yet
these advances have not leveraged modern semiparametric efficiency theory for
more robust and efficient uncertainty quantification. In this paper, we
consider the problem of obtaining distribution-free prediction regions
accounting for a shift in the distribution of the covariates between the
training and test data. Under an explainable covariate shift assumption
analogous to the standard missing at random assumption, we propose three
variants of a general framework to construct well-calibrated prediction
regions for the unobserved outcome in the test sample. Our approach is based
on the efficient influence function for the quantile of the unobserved outcome
in the test population combined with an arbitrary machine learning prediction
algorithm, without compromising asymptotic coverage. We establish that the
resulting prediction sets eventually attain nominal coverage in large samples.
This guarantee is a consequence of the product bias form of our proposal which
implies correct coverage if either the propensity score or the conditional
distribution of the response is estimated sufficiently well. Our results also
provide a framework for construction of doubly robust prediction sets of
individual treatment effects, under the unconfoundedness condition. We further
discuss aggregation of prediction sets from different machine learning
algorithms for optimal prediction and illustrate the performance of our
methods in both synthetic and real data. Finally, inspired by sensitivity
analysis in missing data, we briefly discuss how our proposal could be
extended to account for departures from the explainable covariate shift
setting.

###### Abstract

This supplement contains the proofs to all the main results in the paper and
some supporting lemmas.

##  1 Introduction

Prediction is a major focus of modern machine learning literature. Most
machine learning methods are designed for point prediction, but accurately
quantifying the uncertainty associated with a given point prediction algorithm
remains an important challenge in many applications. Given independent and
identically distributed (i.i.d.) pairs $\left(X_{i},Y_{i}\right)$,
$i=1,\ldots,N$, from a distribution $P=P_{X}\otimes P_{Y|X}$ supported on
$\mathcal{X}\times\mathbb{R}$ (e.g., $\mathcal{X}=\mathbb{R}^{d},d\geq 1$),
and given a desired nominal coverage rate $1-\alpha\in(0,1)$, the goal of
prediction with well-calibrated uncertainty quantification is to build a
prediction set $\widehat{C}_{N,\alpha}$, such that

$\mathbb{P}\bigl{(}Y_{f}\in\widehat{C}_{N,\alpha}(X_{f})\bigr{)}\geq
1-\alpha,$ (1)

where the probability is taken over the marginal distribution of all the
training data along with $(X_{f},Y_{f})$. Note that (1) does not imply
conditional coverage
$\mathbb{P}\bigl{(}Y_{f}\in\widehat{C}_{N,\alpha}(X_{f})|X_{f}=x_{f}\bigr{)}\geq
1-\alpha$, which is known to be impossible without assumptions over the
underlying distribution as shown in Barber et al. (2019). This goal, however,
can be _approximately_ achieved where _approximately_ is meant either
asymptotically or by conditioning on $X_{f}$ belonging to a set $A$ rather
than $X_{f}$ being equal to, say $x_{f}$. Conformal prediction introduced by
Vovk et al. (2005) provides a simple and finite-sample valid solution to (1)
without any assumption on the distribution $P$, requiring only that the
training data and $(X_{f},Y_{f})$ to be exchangeable (jointly). It provides a
valid prediction set by wrapping around any point prediction algorithm,
irrespective of what the point prediction algorithm is.

Recently several works have considered the extension of conformal prediction
methodology to the case of non-exchangeable data; see Section 2 for a review.
Our work lies in this space. One important distinction from the exchangeable
case is that finite sample coverage guarantees are generally impossible for
non-exchangeable data without very restrictive assumptions. Formally, we
consider the problem of prediction with uncertainty quantification under
covariate shift where the completely observed training data is drawn i.i.d.
from $P_{X}\otimes P_{Y|X}$ but the test point that needs to be covered by a
prediction set is drawn from $Q_{X}\otimes P_{Y|X}$. This problem, which we
term explainable covariate shift problem, was first posed by Tibshirani et al.
(2019) who provides a solution called “Weighted Conformal Prediction” (WCP)
which assumes that the covariate shift is known via $dQ_{X}/dP_{X}$. This
method has been extended by Lei and Candès (2020) to allow for unknown
covariate shift.

The prediction problem under covariate shift can be equivalently stated as
predicting the label for a given feature vector observed from a different
covariate distribution. In terms of real-world applications, prediction under
covariate shift is important in semi-supervised and transfer learning
settings. In health care and related problems, it is often the case that the
amount of labeled data is limited in comparison to unlabeled data. This is,
particularly, true with electronic health record (EHR) data where labeling the
response can be costly and/or laborious (Chakrabortty and Cai, 2018). Instead
of assuming the data are missing completely at random (i.e., observations are
chosen to be labeled completely at random), it may be preferable to allow for
the possibility that they are labeled based on observed covariate features.
Prediction in this case is same as prediction under covariate shift as proved
in Section 3.2 and can be useful either for imputation or for understanding
the spread in the response distribution. Recently, Lei and Candès (2020) and
Jin et al. (2021) showed how to use prediction under covariate shift to
construct prediction intervals for individual treatment effect under
traditional causal inference assumptions. This can potentially be more useful
in understanding the impact of a treatment at the individual level than the
standard average treatment effect; please see appendix S.2 for more details.

Therefore, prediction under covariate shift is a building block of several
important prediction problems and an improvement upon the weighted conformal
method of Tibshirani et al. (2019) would lead to advancements in many other
directions as well, which is exactly what we aim to do in this work. We
reconsider the problem of prediction under covariate shift from a missing data
point of view and provide a novel solution that is more robust and
computationally efficient using modern semiparametric efficiency theory.

##  2 Conformal prediction: literature review

There are different forms of guarantees that one might consider for the
validity of a prediction set. A test point $(X,Y)$ is covered by a set $C$ if
and only if $\mathbbm{1}\\{Y\in C(X)\\}=1$. Often in practice, the set $C$ is
constructed based on a training data and the test point $(X,Y)$ is independent
of $C$. If the prediction set denoted as $\widehat{C}$ is computed from a
training data and the test point $(X,Y)$ is drawn from a distribution $P$,
then the (prediction) miscoverage loss with respect to $P$ is given by

$L_{P}(\widehat{C}):=\mathbb{E}_{(X,Y)\sim
P}[\mathbbm{1}\\{Y\notin\widehat{C}(X)\\}\big{|}\widehat{C}]=\int\mathbbm{1}\\{y\notin\widehat{C}(x)\\}dP(x,y).$

Observe that $L_{P}(\widehat{C})$ is a random variable if $\widehat{C}$ is a
random set, which often it is. Further, note that there is no requirement that
the training data used to construct $\widehat{C}$ comes from the same
distribution $P$ as the test data $(X,Y)$ in defining the loss
$L_{P}(\widehat{C})$.

By a $(1-\alpha)$-prediction set, one might expect/ask for
$L_{P}(\widehat{C})\leq\alpha$. But this is too much to expect/ask in general
without very restrictive assumptions on $P$ or $\widehat{C}$. More actionable
goals are to require $L_{P}(\widehat{C})$ to be less than $\alpha$ either in
expectation or with some specific probability with respect to the randomness
of $\widehat{C}$. We refer to the goal of
$\mathbb{E}[L_{P}(\widehat{C})]\leq\alpha$ as joint coverage. This is the
problem solved by conformal prediction (as introduced by Vovk et al. (2005))
and many of its variants in the literature. The goal of
$\mathbb{P}(L_{P}(\widehat{C})\leq\alpha)\geq 1-\delta$ is referred to as
$(\alpha,\delta)$ probably approximately correct (PAC) in the literature,
which belongs to the class of classical tolerance regions — sets that cover a
pre-specified fraction of the population distribution, see e.g. Guttman (1967)
and Krishnamoorthy and Mathew (2009). Notably, the PAC prediction set depends
on an additional parameter $\delta$ which is not required by the conformal
approach. We note that, in general, neither will imply the other. Vovk (2012)
shows that the split conformal prediction sets with $(1-\alpha)$ joint
coverage also satisfy $(\alpha+f_{N}(\delta),\delta)$-PAC guarantee, for an
explicitly computable $f_{N}(\delta)$.

###  2.1 Under Exchangeability

We will now provide a brief introduction to the literature on conformal
prediction with some emphasis on split conformal prediction that will be
important to understanding the current work. Vovk et al. (2005) first
introduced a version of conformal prediction called the transductive conformal
method (later described in Lei et al. (2013) as the full conformal method)
that requires fitting the learning algorithm to samples $(X_{i},Y_{i}),1\leq
i\leq N$ and $X_{N+1}=x$ for all $x\in\mathcal{X}$. While this method makes
full use of the data for prediction, it is computationally intensive in
practice. Papadopoulos et al. (2002) proposed an alternative method called the
inductive conformal method (or the split conformal method in Lei and Wasserman
(2014)) which splits the data into two different parts, and the learning
algorithm is trained only on the first part and the prediction set is
constructed using conformal “scores” on the second part of the data. For
concreteness, we describe the split conformal algorithm here in the regression
setting with $(X_{i},Y_{i})\in\mathcal{X}\times\mathbb{R}$. Let
$\mathcal{A}:\mathcal{X}\to\mathbb{R}$ be any prediction algorithm trained on
the first split of the data, i.e., for any $x\in\mathcal{X}$, $\mathcal{A}(x)$
is the point prediction for $Y$. If $m$ is the number of observations in the
second split, $R_{i}=|Y_{i}-\mathcal{A}(X_{i})|,i=1,\dots,m$ are residuals
computed on the second split of the data, and $\widehat{Q}_{\alpha}$ is the
$\lceil(m+1)(1-\alpha)\rceil$-th largest element of $R_{i}$’s, then for any
$(X_{N+1},Y_{N+1})$ that is exchangeable with $(X_{i},Y_{i}),i=1,\dots,N$, it
holds that
$\mathbb{P}(|Y_{N+1}-\mathcal{A}(X_{N+1})|\leq\widehat{Q}_{\alpha})\geq
1-\alpha.$ In particular, for
$\widehat{C}_{\alpha}(x)=\\{y:\,|y-\mathcal{A}(x)|\leq\widehat{Q}_{\alpha}\\}$,
we have $\mathbb{P}(Y_{N+1}\in\widehat{C}_{\alpha}(X_{N+1}))\geq 1-\alpha$.
Both full and split conformal algorithms require only the assumption of
exchangeable data for coverage validity. To overcome potential statistical
inefficiency due to data splitting, several papers including Barber et al.
(2021), Romano et al. (2020) and Kim et al. (2020) introduced aggregation
techniques such as jackknife+, CV+, bootstrap after jackknife+ to make better
use of the data. All aforementioned prediction set constructions are valid
under exchangeability of the training data as well as the test point $(X,Y)$
to be covered, as discussed in Kuchibhotla (2020). Also, see Solari and
Djordjilović (2022).

The above works consider joint coverage as the criterion for valid prediction
set and focus on exchangeability for coverage validity. There also is a
separate line of research including, but not limited to, Györfi and Walk
(2020) and Yang and Kuchibhotla (2021) on prediction sets based on i.i.d. data
assumption and concentration inequalities. The advantage is that one attains
coverage guarantees conditional on the training data used to construct the
set, which may be more informative and automatically implies unconditional
(PAC) coverage. To elaborate on this line of thought and as an initial glimpse
into the current work, consider the split conformal prediction procedure
described above. In that description, under i.i.d. data assumption, the
residuals computed on the second split of the data
$R_{i}=|Y_{i}-\mathcal{A}(X_{i})|,i=1,\dots,m$ are i.i.d. along with
$R=|Y-\mathcal{A}(X)|$, irrespective of what $\mathcal{A}(\cdot)$ is. We aim
to find $\widehat{r}_{\alpha}$ such that
$\mathbb{P}(|Y-\mathcal{A}(X)|\leq\widehat{r}_{\alpha}|\mathcal{A})\geq
1-\alpha$. The concentration inequalities approach to prediction sets can now
proceed as follows: the well-known DKW inequality (see e.g. Massart (1990))
implies that

$\mathbb{E}\biggl{[}\sup_{t}\Bigl{|}\frac{1}{m}\sum_{i=1}^{m}\mathbbm{1}\\{|Y_{i}-\mathcal{A}(X_{i})|\leq
t\\}-\mathbb{P}\bigl{(}|Y-\mathcal{A}(X)|\leq
t|\mathcal{A}\bigr{)}\Bigr{|}|\mathcal{A}\biggr{]}\leq\sqrt{\frac{2}{m}}.$

Hence, an $\widehat{r}_{\alpha}$ at which the empirical CDF of the residuals
is at least $(1-\alpha)+\sqrt{2/m}$ yields a valid prediction set of joint
coverage of at least $1-\alpha$. This also implies that if we just use the
sample $(1-\alpha)-$quantile of the residuals, then this would give a
prediction set that has an approximate coverage of $1-\alpha$ with a slack at
most $\sqrt{2/m}.$ Moreover, the application of DKW-type inequalities here
also readily yields a PAC guarantee. The well-developed theory of efficient
estimation in semiparametric theory shows that the sample quantile is an
optimal estimator of the population quantile (Pfanzagl, 1976). This optimality
theory is leveraged throughout the manuscript in deriving well-calibrated
prediction sets under covariate shift.

###  2.2 Under Non-exchangeability

Several authors have considered the problem of prediction sets for non-
exchangeable data. The literature can be divided into three parts: (1) works
that consider independent data but with potential non-identical distributions,
(2) works that consider prediction set construction for dependent data
including time series, and (3) works that are agnostic to the randomness
structure in the data. Our current work belongs to the first category. In the
following, we only review the works directly related to our work and leave the
review of the latter two categories to Section S.1 of the apppendix.

Tibshirani et al. (2019) introduced the problem of prediction set construction
under covariate shift. Building on Tibshirani et al. (2019), Lei and Candès
(2020) worked in a covariate shift setting, establishing asymptotic joint
coverage under certain conditions, accounting for estimation of the covariate
density ratio. Cauchois et al. (2020) provided a method that produces
approximate valid prediction set for any test distribution in an
$f$-divergence ball around the training population and also discussed how to
estimate the expected data shift and build robustness to it. Lei and Candès
(2020) and Kivaranovic et al. (2020) constructed prediction sets for
counterfactuals and individual treatment effects (ITE). The former work uses
the classical SUTVA assumption and rewrites the problem in terms of covariate
shift assumption. The latter work assumes covariates are independent of
treatment assignment and therefore implicitly rules out any covariate shift.
These works represent the first connection of conformal prediction to causal
inference.

In this paper, we will focus on the specific form of non-exchangeability
called _covariate shift_ considered in Tibshirani et al. (2019) and Lei and
Candès (2020) where the training data $\mathcal{D}^{\mbox{tr}}$ is composed of
two parts $\mathcal{D}^{\mbox{tr}}_{P}$ and $\mathcal{D}^{\mbox{tr}}_{Q}$, and
random variables $(X_{i},Y_{i})$’s in $\mathcal{D}^{\mbox{tr}}_{P}$ are i.i.d.
from $P_{X}\otimes P_{Y|X}$, while random variables $X_{i}$’s in
$\mathcal{D}^{\mbox{tr}}_{Q}$ are i.i.d. from $Q_{X}$ with missing response.
The goal is to obtain a prediction set such that the probability a new data
pair $(X_{f},Y_{f})$ falls into this prediction set is larger than some
nominal level $1-\alpha$ where $(X_{f},Y_{f})$ is from $Q_{X}\otimes P_{Y|X}$.
As noted, all prior works focused primarily on achieving this goal
asymptotically with the exception of Tibshirani et al. (2019) who assumed the
covariate shift is known (i.e., known $dQ_{X}/dP_{X}$), in which case they
achieved finite sample coverage guarantee. We prove in this paper that it is
impossible to construct a _non-trivial_ finite sample valid prediction set
without complete a priori knowledge of either the covariate shift or the
conditional distribution of $Y_{f}$ given $X_{f}$; see Theorem 1. We thus
resort to the goal of constructing a prediction set with an asymptotic
coverage of at least $1-\alpha$.

Our approach to construction of prediction sets under covariate shift requires
estimation of the quantile of a univariate function of $(X_{f},Y_{f})$ akin to
the conformal score. This task involves nuisance parameters that must be
estimated sufficiently well to ensure the asymptotic coverage guarantee of our
prediction sets. Fortunately, as we establish using modern semiparametric
theory, this task can be accomplished in a robust fashion by using the
efficient influence function as an estimating equation for the quantile. As we
show, the efficient influence function is endowed with a double robustness
property (see e.g. Scharfstein et al. (1999), Robins (2000) and Bang and
Robins (2005)) which ensures that the coverage bias of our prediction sets can
be made negligible even if the nuisance parameters are estimated at
nonparametric rates using highly adaptive machine learning algorithms.

Before completing the review of the relevant literature, we mention a
concurrent work Qiu et al. (2022) that uses similar connection to
semiparametric statistics as ours does but targets a different type of
asymptotic PAC guarantee than we are able to provide in Theorem 4, a detailed
discussion comparing the two goals and corresponding methods is given at the
end of Section 3.1.

#### Organization

The remainder of the paper is organized as follows. In Section 3, we formally
introduce the explainable covariate shift problem, providing the missing data
formulations of the problem that gives notation used throughout the paper and
also makes its connection to causal inference. In Section 4, we formally
establish that it is impossible to construct a well-calibrated prediction set
that is informative without a priori knowledge about the covariate shift such
as knowledge of the ratio of covariate densities, in the sense that any valid
prediction set would with high probability have infinite Lebesgue measure. In
Section 5, we provide our first doubly robust algorithm for the explainable
covariate shift problem, using a sample splitting strategy which attains
nominal asymptotic prediction coverage. In Section 6, we provide a second
doubly robust algorithm that makes more efficient use of the observed data by
avoiding sample splitting, with guaranteed validity under certain regularity
conditions we establish. Section 7 reports simulation studies (both synthetic
and real data) validating the theoretical results of the proposed doubly
robust methods, and comparing them to the weighted conformal prediction method
of Tibshirani et al. (2019). In Section 8, we consider the problem of
aggregating a collection of prediction sets by providing an explicit algorithm
adapted from Yang and Kuchibhotla (2021) in an effort to optimize prediction
accuracy; we illustrate the efficiency of the proposed algorithm through
simulations. In Section 9, we relax the explainable covariate shift assumption
by allowing for the presence of latent covariate shift encoded in a
sensitivity parameter and discuss the efficient influence function. This is
the first step in extending the framework of the current paper to account for
a departure from explainable covariate shift as well as constructing
prediction sets for ITE under unmeasured confounding. Further elaboration of
Section 9 will appear elsewhere. Finally in Section 10, we conclude the paper
with a brief discussion.

Proofs of all results and supporting lemmas are provided in the supplementary
where for convenience, the sections and equations are prefixed with “S.” and
“E.”, respectively.

##  3 Our Problem and Notation

In this section, we provide a formal description of the _explainable_
covariate shift problem, and introduce both missing data and counterfactual
formulations of the problem.

###  3.1 The Covariate Shift Problem

The common assumption that the training and test data follow a common
probability distribution can fail in practice, for example training data may
be collected under stringent laboratory conditions that cannot be met when
deployed in clinical practice; likewise, image training data may be obtained
in one region whereas the test data may be collected in another. Therefore, it
is important to consider situations where training and test distributions are
different, also known as covariate shift, which we formalize below. Assuming
we have training data $\mathcal{D}^{\mbox{tr}}$ composed of two parts
$\mathcal{D}^{\mbox{tr}}_{P}$ and $\mathcal{D}^{\mbox{tr}}_{Q}$, where

$\mathcal{D}^{\mbox{tr}}_{P}~{}:=~{}\\{Z_{i}=(X_{i},Y_{i}):\,1\leq i\leq
n\\}\quad\text{ and
}\quad\mathcal{D}^{\mbox{tr}}_{Q}~{}:=~{}\\{Z_{i}=X_{i}:\,n+1\leq i\leq N\\},$
(2)

and random variables in $\mathcal{D}^{\mbox{tr}}_{P}$ are i.i.d. from
$P_{X}\otimes P_{Y|X}$, while random variables in
$\mathcal{D}^{\mbox{tr}}_{Q}$ are i.i.d. from $Q_{X}$. Note that in
$\mathcal{D}^{\mbox{tr}}_{Q}$ only data on covariates are available and the
outcome/response is missing, and thus the subject of prediction. In such
setting, a covariate shift problem is said to be present as the covariates in
$\mathcal{D}^{\mbox{tr}}_{Q}$ are sampled from $Q_{X}$ which may be different
from $P_{X}$ (the distribution of covariates in
$\mathcal{D}^{\mbox{tr}}_{P}$). This setting readily extends to the more
general case where one also observes samples from $Q_{X}\otimes P_{Y|X}$,
and/or samples from $Q_{X}\otimes Q_{Y|X}$; settings that are closely related
to transfer learning (see e.g. Kpotufe and Martinet (2018) and Reeve et al.
(2021)) and in the special case when $P_{X}$ and $Q_{X}$ are identical, this
reduces exactly to the setting of semi-supervised learning (see e.g. Zhu and
Goldberg (2009) and Zhang et al. (2019)). Though several variants of the
covariate shift problem have been of interest in statistics and ML
literatures, the specific setup considered has recently generated renewed
interest in ML literature. While the focus has usually been on approaches to
account for covariate shift while conducting model selection such as
regression, see e.g. Sugiyama et al. (2007), Quiñonero-Candela et al. (2008),
Bickel et al. (2009), Reddi et al. (2015), Chen et al. (2016), the goal here
is to obtain a prediction set such that

$\mathbb{P}\bigl{(}Y_{f}\in\widehat{C}_{N,\alpha}(X_{f})\bigr{)}\geq
1-\alpha,\quad\mbox{whenever}\quad(X_{f},Y_{f})\sim Q_{X}\otimes P_{Y|X}.$ (3)

This problem was first posed in Tibshirani et al. (2019) where they developed
a weighted version of conformal prediction that produces valid prediction sets
when the likelihood ratio between training and test distributions is known.
The idea of their construction is similar to importance sampling Monte Carlo.
If, as in most practical settings, the likelihood ratio between the two
distributions is unknown and therefore must be estimated, they empirically
demonstrate in a low-dimensional setting that approximate coverage might still
be possible via simulation studies. However, Tibshirani et al. (2019) do not
formally consider the extent to which bias in estimating the likelihood ratio
propagates to impact coverage. In addition, as noted by the authors of that
paper, for every new test point $(x_{f},y_{f})$, their prediction set that
involves a weighted quantile has to be recalculated, and this would be
computationally intensive. Specifically, their approach requires that the test
point $x_{f}$ is specified in advance.

In this work we will construct prediction regions $\widehat{C}_{N}$ that are
determined by one-dimensional functions of $(X_{f},Y_{f})$, i.e., we take an
arbitrary function $(x,y)\mapsto R(x,y)\in\mathbb{R}$ and estimate the
quantile of $R(X_{f},Y_{f})$. For now, we will think of $R(\cdot,\cdot)$ as a
fixed non-stochastic function. In practice, $R(\cdot,\cdot)$ is a conformal
score computed from an independent sample; examples include
$R(x,y)=|y-\mathcal{A}(x)|$ (regression residual, Lei et al., 2018) or
$R(x,y)=\max\\{\widehat{q}_{\alpha/2}(x)-y,y-\widehat{q}_{1-\alpha/2}(x)\\}$
(conformalized quantile residual with estimated conditional quantiles
$\widehat{q}_{\alpha/2}(\cdot)$ and $\widehat{q}_{1-\alpha/2}(\cdot)$, Romano
et al., 2019). If $r_{\alpha}$ is the smallest $(1-\alpha)$-quantile of
$R(X_{f},Y_{f})$ in the target population $Q_{X}\otimes P_{Y|X}$, then

$\mathbb{P}\bigl{(}R(X_{f},Y_{f})\leq r_{\alpha}\bigr{)}\geq 1-\alpha,$

and hence, for $C_{\alpha}=\\{(x,y):\,R(x,y)\leq r_{\alpha}\\}$, we have
$\mathbb{P}\bigl{(}(X_{f},Y_{f})\in C_{\alpha}\bigr{)}\geq 1-\alpha$. This
result holds, irrespective of the choice of function $R(\cdot,\cdot)$. Note
that the assumption that $Y_{i}$ conditional on $X_{i}$ for $1\leq i\leq n$
has the same distribution as $Y_{f}$ conditional on $X_{f}$ essentially means
that the covariate shift problem can be completely accounted for by
conditioning on observed covariates $X$, hence the reference to this setting
as “explainable covariate shift” problem. It implies that the conditional
distribution of $R(X_{i},Y_{i})$ given $X_{i}$ coincides with the conditional
distribution of $R(X_{f},Y_{f})$ given $X_{f}$, i.e.,

$\mathbb{P}\bigl{(}R(X_{i},Y_{i})\in
B\big{|}X_{i}=x\bigr{)}=\mathbb{P}\bigl{(}R(X_{f},Y_{f})\in
B\big{|}X_{f}=x\bigr{)},1\leq i\leq n$ (4)

for all Borel sets $B\subseteq\mathbb{R}$. We borrow results from
semiparametric theory and estimate $r_{\alpha}$ based on its efficient
influence function, which is intimately related to efficient influence
function of the average treatment effect among the treated (ATT) functional;
we then combine this influence function with arbitrary training map
$R(\cdot,\cdot)$ and establish that the resulting prediction set has
asymptotic nominal coverage with coverage bias of a product form which implies
correct coverage if either the likelihood ratio between the two distributions
or the conditional distribution of $R$ given $X$ can be estimated sufficiently
well, also known as double robustness. Note that in our construction of the
prediction set $\widehat{C}_{N,\alpha}$, the mapping $R$ does not depend on
the test point $x_{f}$ at which prediction is needed. Tibshirani et al. (2019)
approached (3) using the equation

$\mathbb{P}_{(X_{f},Y_{f})\sim Q_{X}\otimes
P_{Y|X}}(R(X_{f},Y_{f})\leq\theta)=\mathbb{E}_{(X,Y)\sim P_{X}\otimes
P_{Y|X}}\left[\mathbbm{1}\\{R(X,Y)\leq\theta\\}\frac{dQ_{X}}{dP_{X}}(X)\right].$
(5)

This implies that consistent estimation of $dQ_{X}/dP_{X}$ allows for
consistent estimation of $r_{\alpha}$. But, given (4), consistent estimation
of the conditional distribution of $R(X,Y)$ given $X$ also allows for
consistent estimation of $r_{\alpha}$. This gives a hint at double robustness
in estimating $r_{\alpha}$ which we will formalize later in the paper.

In Lei and Candès (2020), the authors proposed a method that targets the
covariate shift problem under the framework of counterfactual prediction in a
causal inference setting, which we consider in Section S.2. Interestingly, the
coverage of their prediction set has bias of the order of the minimum of two
errors, that of the prediction ML algorithm and that of the estimated
covariate likelihood ratio, a property that appears to hold under the so-
called conformal quantile regression (CQR) which restricts the choice of
conformal score to a quantile regression function for the outcome in view; it
is unclear whether similar robustness extends beyond CQR. In contrast, while
our approach equally applies to the counterfactual prediction framework
considered by Lei and Candès (2020), as we establish, its coverage bias is
guaranteed to be of the order of the product of two errors, that of an
estimated quantile function for $R$ with that of the covariate likelihood
ratio, an immediate consequence of double robustness. Therefore, the bias of
our coverage error rate can be substantially smaller relative to that of Lei
and Candès (2020). In addition, the product bias property of the proposed
method is guaranteed to hold for any ML technique used to empirically
construct $R$, therefore making our approach potentially more general than
theirs. Park et al. (2021) uses probably approximately correct (PAC)
prediction sets (tolerance regions that cover a pre-specified fraction of the
population distribution) for deep learning models including in the presence of
covariate shift that also requires prior knowledge on the shift of the
distributions. Notably, the PAC prediction set depends on an additional
parameter $\delta$ which is not required by our approach.

Concurrently to our paper, Qiu et al. (2022) consider a related prediction
setting, also drawing from modern semiparametric theory, but mainly focusing
on a particular form of asymptotic PAC guarantee with covariate shift. We view
these contributions as complementary. Specifically, in our covariate shift
setting, asymptotic joint coverage is interpreted as:
$\mathbb{P}(Y_{f}\in\widehat{C}_{N,\alpha}(X_{f}))\geq 1-\alpha-o(1)$ as
$N\to\infty$. In contrast, the $(\alpha,\delta)$-PAC guarantee can be met
approximately with negligible errors in either $\alpha$, $\delta$ or both. The
target of Qiu et al. (2022) is to develop $(\alpha,\delta+o(1))$-PAC
prediction set. The methodology developed in the current paper provides
$(\alpha+o_{p}(1),\delta)$-PAC guarantee along with the joint coverage
guarantee. Whichever version of asymptotic PAC guarantee is more useful
depends on the application. Notably, in order to minimize the impact of bias
due to nuisance parameter estimation, Qiu et al. (2022) leverage the product
bias property of the efficient influence function for the coverage probability
of their prediction set, while we leverage the product bias structure of the
efficient influence function for the $(1-\alpha)$-quantile of the statistic
generating the prediction set; these two influence functions are equal up to a
multiplicative constant for a given coverage guarantee, although the
difference in their use leads to a different guarantee. Importantly, their PAC
approach involves the construction of a valid confidence interval for the
coverage probability which in turn requires a regular asymptotic linear
estimator of coverage, in which case, their product bias must be of order
smaller than root-$n$. Our proposed approach does not have this requirement,
and therefore our guarantee is attainable even if the product bias is of order
larger than root-$n$ provided that at least one of our nuisance functions is
consistent.

###  3.2 Reformulation as a Missing Data Problem

In this section, we formulate the covariate shift problem in a missing data
framework, and introduce notation used throughout the remainder of the paper.
This reformulation allows us to make use of the modern theory of
semiparametric statistics. Recall the training data (2). For each
$(X_{i},Y_{i})$ contained in $\mathcal{D}^{\mbox{tr}}_{P}$, define
$R_{i}=R(X_{i},Y_{i})$ and set $T_{i}=0$. For each $(X_{i},Y_{i})$ contained
in $\mathcal{D}_{Q}^{\mbox{tr}}$, $R_{i}$ is unobserved because the
corresponding $Y_{i}$ is unobserved. Hence the observed data
$\mathcal{D}=\mathcal{D}_{P}^{\mbox{tr}}\cup\mathcal{D}_{Q}^{\mbox{tr}}$ can
be succinctly written as
$Z_{i}=\bigl{(}X_{i},T_{i},(1-T_{i})R_{i}\bigr{)},1\leq i\leq N$ such that

$\mathbb{P}(X_{i}\in
A|T_{i}=0)=:P_{X}(A)\quad\mbox{and}\quad\mathbb{P}(X_{i}\in
A|T_{i}=1)=:Q_{X}(A),$ (6)

while

$\mathbb{P}(R_{i}\in B|T_{i}=0,X_{i}=x)=\mathbb{P}(R_{i}\in
B|T_{i}=1,X_{i}=x)=:P_{Y|X=x}(B)\quad\mbox{a.e.}\quad x.$ (7)

Equations in (6) signify that the covariates in $\mathcal{D}_{P}^{\mbox{tr}}$
are distributed as $P_{X}$ and that the covariates in
$\mathcal{D}_{Q}^{\mbox{tr}}$ are distributed as $Q_{X}$. Condition (7), on
the other hand, signifies that the conditional distribution of the response
$Y$ given $X$ is the same for $\mathcal{D}_{P}^{\mbox{tr}}$,
$\mathcal{D}_{Q}^{\mbox{tr}}$, and the future data $(X_{f},Y_{f})$. Condition
(7) restates condition (4) in terms of $T_{i}$.

Condition (7) implies that $R_{i}$ is independent of $T_{i}$ conditional on
$X_{i}$ for all $1\leq i\leq N$; this is denoted by $R_{i}\perp T_{i}|X_{i}$
and is equivalent to the missing at random (MAR) assumption in missing data
literature. This assumption, which is not testable without an additional
condition, essentially states that there is no unmeasured factor that is
related with both $R$ and $T$. For identification, we further assume that for
any Borel set $B$,

$P_{X}(B)=0\text{ implies }Q_{X}(B)=0.$ (8)

This is the same as assuming the measure $Q_{X}$ is absolutely continuous
w.r.t $P_{X}$. In other words, the support of $X|T=1$ is contained in the
support of $X|T=0$. Assumption (8) is needed for (5), which is a crucial
component of the double robustness of our methodology.

Summarizing the above discussion, the covariate shift assumption which states
that only the covariate distributions between $\mathcal{D}_{P}^{\mbox{tr}}$
and $\mathcal{D}_{Q}^{\mbox{tr}}$ can be different but not the conditional
distributions of the response given the covariates is equivalent to the MAR
assumption.

For any $\theta\in\mathbb{R}$, define sets

$\widehat{C}(\theta;x):=\bigl{\\{}y:R(x,y)\leq\theta\bigr{\\}}.$ (9)

In this notation, under condition (8), we aim to find a data-dependent random
variable $\widehat{r}_{\alpha}$ such that for any $(X_{f},Y_{f})\sim
Q_{X}\otimes P_{Y|X}$,

$\displaystyle\mathbb{P}\bigl{(}Y_{f}\in\widehat{C}(\widehat{r}_{\alpha};X_{f})\bigr{)}=\mathbb{P}_{(X,Y)\sim
Q_{X}\times P_{Y\mid
X}}(R(X,Y)\leq\widehat{r}_{\alpha})=\mathbb{P}(R(X,Y)\leq\widehat{r}_{\alpha}|T=1)\geq
1-\alpha,$ (10)

while $r_{\alpha}$, the “target” of $\widehat{r}_{\alpha}$ is defined to be
the smallest real number such that

$\mathbb{P}\bigl{(}R(X,Y)\leq r_{\alpha}|T=1\bigr{)}\geq 1-\alpha.$ (11)

Because we do not observe random variables $R$ when $T=1$, this goal is not
achievable in finite samples without restrictive assumptions such as a known
$\mathbb{P}(T=1|X)/\mathbb{P}(T=0|X)$ (Tibshirani et al., 2019). We provide a
random variable $\widehat{r}_{\alpha}$ so that (10) is achieved with a slack
that converges to zero as $N$ tends to $\infty.$

##  4 Impossibility of finite sample coverage

Recall our aim from (10). Resorting to semiparametric theory, in most cases,
implies that the resulting coverage guarantee is only asymptotic. In our
problem of covariate shift, one can prove that it is impossible to construct a
finite sample valid _non-trivial_ prediction set without the knowledge of
either the covariate shift or the conditional distribution of $Y$ given $X$.
Here, by a non-trivial prediction set, we mean a set with a finite Lebesgue
measure. Lemma S1 of Section S6.1 of Qiu et al. (2022) prove an analogous but
weaker result for PAC guarantee as they establish a result similar to (13) of
the following Theorem, however (14) appears to be an entirely novel
contribution.

######  Theorem 1.

Suppose the observed data consists of $n$ i.i.d. tuples
$(X_{i},T_{i},(1-T_{i})Y_{i})$. Further assume that
$\mathcal{X}\subseteq\mathbb{R}^{d}$ is the support of $X_{i}$ and
$\mathcal{Y}\subseteq\mathbb{R}$ is the support of $Y_{i}$. Let
$\bar{\mathcal{P}}^{0}$ be the set of all distributions $\bar{P}^{0}$ on the
random vector $\bar{O}=(X,T,Y)$ such that $T$ is independent of $Y$ given $X$
($T\perp Y|X$), and the joint distribution of $(X,Y)$ is absolutely continuous
with respect to the Lebesgue measure on $\mathcal{X}\times\mathcal{Y}$.

Suppose that a (possibly randomized) prediction set $\widehat{C}_{\alpha}$ has
finite-sample joint coverage guarantee in the target population, that is,

$\sup_{\bar{P}^{0}\in\bar{\mathcal{P}}^{0}}\,\mathbb{P}_{\bar{P}^{0}}(Y\notin\widehat{C}_{\alpha}(X)\mid
T=1)\leq\alpha,\quad\mbox{for some}\quad\alpha\in(0,1).$ (12)

Then, for any $\bar{P}^{0}\in\bar{\mathcal{P}}^{0}$ and a.e. $y\in\mathcal{Y}$
with respect to the Lebesgue measure,

$\mathbb{P}_{\bar{P}^{0}}\bigl{(}y\notin\widehat{C}_{\alpha}(X)\bigr{)}\leq\alpha.$
(13)

Furthermore, $\widehat{C}_{\alpha}(X)$ would at least cover one of the end
points $\mathcal{Y}$ with probability at least $1-\alpha$, and hence if
$\mathcal{Y}=\mathbb{R}$, then

$\mathbb{E}_{\bar{P}^{0}}[\mbox{Leb}(\widehat{C}_{\alpha}(X))]=\infty.$ (14)

We defer the proof to the Section S.8 of the appendix. The proof is based on
the lack of a non-trivial test for the problem of conditional independence
hypothesis testing proved in Shah and Peters (2020). The connection to
conditional independence testing can be seen from the fact that any prediction
set that is valid under the conditional independence of $T$ and $Y$ given $X$
can be rewritten as a valid test for the hypothesis that $T$ is conditionally
independent of $Y$ given $X$; see the proof in Section S.8 for more details.

Given the lack of finite-sample valid non-trivial prediction, we resort to
finding an efficient asymptotically valid prediction set based on semi-
parametric theory in the following sections.

##  5 Methodology with split data

In this section, we discuss our novel prediction set construction under the
covariate shift setting using semiparametric theory.

Recall our notation $Z=(X,T,(1-T)R)$. Suppose that one is interested in the
$q$th-quantile of a random variable $R|T=1$, denoted
$\theta_{0}:=\inf\\{r:F(r)\geq q\\}$, where $F$ is the CDF of $R|T=1$. An
estimator $\widehat{\theta}$ is said to be asymptotically linear if it
satisfies

$\sqrt{N}\bigl{(}\widehat{\theta}-\theta_{0}\bigr{)}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\psi(Z_{i})+o_{p}(1),\quad\mathbb{E}[\psi(Z)]=0,\mathbb{E}\bigl{[}\psi(Z)^{\top}\psi(Z)\bigr{]}<\infty,$

as $N\to\infty$. Clearly, the asymptotic variance of $\widehat{\theta}$ is
then $\mathbb{E}\bigl{[}\psi(Z)\psi(Z)^{\top}\bigr{]}/N$. The function
$\psi(z)$ is referred to as the influence function, following terminology of
Hampel (1974). Furthermore, as the model is nonparametric in the sense that
the observed data distribution is not restricted and all such distributions
are regular, (see Chapter 2 of Bickel et al. (1993) on regularity), any
estimator satisfying the above expansion is said to attain the semiparametric
efficiency bound, and $\psi(\cdot)$ is said to be the efficient influence
function of $\theta_{0}$ in the nonparametric model. For more on influence
functions and semiparametric theory, see e.g. Newey (1990), Chapter 25 of Van
der Vaart (2000) and van der Vaart (2002).

We now state the efficient influence function for the $(1-\alpha)$-th quantile
of $R|T=1$ under the MAR assumption. For every
$x\in{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}$
and $r\in\mathbb{R}$, define

$\begin{split}\pi^{\star}(x)~{}&:=~{}\mathbb{P}(T=1|X=x)/\mathbb{P}(T=0|X=x),\\\
m^{\star}(r,x)~{}&:=~{}\mathbb{E}[\mathbbm{1}\\{R\leq r\\}|X=x].\end{split}$
(15)

The function $\pi^{\star}(\cdot)$ represents the true density ratio of the
covariates among labeled and unlabeled data. The function
$m^{\star}(\cdot,\cdot)$ represents the true conditional mean function that by
assumption (4) is common for both labeled and unlabeled data.

We now state the efficient influence function for the quantile of interest
under regularity conditions for the data distribution that will motivate our
proposed method, where it should be noted that strictly speaking, the
regularity conditions are actually not needed for the theoretical guarantees
of our proposed methods.

######  Lemma 1.

Suppose $\mathbb{E}[\mathbbm{1}\\{T=0\\}\pi^{\star
2}(X)]=\mathbb{E}{\mathbb{P}^{2}(T=1|X)}/{\mathbb{P}(T=0|X)}$ is finite and
that the density of the conditional distribution of $R|T=1$ at $r_{\alpha}$ is
bounded away from zero. Then the efficient influence function of the
$(1-\alpha)$-quantile of $R|T=1$ in the nonparametric model for $Z$ which
allows the distribution of $Z$ to remain unrestricted under the condition that
it is regular, is given up to a proportionality constant by

$\displaystyle\psi(z)=\mathrm{IF}(r_{\alpha},x,r,t;\pi^{\star},m^{\star})$
$\displaystyle=\mathbbm{1}\\{t=0\\}\pi^{\star}(x)\Big{[}\mathbbm{1}\\{r\leq
r_{\alpha}\\}-m^{\star}(r_{\alpha},x)\Big{]}$ (16)
$\displaystyle\quad+\mathbbm{1}\\{t=1\\}\Big{[}m^{\star}(r_{\alpha},x)-(1-\alpha)\Big{]}.$

###### Proof.

Hahn (1998) gave a derivation of the influence function for the average
treatment effect among the treated (ATT). We adapt their proof to that of the
conditional quantile and give the complete derivation of (16) along with a
basic introduction to semiparametric theory needed to derive the result in
Section S.3 of the supplementary. ∎

For any two functions $\pi(\cdot)$ and $m(\cdot,\cdot)$, let

$\displaystyle\mathrm{IF}(\theta,x,r,t;\pi,m):=\mathbbm{1}\\{t=0\\}\pi(x)\Big{[}\mathbbm{1}\\{r\leq\theta\\}-m(\theta,x)\Big{]}+\mathbbm{1}\\{t=1\\}\Big{[}m(\theta,x)-(1-\alpha)\Big{]},$
(17)

ignoring the scaling factor. Note that $\mathrm{IF}(\theta,x,r,t;\pi,m)$ is
only a function of $(x,t,(1-t)r)$ because the term that depends on $r$ has a
multiplicative factor of $\mathbbm{1}\\{t=0\\}$. Let $P[f]$ denote integration
conditional on the training sample. For example, for any $\theta,\pi,m$ that
are potentially data-dependent,

$P[\mathrm{IF}(\theta,x,r,t;\pi,m)]=\int\mathrm{IF}(\theta,x,r,t;\pi,m)dP_{R|X=x}(r|x)dP_{T|X=x}(t|x)dP_{X}(x).$

This is a random variable if $\theta$ or $\pi$ or $m$ are random.

First, we draw a key connection between the desired coverage and the
aforementioned influence function.

######  Lemma 2.

Let
$\pi:{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}\to\mathbb{R}_{+}$
and
$m:\mathbb{R}\times{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}\to[0,1]$
be any two functions. Then for every (potentially) data-dependent
$\theta\in\mathbb{R}$, the representation

$\displaystyle\mathbb{P}_{(X,Y)\sim Q_{X}\otimes
P_{Y|X}}\bigl{(}Y\in\widehat{C}(\theta;X)\mid\theta\bigr{)}~{}$
$\displaystyle=~{}\mathbb{P}(Y\in\widehat{C}(\theta;X)|\theta,T=1)$ (18)
$\displaystyle=~{}1-\alpha+\frac{P[\mathrm{IF}(\theta,X,R,T;\pi,m)]}{\mathbb{P}(T=1)},$
(19)

holds true, whenever either of the following holds true:

  1. 1.

$\pi(x)=\pi^{\star}(x)$ for all $x$; or

  2. 2.

$m(\gamma,x)=m^{\star}(\gamma,x)$ for all $\gamma$ and $x$.

###### Proof.

See Section S.4 for a proof. ∎

Note that even if data-dependent, an estimate of $\theta$ remains independent
of a future observation $(X,Y)$. Lemma 2 has a key implication that
$\mathrm{IF}(\cdots)$ is a doubly robust influence function. Because

$\mathbb{P}\bigl{(}Y\in\widehat{C}(\theta;X)\mid\theta,T=1\bigr{)}=\mathbb{P}\bigl{(}R(X,Y)\leq\theta\mid\theta,T=1\bigr{)},$

taking $\theta=r_{\alpha}$, the quantile of $R(X,Y)$ conditional on
$T=1$,111We assume here that $\mathbb{P}_{(X,Y)\sim Q_{X}\otimes
P_{Y|X}}(R(X,Y)\leq r_{\alpha})=1-\alpha$, which is mild as one can add small
Gaussian noise to $R(X,Y)$. Lemma 2 implies that

$P[\mathrm{IF}(r_{\alpha},X,R,T;\pi,m)]=0,\quad\mbox{if either
}\quad\pi\equiv\pi^{\star}\mbox{ or }m\equiv m^{\star}.$ (20)

Because $r_{\alpha}$ is a constant,
$P[\mathrm{IF}(r_{\alpha},X,R,T;\pi,m)]=\mathbb{E}[\mathrm{IF}(r_{\alpha},X,R,T;\pi,m)]$.

######  Remark.

Note that our results do not actually require uniqueness of a solution to
$\mathbb{E}\big{[}\mathrm{IF}(\theta,x,r,t;\pi^{\star},m^{\star})\bigl{]}=0$.
This could arise for example in settings when $R$ is discrete. In principle,
our result would continue to hold for any element $\theta$ of a solution set.

Property (20) is a double robustness property in that the expectation is zero,
as long as one of $\pi$ and $m$ is the true function. This property implies
that when, as would generally be the case in practice, $\pi^{\star}$ and
$m^{\star}$ are estimated, the resulting bias is of the following product
form, where we introduce the notation $\|f(\cdot)\|_{2}$ as the $L_{2}$-norm
of $f(\cdot)$, where $\|f(\cdot)\|_{2}:=[\int f^{2}(x)M(\mathrm{d}x)]^{1/2}$,
and $M(\cdot)$ is the probability measure of $X$ that for any Borel set $B$,
$M(B)$ is given by

$\displaystyle M(B)$ $\displaystyle=\mathbb{P}(X\in B)$
$\displaystyle=\mathbb{P}(X\in B|T=1)\mathbb{P}(T=1)+\mathbb{P}(X\in
B|T=0)\mathbb{P}(T=0)$
$\displaystyle=Q_{X}(B)\mathbb{P}(T=1)+P_{X}(B)\mathbb{P}(T=0).$

######  Theorem 2.

For any functions $\widehat{\pi}(\cdot)$ and $\widehat{m}(\cdot,\cdot)$, it
holds that

$\displaystyle\sup_{\gamma\in\mathbb{R}}\Bigl{|}P\bigl{[}\mathrm{IF}(\gamma,X,R,T;\widehat{\pi},\widehat{m})-\mathrm{IF}(\gamma,X,R,T;{\pi}^{\star},m^{\star})\bigr{]}\Bigr{|}$
(21)
$\displaystyle\qquad\quad\leq\|\widehat{\pi}-\pi^{\star}\|_{2}\sup_{\gamma}\|\widehat{m}(\gamma,\cdot)-m^{\star}(\gamma,\cdot)\|_{2}.$
(22)

###### Proof.

See Section S.5 for a proof. ∎

Theorem 2 implies that
$P[\mathrm{IF}(\gamma,X,R,T;\widehat{\pi},\widehat{m})]$ converges to
$P[\mathrm{IF}(\gamma,X,R,T;\pi^{\star},m^{\star})]$ as long as one of
$\pi^{\star}$ and $m^{\star}$ is estimated consistently.

Lemma 2 is the main building block for our methodology. In order to ensure
approximately correct coverage of $1-\alpha$, we need to find $\theta$ such
that $P[\mathrm{IF}(\theta,X,R,T;\pi,m)]$ is approximately zero, with either
$\pi\equiv\pi^{\star}$ or $m\equiv m^{\star}$. In practice where we do not
have access to either $\pi^{\star}$ or $m^{\star}$ and even if we know either
of them, one cannot compute $P[\mathrm{IF}(\theta,X,R,T;\pi,m)]$ without
access to the true distribution of $(X,(1-T)R,T)$. Our methodology, hence, is
as follows. We construct estimators $\widehat{\pi}(\cdot)$ and
$\widehat{m}(\cdot,\cdot)$ such that $\|\widehat{\pi}-\pi\|_{2}=o_{p}(1)$ and
$\|\widehat{m}-m\|_{2}=o_{p}(1)$ for some $\pi(\cdot)$ and $m(\cdot,\cdot)$,
and either $\pi\equiv\pi^{\star}$ or $m\equiv m^{\star}$. Then we find the
smallest $\widehat{\theta}$ such that

$\mathbb{P}_{N}\bigl{[}\mathrm{IF}(\widehat{\theta},X,R,T;\widehat{\pi},\widehat{m})\bigr{]}:=\frac{1}{N}\sum_{1\leq
i\leq
N}\mathrm{IF}(\widehat{\theta},X_{i},R_{i},T_{i};\widehat{\pi},\widehat{m})\geq
0.$ (23)

We can prove under certain regularity conditions on $\widehat{\pi}$ and
$\widehat{m}$ that

$\mathbb{P}_{N}\bigl{[}\mathrm{IF}(\widehat{\theta},X,R,T;\widehat{\pi},\widehat{m})\bigr{]}-P\bigl{[}\mathrm{IF}(\widehat{\theta},X,R,T;\pi,m)\bigr{]}=o_{p}(1),$
(24)

even for a data-dependent $\widehat{\theta}$. Then, Lemma 2 implies that
$\widehat{\theta}$ satisfying (23) also satisfies

$\mathbb{P}_{(X,Y)\sim Q_{X}\otimes
P_{Y|X}}\bigl{(}Y\in\widehat{C}(\widehat{\theta};X)\mid\widehat{\theta}\bigr{)}\geq(1-\alpha)+o_{p}(1).$
(25)

This yields the desired coverage guarantee (10). In finding $\widehat{\theta}$
and proving (24), one can avoid restrictive regularity conditions (such as
smoothness or Donsker class) on $\widehat{\pi},\widehat{m}$ by splitting the
data into two parts, using the first part to determine
$\widehat{\pi},\widehat{m}$ and using the second part to compute
$\mathbb{P}_{\mathcal{I}_{2}}[\mathrm{IF}(\cdots)]$. The detailed split sample
procedure is succinctly described in Algorithm 1.

Input: Training data
$\mathcal{D}^{\text{tr}}=\mathcal{D}^{\text{tr}}_{P}\cup\mathcal{D}^{\text{tr}}_{Q}$;
Coverage probability $1-\alpha$, a training method $\mathcal{A}$ and
estimators $\widehat{\pi},\widehat{m}$, the point for prediction $x$.

Output: A valid prediction set $\widehat{C}_{\alpha}(x)$.

1 Split training data $\mathcal{D}^{\text{tr}}$ randomly into
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$, where
$\mathcal{D}_{1}=\\{Z_{i}\in\mathcal{D}^{\text{tr}},i\in\mathcal{I}_{1}\\}$
and
$\mathcal{D}_{2}=\\{Z_{i}\in\mathcal{D}^{\text{tr}},i\in\mathcal{I}_{2}\\}$.;

2 Fit the training method $\mathcal{A}$ on $\mathcal{D}_{1}$ and using fitted
method $\mathcal{A}$, construct an increasing (nested) sequence of sets
$\\{\mathcal{F}_{t}\\}_{t\in\mathcal{T}}$. Here $\mathcal{T}$ is a subset of
$\mathbb{R}$. The nested sets $\\{\mathcal{F}_{t}\\}_{t\in\mathcal{T}}$ can
depend arbitrarily on $\mathcal{D}_{1}$.;

3 For each $i\in\mathcal{I}_{2}$ that satisfies $T_{i}=0$, define the
conformal score
$r_{i}~{}=~{}r(Z_{i}):=\inf\\{t\in\mathcal{T}:\,Z_{i}\in\mathcal{F}_{t}\\}.$ ;

4 Fit estimators $\widehat{\pi},\widehat{m}$ on $\mathcal{D}_{1}$ and find the
smallest $\widehat{\theta}=\widehat{r}_{\alpha}$ such that
$\mathbb{P}_{\mathcal{I}_{2}}[\mathrm{IF}(\widehat{\theta},X,R,T;\widehat{\pi},\widehat{m})]\geq
0$, where
$\displaystyle\mathbb{P}_{\mathcal{I}_{2}}[\mathrm{IF}(\widehat{\theta},X,R,T;\widehat{\pi},\widehat{m})]~{}$
$\displaystyle=~{}\frac{1}{|\mathcal{I}_{2}|}\sum_{i\in\mathcal{I}_{2}}\mathbbm{1}\\{t_{i}=0\\}\widehat{\pi}(x_{i})\Big{[}\mathbbm{1}\\{r_{i}\leq\widehat{\theta}\\}-\widehat{m}(\widehat{\theta},x_{i})\Bigr{]}$
$\displaystyle\quad+~{}\frac{1}{|\mathcal{I}_{2}|}\sum_{i\in\mathcal{I}_{2}}\mathbbm{1}\\{t_{i}=1\\}\Bigl{[}\widehat{m}(\widehat{\theta},x_{i})-(1-\alpha)\Bigr{]}.$

return the prediction set
$\widehat{C}_{\alpha}(x):=\bigl{\\{}y:R(x,y)\leq\widehat{r}_{\alpha}\bigr{\\}}$.

Algorithm 1 Split doubly robust prediction

In step 2 of Algorithm 1, the training method $\mathcal{A}$ can be a
regression estimator of $Y$ on $X$ leading to $\widehat{\mu}$ and the nested
sets, for example, could be
$\mathcal{F}_{t}=\\{(x,y):\,|y-\widehat{\mu}(x)|\leq t,t\geq 0\\}$. This
corresponds to using $R=R(x,y)=|y-\widehat{\mu}(x)|$. Alternatively, one can
also consider a training method that leads to conditional quantile estimator
$\widehat{q}_{\alpha/2}(\cdot)$ and $\widehat{q}_{1-\alpha/2}(\cdot)$. The
nested sets, for example, could be
$\mathcal{F}_{t}=\\{(x,y):\,y\in[\widehat{q}_{\alpha/2}(x)-t,\widehat{q}_{1-\alpha/2}(x)+t]\\}.$
This corresponds to the map
$R=R(x,y)=\max\\{\widehat{q}_{\alpha/2}(x)-y,y-\widehat{q}_{1-\alpha/2}(x)\\}$
which is the conformal score of the conformalized quantile regression (CQR)
method of Romano et al. (2019). For discrete/categorial response $Y$,
conditional probability of each class $\mathbb{P}(Y=j|X=x)$ could be used to
construct the nested sets and the map $R$, see for example, Section 4 of
Kuchibhotla and Berk (2021).

Now we list some assumptions that will be used in the following theorems.

  1. (A1)

$(X_{i},T_{i},(1-T_{i})R_{i}),i\in\mathcal{I}_{2}$ are independent and
identically distributed random vectors satisfying condition (7).

  2. (A2)

The functions $(\theta,x)\mapsto\widehat{m}(\theta,x)$ and
$x\mapsto\widehat{\pi}(x)$ are bounded, i.e., there exist $m_{0}$ and
$\pi_{0}$ such that for all $\theta\in\mathbb{R}$ and $x\in\mathbb{R}^{d}$,
$|\widehat{m}(\theta,x)|\leq m_{0}$ and $|\widehat{\pi}(x)|\leq\pi_{0}$.

  3. (A3)

The estimator $\widehat{m}(\theta,x)$ is a nondecreasing function of $\theta$.

Assumptions (A1) and (A2) are both standard conditions where we note that
$m^{\star}(\cdot,\cdot)$ is a conditional CDF contained in the unit interval
[0,1]. For assumption (A3), because $m^{\star}(\theta,x)$ is a conditional CDF
which must be monotonically nondecreasing in $\theta$, any given estimator
$\widetilde{m}$ can be improved upon by transforming it into a monotone
estimator $\widehat{m}^{*}$ such that
$\|\widehat{m}^{*}-m\|\leq\|\widetilde{m}-m\|$, see e.g. the first two
properties of Proposition 2 of Chernozhukov et al. (2009). We state their
result in Proposition 1 in S.9 of the supplementary for completeness. Given
this proposition, it is natural that we restrict the estimator
$\widehat{m}(\cdot,\cdot)$ to the class of functions that are non-decreasing
in their first argument. Hence we impose assumption (A3).

Under assumptions (A1)–(A3), we now provide a coverage guarantee for the
prediction set $\widehat{C}_{\alpha}$ returned by Algorithm 1. Observe that
following Lemma 2 and the discussion surrounding (23) and (24), we obtain

$\begin{split}\mathbb{P}\Bigl{(}Y\in\widehat{C}_{\alpha}(X)\mid\mathcal{D}^{\mathrm{tr}},T=1\Bigr{)}-(1-\alpha)~{}&=~{}\frac{\mathbb{P}_{\mathcal{I}_{2}}[\mathrm{IF}(\widehat{r}_{\alpha},X,R,T;\widehat{\pi},\widehat{m})]}{\mathbb{P}(T=1)}\\\
~{}&\quad+\frac{P[\mathrm{IF}(\widehat{r}_{\alpha},X,R,T;\widehat{\pi},\widehat{m})]-\mathbb{P}_{\mathcal{I}_{2}}[\mathrm{IF}(\widehat{r}_{\alpha},X,R,T;\widehat{\pi},\widehat{m})]}{\mathbb{P}(T=1)}\\\
~{}&\quad+\frac{P[\mathrm{IF}(\widehat{r}_{\alpha},X,R,T;\pi^{\star},m^{\star})]-P[\mathrm{IF}(\widehat{r}_{\alpha},X,R,T;\widehat{\pi},\widehat{m})]}{\mathbb{P}(T=1)}\\\
~{}&\geq~{}0+\mathbf{I}+\mathbf{II}.\end{split}$ (26)

Here we use the fact that $\widehat{r}_{\alpha}$ satisfies
$\mathbb{P}_{\mathcal{I}_{2}}[\mathrm{IF}(\widehat{r}_{\alpha},X,R,T;\widehat{\pi},\widehat{m})]\geq
0$. The term $\mathbf{II}$ in (26) can be bounded in absolute value using
Theorem 2. We now provide a bound on $\mathbf{I}$ in Theorem 3 below under
assumptions (A1)–(A3). Theorem 3 actually proves the tail and expectation
bound for

$\sup_{\theta\in\mathbb{R}}|\mathbb{P}_{\mathcal{I}_{2}}[\mathrm{IF}(\theta,X,R,T;\widehat{\pi},\widehat{m})]-P[\mathrm{IF}(\theta,X,R,T;\widehat{\pi},\widehat{m})]|.$

######  Theorem 3.

Under assumption (A1), for any estimators $\widehat{\pi},\widehat{m}$
satisfying assumptions (A2) and (A3), there exists a universal constant
$\mathfrak{C}$ such that for any $\delta>0$,

$\displaystyle\mathbb{P}\left(|\mathbf{I}|\leq\frac{\mathfrak{C}}{\mathbb{P}(T=1)}\sqrt{\frac{(m_{0}+\pi_{0}+1-\alpha)^{2}\log\bigl{(}{1}/{\delta}\bigr{)}+(m_{0}+\pi_{0})^{2}}{|\mathcal{I}_{2}|}}\;\bigg{|}\mathcal{D}_{1}\right)~{}\geq~{}1-\delta.$
(27)

Moreover, there exists a universal constant $\mathfrak{C}^{\prime}$ such that

$\mathbb{E}\left[|\mathbf{I}|\big{|}\mathcal{D}_{1}\right]~{}\leq~{}\frac{\mathfrak{C}^{\prime}}{\mathbb{P}(T=1)}\sqrt{\frac{(m_{0}+\pi_{0}+1-\alpha)^{2}+(m_{0}+\pi_{0})^{2}}{|\mathcal{I}_{2}|}}\leq\frac{\mathfrak{C}^{\prime}}{\mathbb{P}(T=1)}\frac{m_{0}+\pi_{0}+1}{\sqrt{|\mathcal{I}_{2}|}}.$

###### Proof.

See Section S.6 of the appendix. ∎

This theorem is proved using techniques from empirical process theory by
bounding
$\sup_{\theta}\bigl{|}\mathbb{P}_{N}[\mathrm{IF}(\widehat{r}_{\alpha},\dots)]-\mathbb{P}[\mathrm{IF}(\widehat{r}_{\alpha},\dots)]\bigr{|}$.
It gives a convergence rate of
$\mathbb{P}_{N}[\mathrm{IF}(\widehat{r}_{\alpha},\dots)]$ to
$\mathbb{P}[\mathrm{IF}(\widehat{r}_{\alpha},\dots)]$ that scales as
$O(N^{-1/2})$, if $|\mathcal{I}_{1}|\asymp|\mathcal{I}_{2}|\asymp N$.

Using the definition of $\widehat{r}_{\alpha}$ and then combining Theorems 2
and 3 together with (26) yields the following main result.

######  Theorem 4.

Under assumption (A1), for any estimators $\widehat{\pi},\widehat{m}$
satisfying assumptions (A2) and (A3), there exists a universal constant
$\mathfrak{C}$ such that for any $\delta>0$ with probability at least
$1-\delta$,

$\displaystyle\mathbb{P}_{(X,Y)\sim Q_{X}\times P_{Y\mid
X}}\Bigl{(}Y\in\widehat{C}(\widehat{r}_{\alpha};X)\mid\mathcal{D}^{\mathrm{tr}}\Bigr{)}$
$\displaystyle\geq 1-\alpha$ (28)
$\displaystyle\quad-\frac{\|\widehat{\pi}-\pi^{\star}\|_{2}}{\mathbb{P}(T=1)}\sup_{\theta}\|\widehat{m}(\theta,\cdot)-m^{\star}(\theta,\cdot)\|_{2}$
$\displaystyle\quad-\mathfrak{C}\frac{(m_{0}+\pi_{0}+1)}{\mathbb{P}(T=1)}\sqrt{\frac{\log\bigl{(}{1}/{\delta}\bigr{)}+1}{|\mathcal{I}_{2}|}}.$

Moreover,

$\begin{split}\mathbb{P}_{(X,Y)\sim Q_{X}\otimes
P_{Y|X}}\Bigl{(}Y\in\widehat{C}(\widehat{r}_{\alpha};X)\Bigr{)}~{}&\geq~{}(1-\alpha)-\mathbb{E}\left[\frac{\|\widehat{\pi}-\pi^{\star}\|_{2}}{\mathbb{P}(T=1)}\sup_{\theta}\|\widehat{m}(\theta,\cdot)-m^{\star}(\theta,\cdot)\|_{2}\right]\\\
&\qquad-\frac{\mathfrak{C}}{\mathbb{P}(T=1)}\frac{m_{0}+\pi_{0}+1}{\sqrt{|\mathcal{I}_{2}|}}.\end{split}$
(29)

###### Proof.

This is a direct result of Lemma 2, (26), and Theorem 3. ∎

Equation (28) of Theorem 4 provides a coverage guarantee conditional on the
training data and (29) provides a bound on the unconditional coverage. Note
that the slack for the coverage is the sum of two terms: the product bias from
the estimation of $\pi^{\star}$ and $m^{\star}$ and a term of order
$O(N^{-1/2})$ if we assume the two splits $\mathcal{I}_{1}$ and
$\mathcal{I}_{2}$ are of similar size. Proposition 2a of Vovk (2012)
establishes a conditional prediction coverage guarantee which also involves a
slack analogous to the last term of (28). The miscoverage error slacks in (28)
and (29) have clear meaning. The first slack (product of errors) comes from
the double robustness property of $\mathrm{IF}$ and the second slack (of order
$|\mathcal{I}_{2}|^{-1/2}$) comes from approximating $P[\mathrm{IF}(\cdots)]$
with $\mathbb{P}_{\mathcal{I}_{2}}[\mathrm{IF}(\cdots)]$.

######  Remark.

Theorem 4 only provides lower bounds on the (conditional or unconditional on
$\mathcal{D}^{\mathrm{tr}}$) coverage probability. Without a continuity
assumption on the distribution of $R(X,Y)$ when $(X,Y)\sim Q_{X}\otimes
P_{Y|X}$, it is not possible to provide an upper bound. If there is an
$r_{\alpha}$ such that
$P[\mathrm{IF}(r_{\alpha},X,R,T;\pi^{\star},m^{\star})]=0$, then the
conclusions of Theorem 4 can be made two-sided. The condition of existence of
$r_{\alpha}$ which makes $P[\mathrm{IF}(\cdots)]=0$ is same as saying that
there exists an $r_{\alpha}$ such that $\mathbb{P}(R(X,Y)\leq
r_{\alpha}|T=1)=1-\alpha$, i.e., there are no jumps in the distribution of
$R(X,Y)|T=1$ at $(1-\alpha)$-th quantile, or equivalently, the distribution
function takes the value of $1-\alpha$. Under this condition, inequality (29),
for instance, can be strengthened to

$\begin{split}\Bigl{|}\mathbb{P}_{(X,Y)\sim Q_{X}\otimes
P_{Y|X}}\Bigl{(}Y\in\widehat{C}(\widehat{r}_{\alpha};X)\Bigr{)}-(1-\alpha)\Bigr{|}&\leq\mathbb{E}\left[\frac{\|\widehat{\pi}-\pi^{\star}\|_{2}}{\mathbb{P}(T=1)}\sup_{\theta}\|\widehat{m}(\theta,\cdot)-m^{\star}(\theta,\cdot)\|_{2}\right]\\\
&\qquad+\frac{\mathfrak{C}}{\mathbb{P}(T=1)}\frac{m_{0}+\pi_{0}+1}{\sqrt{|\mathcal{I}_{2}|}}.\end{split}$
(30)

Similar strengthening also holds for (28). See the proof of Theorem 4 for
details.

##  6 Methodology without sample splitting

In this section, we provide an alternative methodology that builds upon
Algorithm 1 but is potentially more efficient by avoiding sample splitting.
The procedure is summarized in Algorithm 2.

Input: Training data
$\mathcal{D}^{\text{tr}}=\mathcal{D}^{\text{tr}}_{P}\cup\mathcal{D}^{\text{tr}}_{Q}$;
Coverage probability $1-\alpha$, a training method $\mathcal{A}$ and
estimators $\widehat{\pi},\widehat{m}$, the point for prediction $x$.

Output: A valid prediction set $\widehat{C}_{\alpha}(x)$.

1 Fit the training method $\mathcal{A}$ on $\mathcal{D}^{\text{tr}}$ and using
fitted method $\mathcal{A}$, construct an increasing (nested) sequence of sets
$\\{\mathcal{F}_{t}\\}_{t\in\mathcal{T}}$. Here $\mathcal{T}$ is a subset of
$\mathbb{R}$. The nested sets $\\{\mathcal{F}_{t}\\}_{t\in\mathcal{T}}$ can
depend arbitrarily on $\mathcal{D}^{\text{tr}}$.;

2 For each $i\in[N]$ that satisfies $T_{i}=0$, define the conformal score
$r_{i}~{}=~{}r(Z_{i}):=\inf\\{t\in\mathcal{T}:\,Z_{i}\in\mathcal{F}_{t}\\}.$ ;

3 Fit estimators $\widehat{\pi},\widehat{m}$ on $\mathcal{D}^{\text{tr}}$ and
find the smallest $\widehat{\theta}=\widehat{r}_{\alpha}$ such that
$\mathbb{P}_{N}[\mathrm{IF}({\theta},X,R,T;\widehat{\pi},\widehat{m})]\geq 0$,
where
$\displaystyle\mathbb{P}_{N}\bigl{[}\mathrm{IF}({\theta},X,R,T;\widehat{\pi},\widehat{m})\bigr{]}~{}$
$\displaystyle=~{}\frac{1}{N}\sum_{i=1}^{N}\mathbbm{1}\\{t_{i}=0\\}\widehat{\pi}(x_{i})\big{[}\mathbbm{1}\\{r_{i}\leq{\theta}\\}-\widehat{m}({\theta},x_{i})\bigr{]}$
$\displaystyle\quad+~{}\frac{1}{N}\sum_{i=1}^{N}\mathbbm{1}\\{t_{i}=1\\}[\widehat{m}({\theta},x_{i})-(1-\alpha)].$

return the prediction set
$\widehat{C}_{\alpha}(x):=\\{y:R(x,y)\leq\widehat{r}_{\alpha}\\}$.

Algorithm 2 Full doubly robust prediction

Similar as in (26), let
$\mathbf{I}:=\\{P[\mathrm{IF}(\widehat{r}_{\alpha},X,R,T;\widehat{\pi},\widehat{m})]-\mathbb{P}_{N}[\mathrm{IF}(\widehat{r}_{\alpha},X,R,T;\widehat{\pi},\widehat{m})]\\}/\mathbb{P}(T=1)$.
We now provide a bound on $\mathbf{I}$ in Theorem 5. Note that in addition to
assumptions (A1)–(A3), because we are doing the training and evaluating on the
same dataset, we need some additional assumptions on the classes of estimators
in order to apply results from empirical processes to ensure $\mathbf{I}$
converges to zero.

  1. (A4)

Assume that $\widehat{\pi}(\cdot)$ and its limit $\pi(\cdot)$ are in function
classes $\mathcal{F}_{\pi}$ and $\widehat{m}(\cdot,\cdot)$ and its limit
${m}(\cdot,\cdot)$ are in $\mathcal{F}_{m}$, such that for some
$\alpha_{\pi},\alpha_{m}\geq 0$, the covering numbers satisfy
$\forall\varepsilon>0,$

$\log N(\varepsilon,\mathcal{F}_{\pi},L_{2}(Q))\leq
C\varepsilon^{-\alpha_{\pi}},\text{ and }\log
N(\varepsilon,\mathcal{F}_{m},L_{2}(Q))\leq C\varepsilon^{-\alpha_{m}},$

where $C$ is some constant and $Q$ is any discrete probability measure, and
the covering number $N\bigl{(}\varepsilon,\mathcal{F},L_{r}(Q)\bigr{)}$ is
defined in the same way as in Kearns et al. (1994).

######  Theorem 5.

Under assumption (A1), for any estimators $\widehat{\pi},\widehat{m}$
satisfying assumptions (A2), (A3) and (A4), there exists a universal constant
$\mathfrak{C}$ such that for any $\delta>0$,

$\displaystyle\mathbb{P}\Biggl{(}|\mathbf{I}|\leq\mathfrak{C}\biggl{\\{}$
$\displaystyle N^{-1/(\alpha_{m}\vee
2)}\bigl{(}1+\mathbbm{1}{\\{\alpha_{m}=2\\}}\log
N\bigr{)}+N^{-1/(\alpha_{\pi}\vee
2)}\bigl{(}1+\mathbbm{1}{\\{\alpha_{\pi}=2\\}}\log N\bigr{)}+$ (31)
$\displaystyle\sqrt{\frac{(m_{0}+\pi_{0}+1-\alpha)^{2}\log\bigl{(}\frac{1}{\delta}\bigr{)}}{N}}\biggr{\\}}/\mathbb{P}(T=1)\Biggr{)}\geq
1-\delta,$

where $\vee$ is the maximum operator. Moreover, there exists a universal
constant $\mathfrak{C}^{\prime}$ such that

$\mathbb{E}\left[|\mathbf{I}|\right]~{}\leq~{}\frac{\mathfrak{C}^{\prime}}{\mathbb{P}(T=1)}\biggl{(}N^{-1/(\alpha_{m}\vee
2)}\bigl{(}1+\mathbbm{1}{\\{\alpha_{m}=2\\}}\log
N\bigr{)}+N^{-1/(\alpha_{\pi}\vee
2)}\bigl{(}1+\mathbbm{1}{\\{\alpha_{\pi}=2\\}}\log
N\bigr{)}+\frac{m_{0}+\pi_{0}+1-\alpha}{\sqrt{N}}\biggr{)}.$

The proof of this theorem is in Section S.7 of the appendix.

Using the definition of $\widehat{r}_{\alpha}$ and combining Theorems 2 and 3
together with (26) yields the following result for the final coverage of the
prediction region from Algorithm 2.

######  Theorem 6.

Under assumption (A1), for any estimators $\widehat{\pi},\widehat{m}$
satisfying assumptions (A2), (A3) and (A4), there exists a universal constant
$\mathfrak{C}$ such that for any $\delta>0$ with probability at least
$1-\delta$,

$\displaystyle\mathbb{P}_{(X,Y)\sim Q_{X}\times P_{Y\mid
X}}\Bigl{(}Y\in\widehat{C}(\widehat{r}_{\alpha};X)\mid\mathcal{D}^{\mathrm{tr}}\Bigr{)}\geq
1-\alpha-\frac{1}{\mathbb{P}(T=1)}\biggl{(}\|\widehat{\pi}-\pi^{\star}\|_{2}\sup_{\theta}\|\widehat{m}(\theta,\cdot)-m^{\star}(\theta,\cdot)\|_{2}$
(32) $\displaystyle\quad+N^{-1/(\alpha_{m}\vee
2)}\bigl{(}1+\mathbbm{1}{\\{\alpha_{m}=2\\}}\log
N\bigr{)}+N^{-1/(\alpha_{\pi}\vee
2)}\bigl{(}1+\mathbbm{1}{\\{\alpha_{\pi}=2\\}}\log
N\bigr{)}+\mathfrak{C}(m_{0}+\pi_{0}+1)\sqrt{\frac{\log\bigl{(}{1}/{\delta}\bigr{)}+1}{N}}\biggr{)}.$

Moreover,

$\begin{split}&\mathbb{P}_{(X,Y)\sim Q_{X}\otimes
P_{Y|X}}\Bigl{(}Y\in\widehat{C}(\widehat{r}_{\alpha};X)\Bigr{)}-(1-\alpha)\\\
&\qquad\geq-\frac{\mathfrak{C}}{\mathbb{P}(T=1)}\biggl{(}N^{-1/(\alpha_{m}\vee
2)}\bigl{(}1+\mathbbm{1}{\\{\alpha_{m}=2\\}}\log
N\bigr{)}+N^{-1/(\alpha_{\pi}\vee
2)}\bigl{(}1+\mathbbm{1}{\\{\alpha_{\pi}=2\\}}\log
N\bigr{)}+\frac{m_{0}+\pi_{0}+1}{\sqrt{N}}\biggr{)}\\\
&\qquad\quad-\mathbb{E}\left[\frac{\|\widehat{\pi}-\pi^{\star}\|_{2}}{\mathbb{P}(T=1)}\sup_{\theta}\|\widehat{m}(\theta,\cdot)-m^{\star}(\theta,\cdot)\|_{2}\right].\end{split}$
(33)

###### Proof.

This is a direct result of Lemma 2, (26), and Theorem 5. ∎

Comparing the results of Theorem 6 with Theorem 4 in terms of the slack in the
miscoverage probability, we notice that the full data based prediction set can
have larger miscoverage error than the split data version. However, in terms
of how close $\widehat{r}_{\alpha}$ is to $r_{\alpha}$, we expect the full
data version to perform better compared to the split data version. This is
expected because with full data version the quantile estimator is based on $N$
observations rather than $|\mathcal{I}_{2}|$ observations which is smaller
than $N$. The variance of the quantile estimator can be better up to a
constant of $N/|\mathcal{I}_{2}|$. We do not pursue these variance comparisons
for $\widehat{r}_{\alpha}$ as that is not our goal.

##  7 Simulation Studies

In practice, unless $R(\cdot,\cdot)$ is a map that is independent of data, we
further split $\mathcal{D}_{1}$ of Algorithm 1 and use the first split to
train $R(\cdot,\cdot)$, while the second split is used to estimate the two
nuisance parameters $\pi^{\star}(\cdot)$ and $m^{\star}(\cdot,\cdot)$. We
include the two algorithms in Algorithm 1 and Algorithm 2 under the name “DRP
w. three splits” and “DRP w. full data” respectively. We also include
simulation results where $R$, $\widehat{m}$ and $\widehat{\pi}$ are trained on
the same split and the remainder of the data is used for validation under the
name “DRP w. two splits” in Section S.10.1 of the appendix. For both synthetic
data and the real data set we use two kinds of score functions to estimate
$R(\cdot,\cdot)$,

  1. 1.

absolute residual score $|y-\widehat{\mu}(x)|$ where
$\mu(x):=\mathbb{E}(Y|X=x)$ is estimated by ridge regression,

  2. 2.

conditional quantile regression 222the cutting edge algorithm employed in Lei
and Candès (2020) paper;

The nuisance parameters $\pi^{\star}$ and $m^{\star}(\cdot,\cdot)$ are
estimated through SuperLearner333SuperLearner uses cross-validation to
estimate the performance of multiple machine learning models and then creates
an optimal weighted average of those models using the test data. This approach
has been proven to be asymptotically as accurate as the best possible
prediction algorithm that is tested. For details please refer to Polley and
van der Laan (2010). that includes both RandomForest and generalized linear
model (GLM). To avoid numerical issues, we clip the propensity score at 0.99
to prevent $\widehat{\pi}$ from becoming unbounded. The proposed methods are
compared against the weighted conformal prediction (WCP) method proposed in
Tibshirani et al. (2019). Note that in the weighted conformal prediction
method, the prediction interval is given by

$\widehat{C}_{n}(x)=\mu_{0}(x)\pm\text{ Quantile
}\biggl{(}1-\alpha;\sum_{i=1}^{n}p_{i}^{w}(x)\delta_{\left|Y_{i}-\mu_{0}\left(X_{i}\right)\right|}+p_{n+1}^{w}(x)\delta_{\infty}\biggr{)},$
(34)

where $p^{w}(x)$ is a function which depends on the likelihood ratio between
the two covariate distributions, or $\pi^{\star}(x)$. Therefore, when the
distribution shift is too “large”, i.e. $p^{w}_{n+1}(x)$ is larger than
$\alpha$, the width becomes $\infty$. And indeed we observe infinite widths
for this method over $50\%$ of the time across all the predicted points and
the Monte Carlo replications on synthetic data and over $90\%$ of the time on
real data. For illustrative purposes, we truncate the width for WCP at 10 and
50 on synthetic data and real data respectively, and demonstrate the mean
width from 500 Monte Carlo simulations.

It is shown in Lei and Candès (2020) that the CQR score would guarantee
asymptotic conditional coverage
$\mathbb{P}\bigl{(}Y_{f}\in\widehat{C}_{N,\alpha}(X_{f})|X_{f}=x_{f}\bigr{)}$,
we also conduct an experiment using this score for our method under the
setting of Section 7.2, where we specify 200 test points of $X_{f}$ (generated
from standard normal distributions $N(0,I_{4})$) and for each test point, 100
$Y_{f}$’s following the distribution (35) are generated to test if they fall
into the prediction sets trained by our method and WCP, where we report the
average (reflected points) for each test point in Figure 3,which shows the
coverage and width for each test point that is represented by the $L_{2}$ norm
(the X axis). We also fit a smoothing spline (with default parameters of the R
function smooth.spline) with these points.

###  7.1 Real data

We demonstrate the use of conformal prediction in the covariate shift setting
in an empirical example. We re-analyze the data set used in Tibshirani et al.
(2019) which is the airfoil data set from the UCI Machine Learning Repository
which contains $N=1503$ observations of a response $Y$ (scaled sound pressure
level of NASA airfoils), and a vector of covariates $X$ with $d=5$ dimension
(log frequency, angle of attack, chord length, free-stream velocity, and
suction side log displacement thickness). Label missingness was then generated
as in Tibshirani et al. (2019) using a propensity score model analogous to the
one specified in Section 7.2 below.

Figure 1: Histograms of coverage and width of Double Robust Prediction (DRP)
and Weighted Conformal Prediction (WCP) on real data through either absolute
residual score or the CQR score, where the width is truncated at 50 when WCP
produces infinity width. Figure 2: Histograms of coverage and width of Doubly
Robust Prediction (DRP) and Weighted Conformal Prediction (WCP) on synthetic
data using the absolute residual score. The width is truncated at 10 for WCP.

###  7.2 Synthetic data

Here we use the setting from Kang and Schafer (2007) where for each unit
$i=1,\ldots,n=2000$, suppose that
$\left(x_{i1},x_{i2},x_{i3},x_{i4}\right)^{\top}$ is independently distributed
as $N(0,I_{4})$ where $I_{4}$ is the $4\times 4$ identity matrix. The
$y_{i}$’s are generated as

$y_{i}=210+27.4x_{i1}+13.7x_{i2}+13.7x_{i3}+13.7x_{i4}+\varepsilon_{i},$ (35)

where $\varepsilon_{i}\sim N(0,1)$, and the true propensity scores are

$\mathbb{P}(T=1|x_{i})=\operatorname{expit}\left(-x_{i1}+0.5x_{i2}-0.25x_{i3}-0.1x_{i4}\right),\text{
where
}\mathrm{expit}(x_{i}^{\top}\alpha)=\frac{\exp(x_{i}^{\top}\alpha)}{1+\exp(x_{i}^{\top}\alpha)}.$

###  7.3 Simulation results

The mean coverage and width from 500 Monte Carlo simulations using the
absolute residual score are shown in Table 1, where the middle column
corresponds to synthetic data results and the rightmost column to real data
results. Figure 1 display histograms of coverage and width from 500 runs on
real data using either the absolute residual score and the CQR score while 2
shows the use of absolute residual score. For synthetic data, we keep the
results using the two scores separate, with the absolute residual score in
Figure 2 and the CQR score to Appendix S.10.2 for easier comparison between
our DRP methods and the WCP method. Below we make a few observations on
simulation results:

  * •

WCP produces wider width and therefore, tends to over-cover by a considerable
amount (by more than $7\%$ over the nominal coverage of $90\%$).

  * •

Doubly robust prediction with full data and multiple splits have similar
performance with valid coverage. In practice, because using three splits
guaranteed nominal coverage in sufficiently large sample with fewest
assumptions, we recommend this approach as it does not appear to suffer much
efficiency loss.

  * •

For the conditional coverage simulation results described in the last
paragraph at the start of Section 7, we see that DRP has similar coverage and
much smaller width compared to WCP. Note that they both attain desired
coverage when the norm of the test data is less than 2, which is what the
$L_{2}$ norm of a standard 4-dimensional normal r.v. concentrates to. As the
norm gets past 2, there are less data points and hence we get fewer
observations.

Mean coverage and width | Synthetic data | Real data
---|---|---
from 500 monte carlo runs | Coverage | Width | Coverage | Width
DRP w. full data | 0.90 | 3.29 | 0.94 | 27.85
DRP w. three splits | 0.90 | 3.30 | 0.90 | 25.79
DRP w. two splits | 0.90 | 3.29 | 0.88 | 25.19
WCP | 0.97 | 7.41 | 0.99 | 47.71
Table 1: Coverage and width of DRP and WCP on synthetic and real data.
Clearly, DRP improves on WCP in terms of width while maintaining coverage
close to the nominal level of 0.9. Figure 3: Estimated conditional coverage
and length at 200 test points with Split DRP and WCP using the CQR score. The
setting is described in the last paragraph at the start of Section 7. The
points are the average of the coverage and width from 100 Monte Carlo
simulations and the lines are drawn by fitting a smoothing splines for those
points.

##  8 Aggregation of prediction sets

When, as typically the case in practice, multiple training methods are
available, we propose to combine the proposed approach with that of Yang and
Kuchibhotla (2021) in order to construct a better prediction set with smaller
width. The detailed procedure is given in Algorithm 3 and can be shown by
arguments given in Yang and Kuchibhotla (2021) to retain coverage validity
while attaining the smallest width in large samples with high probability.

Input: Training data
$\mathcal{D}^{\text{tr}}=\mathcal{D}^{\text{tr}}_{P}\cup\mathcal{D}^{\text{tr}}_{Q}$
split into $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$, where
$\mathcal{D}_{1}=\\{Z_{i}\in\mathcal{D}^{\text{tr}},i\in\mathcal{I}_{1}\\}$
and
$\mathcal{D}_{2}=\\{Z_{i}\in\mathcal{D}^{\text{tr}},i\in\mathcal{I}_{2}\\}$;
Estimators $\widehat{\pi},\widehat{m}$, and training methods
$\mathcal{A}_{k},k\in[K]$; The prediction point $x$.

Output: A valid prediction set $\widehat{C}^{\mathrm{EFCP}}_{\alpha}(x)$ with
smallest width.

1 Fit training methods $\mathcal{A}_{1},\ldots,\mathcal{A}_{K}$ on
$\mathcal{D}_{1}$ and for each fitted method $\mathcal{A}_{k}$, construct an
increasing (nested) sequence of sets
$\\{\mathcal{F}_{t}^{(k)}\\}_{t\in\mathcal{T}}$. Here $\mathcal{T}$ is a
subset of $\mathbb{R}$.;

2 For each $i\in\mathcal{I}_{2}$ that satisfies $T_{i}=0$, define the
conformal score
$r_{k}(Z_{i}):=\inf\\{t\in\mathcal{T}:\,Z_{i}\in\mathcal{F}_{t}^{(k)}\\}.$ ;

3 Fit estimators $\widehat{\pi},\widehat{m}$ on $\mathcal{D}_{1}$ and solve
for $\widehat{\theta}=\widehat{r}_{\alpha,k}$ as the solution to
$\mathbb{P}_{\mathcal{I}_{2}}\big{[}\mathrm{IF}(\widehat{\theta},\mathcal{D}_{2},r_{k};\widehat{\pi},\widehat{m})\bigr{]}=0$,
where
$\displaystyle\mathbb{P}_{\mathcal{I}_{2}}[\mathrm{IF}(\widehat{\theta},X,R,T;\widehat{\pi},\widehat{m})]~{}$
$\displaystyle=~{}\frac{1}{|\mathcal{I}_{2}|}\sum_{i\in\mathcal{I}_{2}}\mathbbm{1}\\{t_{i}=0\\}\widehat{\pi}(x_{i})\big{[}\mathbbm{1}\\{r_{i}\leq\widehat{\theta}\\}-\widehat{m}(\widehat{\theta},x_{i})\bigr{]}$
$\displaystyle\quad+~{}\frac{1}{|\mathcal{I}_{2}|}\sum_{i\in\mathcal{I}_{2}}\mathbbm{1}\\{t_{i}=1\\}[\widehat{m}(\widehat{\theta},x_{i})-(1-\alpha)].$
Compute the corresponding conformal prediction set as
$\widehat{C}_{k}(x)~{}:=~{}\\{y:\,r_{k}(x,y)\leq\widehat{r}_{\alpha,k}\\}.$;

4 Set  $\widehat{k}:=\operatorname*{arg\,min}_{1\leq k\leq
K}\,\mbox{Width}(\widehat{C}_{k}(x)).$ Here $\mbox{Width}(\cdot)$ can be any
measure of width or volume of a prediction set. The quantity $\widehat{k}$
need not be unique and any minimizer can be chosen.;

return the prediction set $\widehat{C}_{\widehat{k}}$ as
$\widehat{C}^{\mathrm{EFCP}}_{\alpha}$.

Algorithm 3 Efficient doubly robust prediction

Table 2 reports results from a simulation study comparing the proposed doubly
robust prediction algorithm to select an optimal tuning parameter for ridge
regression-based conformal score, with cross-validation aimed at minimizing
the ridge regression MSE. The simulation results confirm the proposed
algorithm’s ability to preserve marginal prediction coverage while optimizing
prediction interval width, both in synthetic and real data sets.

Mean coverage and width | Synthetic data | Real data
---|---|---
from 500 monte carlo runs | Coverage | Width | Coverage | Width
Efficient DRP | 0.90 | 3.30 | 0.89 | 27.30
DRP w. CV | 0.89 | 3.32 | 0.84 | 18.61
Table 2: Coverage and width of efficient doubly robust prediction and doubly
robust prediction with cross-validation on synthetic and real data. Efficient
DRP improves on CV in terms of width while maintaining coverage close to the
nominal level of 0.9.

##  9 Sensitivity analysis for latent covariate shift

Thus far, we have assumed that the covariate shift problem is primarily due to
observed covariates, which we have denoted explainable covariate shift (or
equivalently that outcomes for the target population are missing at random),
an assumption that cannot be confirmed empirically without invoking a
different non-testable assumption. In this section, we relax this assumption
$P_{Y|X}=Q_{Y|X}$ and propose a sensitivity analysis for obtaining doubly
robust calibrated prediction sets accounting for a latent covariate shift
problem encoded in a sensitivity parameter. We define a sensitivity function
as

$\displaystyle\gamma^{\star}(x,y)~{}=~{}\log\left(\frac{P_{Y=y|X}}{Q_{Y=y|X}}\bigg{/}\frac{P_{Y=y_{0}|X}}{Q_{Y=y_{0}|X}}\right),$

where $y_{0}$ is any baseline value for $Y$. Here $\gamma^{\star}(x,y)$ is the
sensitivity analysis function encoding a hypothetical departure from the
assumption that $X$ suffices to account for the covariate shift problem with
$\gamma^{\star}(x,y_{0})=0$, and $\gamma^{\star}(x,y)=0$ for all $y$ recovers
the standard assumption of explainable covariate shift. For simplicity we take
$y_{0}=0$. Our sensitivity analysis is inspired by a semiparametric approach
for accounting for data missing not at random (MNAR), due to Chapter 5 of
Robins et al. (2000) in the missing data literature. Formally, the sensitivity
function $\gamma^{\star}$ can be represented in the missing data notation as

$\displaystyle\gamma^{\star}(x,y)~{}=~{}\log\frac{\mathbb{P}(T=0|X=x,Y=y)\mathbb{P}(T=1|X,Y=0)}{\mathbb{P}(T=0|X,Y=0)\mathbb{P}(T=1|X=x,Y=y)}.$

For any three functions $\eta(\cdot)$, $m(\cdot,\cdot)$, and
$\gamma(\cdot,\cdot)$, let

$\displaystyle\mathrm{IF}(\theta,x,y,r,t;\eta,m,\gamma):$
$\displaystyle=\mathbbm{1}\\{t=0\\}\exp\bigl{\\{}-\eta(x)-\gamma(x,y)\bigr{\\}}\Big{[}\mathbbm{1}\\{r\leq\theta\\}-m(\theta,x)\Big{]}$
(36)
$\displaystyle\quad+\mathbbm{1}\\{t=1\\}\Big{[}m(\theta,x)-(1-\alpha)\Big{]},$

where the two nuisance functions are

$\begin{split}\eta^{\star}(x)~{}&:=~{}\log\frac{\mathbb{P}(T=0|X=x,Y=0)}{\mathbb{P}(T=1|X=x,Y=0)}\\\
m^{\star}(\theta,x)~{}&:=~{}\mathbb{P}(R\leq\theta|X=x,T=1).\end{split}$ (37)

Per the sensitivity analysis framework, we assume
$\gamma^{\star}(x,y)=\gamma(x,y)$ is known. See also Example 2 of Tsiatis
(2014) for some explanations on why this is a nonparametric identified model.
In this framework we have the following theorem.

######  Theorem 7.

Under the assumption that
$\mathbb{E}\frac{\mathbb{P}^{2}(T=1|X,Y)}{\mathbb{P}(T=0|X,Y)}$ is finite and
that the density of the conditional distribution of $R|T=1$ at $r_{\alpha}$ is
bounded away from zero, the efficient influence function of $r_{\alpha}$, the
$(1-\alpha)-$quantile of $R|T=1$ in the nonparametric model which allows the
distribution of $Z$ to remain unrestricted, and the conditional log odds ratio
function relating $T$ to $R$ given $X$ is known to equal $\gamma^{\star}$ is
given up to a proportionality constant by

$\displaystyle\psi(z)=\mathrm{IF}(r_{\alpha},x,y,r,t;\eta^{\star},m^{\star},\gamma^{\star})$
$\displaystyle=\mathbbm{1}\\{t=0\\}\exp\bigl{\\{}-\eta^{\star}(x)-\gamma^{\star}(x,y)\bigr{\\}}\Big{[}\mathbbm{1}\\{r\leq
r_{\alpha}\\}-m^{\star}(r_{\alpha},x)\Big{]}$ (38)
$\displaystyle\quad+\mathbbm{1}\\{t=1\\}\Big{[}m^{\star}(r_{\alpha},x)-(1-\alpha)\Big{]}$

Furthermore, the moment function
$\mathrm{IF}(r_{\alpha},x,y,r,t;\eta,m,\gamma^{\star})$ satisfies the double
robustness property that

$\mathbb{E}\bigl{[}\mathrm{IF}(r_{\alpha},x,y,r,t;\eta,m,\gamma^{\star})\bigr{]}=0,$
(39)

if either $\eta=\eta^{\star}$ or $m=m^{\star}$.

###### Proof.

See Section S.11 for a proof. ∎

Theorem 7 gives the efficient influence function of $r_{\alpha}$, which
provides a moment equation for $r_{\alpha}$. And this can be leveraged as an
estimating equation similar to Lemma 1 and Theorem 4 to yield a valid
prediction set with the product bias from estimating the nuisance functions.
In this case, estimation of nuisance functions is not straightforward and not
further pursued in this paper. We also include a review of existing literature
on this topic. The recent paper Jin et al. (2021) generalized the covariate
shift setting to distributional shift where the joint distribution of
covariates and response can be different between the test and training data.
The motivation is sensitivity analysis for individual treatment effects. They
proposed a robust conformal prediction algorithm that builds upon the weighted
conformal inference method from Tibshirani et al. (2019). The paper derives
prediction sets that achieve marginal coverage if the propensity score is
known exactly, allowing for latent covariate shift of a magnitude bounded by a
known sensitivity parameter. They also proposed a second algorithm where the
coverage is attained with high probability $1-\delta$ conditional on the
training data, and showed their results provide tight prediction sets in some
cases that cannot be improved upon under their assumptions. The coverage bias
of both methods is first order. Both their methods require the test point to
be specified in advance so that they can compute a different quantile for
every new $x$ and thus may be computationally intensive. In terms of
sensitivity analysis for unmeasured confounding, we note that their paper
considers a different sensitivity framework than ours, as they posit the
existence of an unobserved confounder $U$ which together with $X$ completely
accounts for confounding. They encode the magnitude of unmeasured confounding
in terms of upper and lower bounds for the likelihood ratio of density of $U$
in the treated and untreated samples over the support of $U$. Thus, their
sensitivity parameter appears to capture both the extent of residual
confounding, but also reflects aspects of the density of $U$ that may not be
of scientific interest, and for which an investigator may not have any prior
information. An implication of this choice of parameterization is that for a
given sensitivity bound, the bound may reflect small amount of confounding
over a wide support of $U$, or a large amount of confounding, over a narrow
support of $U$, rendering such sensitivity analysis difficult to interpret. In
contrast, we favor the approach of Robins et al. (2000) with a more direct
sensitivity analysis, which in the counterfactual setting, encodes the
departure from unconfoundedness in terms of a likelihood ratio for the
counterfactual outcome in view of the treated and control arms conditional on
observed covariates. Yin et al. (2021) also studied the sensitivity analysis
of ITE using a similar approach as the first method proposed in Jin et al.
(2021) while their method of analysis offers a different perspective.

##  10 Discussion

This paper has proposed three separate algorithms to construct prediction
regions which are adaptive to unknown covariate shift between a population
from which labeled data are available and an unlabeled population for which
outcome prediction is in view. Our three proposed methods have been described
as “Split doubly robust prediction”, “Full doubly robust prediction” and
“Efficient doubly robust prediction”. The paper provided a rigorous analysis
of the coverage properties of these algorithms, notably establishing that all
have coverage bias of a product form and providing formal conditions under
which all are asymptotically well-calibrated, in the sense of attaining the
nominal coverage rate in large samples. “Split doubly robust prediction” has
coverage guarantees in large samples under minimal conditions, but requires
one to use a non-negligible subset of the data for training; in contrast,
“Full doubly robust prediction” uses the entire data set both for training and
prediction, and attains the nominal coverage in large samples under relatively
stronger conditions. “Efficient doubly robust prediction” combines “Split
doubly robust prediction” and the EFCP algorithm from Yang and Kuchibhotla
(2021), which is empirically shown to potentially outperform a standard cross-
validation approach. We conjecture that the proposed efficient doubly robust
prediction algorithm is nearly as efficient as an oracle with a priori
knowledge of the optimal prediction interval, although formally proving this
result is left to future work. An important advantage of our framework is that
its large sample efficiency and validity guarantees hold for any collection of
machine learning techniques and their respective tuning parameters, under the
relatively mild requirement that at least one of two estimated nuisance
functions is consistent, without necessarily requiring fast convergence rates
for the latter.

Another important contribution of the paper is to draw upon a key equivalence
between the explainable covariate shift problem, the MAR assumption in the
missing data literature, and the notion of unconfoundedness in the causal
inference literature, to develop a sensitivity analysis approach to evaluate
the extent to which prediction regions might be impacted by hypothetical
departures from this assumption. Notably, the proposed methods readily extend
to accommodate such sensitivity analysis via a slight modification of our
procedure to incorporate a sensitivity parameter, without compromising the
product bias or double robustness property of the approach. Fully developing
prediction inference for this sensitivity analysis framework however requires
care in estimation of nuisance functions which, due to space limitation, we
plan to consider in future work. Overall, this paper reveals and leverages
deep connections between modern literatures of semiparametric theory, missing
data and causal inference, and emerging methods for well-calibrated prediction
inference. To the best of our knowledge, such connections have previously not
been drawn upon as deliberately as shown to be possible in this work which we
hope will generate both interest and further developments towards even more
robust and efficient well-calibrated prediction. One possible line of future
investigation might be to build on recent theory of higher order influence
functions due to Robins and colleagues, see e.g. Robins et al. (2008) and
Robins et al. (2017), which in principle could be used to reduce the second
order product bias obtained in this paper to a higher order product bias,
therefore potentially improving finite sample coverage over a wider range of
regimes.

## References

  * Bang and Robins (2005) Heejung Bang and James M Robins.  Doubly robust estimation in missing data and causal inference models.  _Biometrics_ , 61(4):962–973, 2005.
  * Barber et al. (2019) Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani.  The limits of distribution-free conditional predictive inference.  _arXiv:1903.04684_ , 2019.
  * Barber et al. (2021) Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani.  Predictive inference with the jackknife+.  _The Annals of Statistics_ , 49(1):486–507, 2021\.
  * Barber et al. (2022) Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani.  Conformal prediction beyond exchangeability.  _arXiv preprint arXiv:2202.13415_ , 2022.
  * Bickel et al. (1993) Peter J Bickel, Chris AJ Klaassen, Peter J Bickel, Ya’acov Ritov, J Klaassen, Jon A Wellner, and YA’Acov Ritov.  _Efficient and adaptive estimation for semiparametric models_ , volume 4.  Springer, 1993.
  * Bickel et al. (2009) Steffen Bickel, Michael Brückner, and Tobias Scheffer.  Discriminative learning under covariate shift.  _Journal of Machine Learning Research_ , 10(9), 2009.
  * Candès et al. (2021) Emmanuel J Candès, Lihua Lei, and Zhimei Ren.  Conformalized survival analysis.  _arXiv preprint arXiv:2103.09763_ , 2021.
  * Cattaneo et al. (2021) Matias D Cattaneo, Yingjie Feng, and Rocio Titiunik.  Prediction intervals for synthetic control methods.  _Journal of the American Statistical Association_ , 116(536):1865–1880, 2021.
  * Cauchois et al. (2020) Maxime Cauchois, Suyash Gupta, Alnur Ali, and John C Duchi.  Robust validation: Confident predictions even when distributions shift.  _arXiv preprint arXiv:2008.04267_ , 2020.
  * Chakrabortty and Cai (2018) Abhishek Chakrabortty and Tianxi Cai.  Efficient and adaptive linear regression in semi-supervised settings.  _The Annals of Statistics_ , 46(4):1541–1572, 2018.
  * Chen et al. (2016) Xiangli Chen, Mathew Monfort, Anqi Liu, and Brian D Ziebart.  Robust covariate shift regression.  In _Artificial Intelligence and Statistics_ , pages 1270–1279. PMLR, 2016.
  * Chernozhukov et al. (2009) Victor Chernozhukov, Ivan Fernandez-Val, and Alfred Galichon.  Improving point and interval estimators of monotone functions by rearrangement.  _Biometrika_ , 96(3):559–575, 2009.
  * Chernozhukov et al. (2018) Victor Chernozhukov, Kaspar Wüthrich, and Zhu Yinchu.  Exact and robust conformal inference methods for predictive machine learning with dependent data.  In _Conference On Learning Theory_ , pages 732–749. PMLR, 2018.
  * Chernozhukov et al. (2021a) Victor Chernozhukov, Kaspar Wüthrich, and Yinchu Zhu.  Distributional conformal prediction.  _Proceedings of the National Academy of Sciences_ , 118(48), 2021a.
  * Chernozhukov et al. (2021b) Victor Chernozhukov, Kaspar Wüthrich, and Yinchu Zhu.  An exact and robust conformal inference method for counterfactual and synthetic controls.  _Journal of the American Statistical Association_ , 116(536):1849–1864, 2021b.
  * Fan and Park (2007) Y Fan and S Park.  Sharp bounds on the distribution of the treatment effect and their statistical inference,% forthcoming in econometric theory.  2007\.
  * Gibbs and Candès (2021) Isaac Gibbs and Emmanuel Candès.  Adaptive conformal inference under distribution shift.  _arXiv preprint arXiv:2106.00170_ , 2021.
  * Guttman (1967) Irwin Guttman.  Statistical tolerance regions.  _Classical and Bayesian_ , 1967.
  * Györfi and Walk (2020) László Györfi and Harro Walk.  Nearest neighbor based conformal prediction.  _Pub. Inst. Stat. Univ. Paris_ , Special issue in honour of Denis Bosq’s 80th birthday(63):173–190, 2020.
  * Hahn (1998) Jinyong Hahn.  On the role of the propensity score in efficient semiparametric estimation of average treatment effects.  _Econometrica_ , pages 315–331, 1998.
  * Hampel (1974) Frank R Hampel.  The influence curve and its role in robust estimation.  _Journal of the american statistical association_ , 69(346):383–393, 1974.
  * Jin et al. (2021) Ying Jin, Zhimei Ren, and Emmanuel J Candès.  Sensitivity analysis of individual treatment effects: A robust conformal inference approach.  _arXiv preprint arXiv:2111.12161_ , 2021.
  * Kang and Schafer (2007) Joseph DY Kang and Joseph L Schafer.  Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data.  _Statistical science_ , pages 523–539, 2007.
  * Kearns et al. (1994) Michael J Kearns, Umesh Virkumar Vazirani, and Umesh Vazirani.  _An introduction to computational learning theory_.  MIT press, 1994.
  * Kim et al. (2020) Byol Kim, Chen Xu, and Rina Foygel Barber.  Predictive inference is free with the jackknife+-after-bootstrap.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Kivaranovic et al. (2020) Danijel Kivaranovic, Robin Ristl, Martin Posch, and Hannes Leeb.  Conformal prediction intervals for the individual treatment effect.  _arXiv preprint arXiv:2006.01474_ , 2020.
  * Kpotufe and Martinet (2018) Samory Kpotufe and Guillaume Martinet.  Marginal singularity, and the benefits of labels in covariate-shift.  In _Conference On Learning Theory_ , pages 1882–1886. PMLR, 2018\.
  * Krishnamoorthy and Mathew (2009) Kalimuthu Krishnamoorthy and Thomas Mathew.  _Statistical tolerance regions: theory, applications, and computation_ , volume 744.  John Wiley & Sons, 2009.
  * Kuchibhotla and Berk (2021) Arun K Kuchibhotla and Richard A Berk.  Nested conformal prediction sets for classification with applications to probation data.  _arXiv preprint arXiv:2104.09358_ , 2021.
  * Kuchibhotla (2020) Arun Kumar Kuchibhotla.  Exchangeability, conformal prediction, and rank tests.  _arXiv preprint arXiv:2005.06095_ , 2020.
  * Lei and Wasserman (2014) Jing Lei and Larry Wasserman.  Distribution-free prediction bands for non-parametric regression.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76(1):71–96, 2014.
  * Lei et al. (2013) Jing Lei, James Robins, and Larry Wasserman.  Distribution-free prediction sets.  _Journal of the American Statistical Association_ , 108(501):278–287, 2013.
  * Lei et al. (2018) Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J Tibshirani, and Larry Wasserman.  Distribution-free predictive inference for regression.  _Journal of the American Statistical Association_ , 113(523):1094–1111, 2018.
  * Lei and Candès (2020) Lihua Lei and Emmanuel J Candès.  Conformal inference of counterfactuals and individual treatment effects.  _arXiv preprint arXiv:2006.06138_ , 2020.
  * Massart (1990) Pascal Massart.  The tight constant in the dvoretzky-kiefer-wolfowitz inequality.  _The annals of Probability_ , pages 1269–1283, 1990.
  * Newey (1990) Whitney K Newey.  Semiparametric efficiency bounds.  _Journal of applied econometrics_ , 5(2):99–135, 1990.
  * Oliveira et al. (2022) Roberto I Oliveira, Paulo Orenstein, Thiago Ramos, and João Vitor Romano.  Split conformal prediction for dependent data.  _arXiv preprint arXiv:2203.15885_ , 2022.
  * Papadopoulos et al. (2002) Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman.  Inductive confidence machines for regression.  In _European Conference on Machine Learning_ , pages 345–356. Springer, 2002.
  * Park et al. (2021) Sangdon Park, Edgar Dobriban, Insup Lee, and Osbert Bastani.  Pac prediction sets under covariate shift.  _arXiv preprint arXiv:2106.09848_ , 2021.
  * Pfanzagl (1976) J Pfanzagl.  Investigating the quantile of an unknown distribution.  In _Contribution to Applied Statistics_ , pages 111–126. Springer, 1976.
  * Politis (2015) Dimitris N Politis.  _Model-Free Prediction and Regression: A Transformation-Based Approach to Inference_.  Springer, 2015.
  * Polley and van der Laan (2010) Eric C. Polley and Mark J. van der Laan.  Super learner in prediction.  2010\.
  * Qiu et al. (2022) Hongxiang Qiu, Edgar Dobriban, and Eric Tchetgen Tchetgen.  Distribution-free prediction sets adaptive to unknown covariate shift.  _arXiv preprint arXiv:2203.06126_ , 2022.
  * Quiñonero-Candela et al. (2008) Joaquin Quiñonero-Candela, Masashi Sugiyama, Anton Schwaighofer, and Neil D Lawrence.  _Dataset shift in machine learning_.  Mit Press, 2008.
  * Reddi et al. (2015) Sashank Reddi, Barnabas Poczos, and Alex Smola.  Doubly robust covariate shift correction.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 29, 2015.
  * Reeve et al. (2021) Henry WJ Reeve, Timothy I Cannings, and Richard J Samworth.  Adaptive transfer learning.  _The Annals of Statistics_ , 49(6):3618–3649, 2021.
  * Robins et al. (2008) James Robins, Lingling Li, Eric Tchetgen, and Aad van der Vaart.  Higher order influence functions and minimax estimation of nonlinear functionals.  In _Probability and statistics: essays in honor of David A. Freedman_ , pages 335–421. Institute of Mathematical Statistics, 2008.
  * Robins (2000) James M Robins.  Robust estimation in sequentially ignorable missing data and causal inference models.  In _Proceedings of the American Statistical Association_ , volume 1999, pages 6–10. Indianapolis, IN, 2000.
  * Robins et al. (2000) James M Robins, Andrea Rotnitzky, and Daniel O Scharfstein.  Sensitivity analysis for selection bias and unmeasured confounding in missing data and causal inference models.  In _Statistical models in epidemiology, the environment, and clinical trials_ , pages 1–94. Springer, 2000.
  * Robins et al. (2017) James M Robins, Lingling Li, Rajarshi Mukherjee, Eric Tchetgen Tchetgen, and Aad van der Vaart.  Minimax estimation of a functional on a structured high-dimensional model.  _The Annals of Statistics_ , 45(5):1951–1987, 2017.
  * Romano et al. (2019) Yaniv Romano, Evan Patterson, and Emmanuel Candes.  Conformalized quantile regression.  In _Advances in Neural Information Processing Systems_ , pages 3543–3553, 2019.
  * Romano et al. (2020) Yaniv Romano, Matteo Sesia, and Emmanuel J Candès.  Classification with valid and adaptive coverage.  In _Advances in Neural Information Processing Systems_ , pages 3538–3548, 2020.
  * Rubin (1974) Donald B Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of educational Psychology_ , 66(5):688, 1974.
  * Rubin (1990) Donald B Rubin.  Formal mode of statistical inference for causal effects.  _Journal of statistical planning and inference_ , 25(3):279–292, 1990.
  * Scharfstein et al. (1999) Daniel O Scharfstein, Andrea Rotnitzky, and James M Robins.  Adjusting for nonignorable drop-out using semiparametric nonresponse models.  _Journal of the American Statistical Association_ , 94(448):1096–1120, 1999.
  * Shah and Peters (2020) Rajen D Shah and Jonas Peters.  The hardness of conditional independence testing and the generalised covariance measure.  _The Annals of Statistics_ , 48(3):1514–1538, 2020.
  * Solari and Djordjilović (2022) Aldo Solari and Vera Djordjilović.  Multi split conformal prediction.  _Statistics & Probability Letters_, 184:109395, 2022.
  * Splawa-Neyman et al. (1990) Jerzy Splawa-Neyman, Dorota M Dabrowska, and TP Speed.  On the application of probability theory to agricultural experiments. essay on principles. section 9.  _Statistical Science_ , pages 465–472, 1990.
  * Srebro et al. (2012) Nathan Srebro, Karthik Sridharan, and Ambuj Tewari.  Optimistic rates for learning with a smooth loss, 2012.
  * Sugiyama et al. (2007) Masashi Sugiyama, Matthias Krauledat, and Klaus-Robert Müller.  Covariate shift adaptation by importance weighted cross validation.  _Journal of Machine Learning Research_ , 8(5), 2007.
  * Teng et al. (2021) Jiaye Teng, Zeren Tan, and Yang Yuan.  T-sci: A two-stage conformal inference algorithm with guaranteed coverage for cox-mlp.  In _International Conference on Machine Learning_ , pages 10203–10213. PMLR, 2021.
  * Tibshirani et al. (2019) Ryan J Tibshirani, Rina Foygel Barber, Emmanuel Candes, and Aaditya Ramdas.  Conformal prediction under covariate shift.  _Advances in neural information processing systems_ , 32, 2019.
  * Tsiatis (2014) Anastasios Tsiatis.  Sensitivity analysis: A semi-parametric perspective.  In _Handbook of Missing Data Methodology_ , pages 403–428. Chapman and Hall/CRC, 2014.
  * Van der Vaart (2000) Aad W Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * van der Vaart and Wellner (1996) Aad W van der Vaart and Jon A Wellner.  _Weak convergence and empirical processes_.  Springer, 1996.
  * van der Vaart (2002) AW van der Vaart.  Semiparametric statistics, ecole d’ete de saint-flour 1999. in “lectures on probability theory and statistics”, 2002.
  * Vovk (2012) Vladimir Vovk.  Conditional validity of inductive conformal predictors.  In _Asian conference on machine learning_ , pages 475–490, 2012.
  * Vovk et al. (2005) Vladimir Vovk, Alex Gammerman, and Glenn Shafer.  _Algorithmic learning in a random world_.  Springer Science & Business Media, 2005.
  * Yang and Kuchibhotla (2021) Yachong Yang and Arun Kumar Kuchibhotla.  Finite-sample efficient conformal prediction.  _arXiv preprint arXiv:2104.13871_ , 2021.
  * Yin et al. (2021) Mingzhang Yin, Claudia Shi, Yixin Wang, and David M. Blei.  Conformal sensitivity analysis for individual treatment effects, 2021\.
  * Zaffran et al. (2022) Margaux Zaffran, Aymeric Dieuleveut, Olivier Féron, Yannig Goude, and Julie Josse.  Adaptive conformal predictions for time series.  _arXiv preprint arXiv:2202.07282_ , 2022.
  * Zhang et al. (2019) Anru Zhang, Lawrence D Brown, and T Tony Cai.  Semi-supervised inference: General theory and estimation of means.  _The Annals of Statistics_ , 47(5):2538–2566, 2019.
  * Zhu and Goldberg (2009) Xiaojin Zhu and Andrew B Goldberg.  Introduction to semi-supervised learning.  _Synthesis lectures on artificial intelligence and machine learning_ , 3(1):1–130, 2009.

Supplement to “Doubly robust calibration of prediction sets under covariate
shift”

##  S.1 A review of conformal inference with non-exchangeable data

Politis (2015) considered a transformation based approach that transforms
non-i.i.d. data to i.i.d. data and applies i.i.d. data prediction sets on the
transformations. For example, in an AR(1) model, transform the observed data
to obtain estimated innovations which are assumed i.i.d. and predict the
future innovation, then observed data forecast is obtained using the AR(1)
model estimates and the prediction of future innovation. Chernozhukov et al.
(2018) dealt with time series data where they developed a randomization method
by including block structure in the permutation scheme and showed asymptotic
validity under some modeling assumptions on the conformity score when
exchangeability fails, see also Chernozhukov et al. (2021a) and Chernozhukov
et al. (2021b) from the same authors under more general settings. Cattaneo et
al. (2021) recently described a different approach to obtain prediction
regions with time series data in a synthetic control framework. Recently,
Oliveira et al. (2022) proved that the split conformal prediction methodology
retains asymptotic coverage guarantee for several dependent data settings.

In a separate strand of work, Gibbs and Candès (2021) developed an adaptive
approach based on ideas from conformal inference that builds predictions sets
in an online setting where the data generating distribution is allowed to vary
over time and established coverage validity over the long term; see also
follow-up work in Zaffran et al. (2022) for an adaptive, tuning-free method.
These works partly combine ideas from online learning and sequential
prediction literature.

Candès et al. (2021) developed a method for survival analysis subject to
administrative censoring that has approximate marginal coverage if the
censoring mechanism or the conditional survival function is estimated well,
and Teng et al. (2021) focused on a similar censoring scenario under a Cox
proportional hazards model under the strong ignorability condition. A recent
paper Barber et al. (2022) designed a new technique for non-exchangeable data
that does not treat data points symmetrically and is robust against
distribution shift. In almost all of these cases, the coverage guarantee is
attained asymptotically as the number of training samples diverges to
infinity.

##  S.2 Connection to Causal Inference

In this section, we briefly discuss how the goal of prediction in the
covariate shift setting connects with that of prediction of counterfactuals
and individual treatment effects in a potential outcome framework (Rubin,
1974; Splawa-Neyman et al., 1990). Lei and Candès (2020) were the first to
formally draw this connection. Given $N$ subjects, let $A_{i}\in\\{0,1\\}$
denote a binary treatment indicator, $(Y_{i}(1),Y_{i}(0))$ be the pair of
potential outcomes for unit $i$, and $X_{i}$ be the corresponding vector of
measured covariates needed to control of confounding. We assume that

$\bigl{(}Y_{i}(1),Y_{i}(0),A_{i},X_{i}\bigr{)}\stackrel{{\scriptstyle\text{
i.i.d. }}}{{\sim}}(Y(1),Y(0),A,X),$

where $(Y(1),Y(0),A,X)$ is a random vector. Under the stable unit treatment
value assumption (SUTVA) commonly assumed in the literature (see e.g. Rubin
(1990)), the observed dataset consists of triples
$(Y_{i}^{\text{obs}},A_{i},X_{i})$ where

$Y_{i}^{\mathrm{obs}}=\begin{cases}Y_{i}(1),&A_{i}=1,\\\
Y_{i}(0),&A_{i}=0.\end{cases}$

For regularity, we assume the distributions of $X|A=1$ and $X|A=0$ are
absolutely continuous with respect to each other. In this counterfactual
setting, we wish to predict the individual treatment effect (ITE)
$\tau_{i}:=Y_{i}(1)-Y_{i}(0)$, which cannot be observed because for every unit
only one potential outcome is observed while the other is missing. In order to
obtain such prediction, we make the standard unconfoundedness assumption (also
known as strong ignorability condition) that $(Y(1),Y(0))\perp A\,|\,X$, the
counterfactual analog of MAR. Our approach thus yields prediction intervals
for ITE with valid coverage for subjects in the study, for whom one potential
outcome is observed; below, we briefly also discuss how one might obtain
prediction intervals for subjects not in the study, for whom both potential
outcomes are missing, but a covariate $X$ is available.

For any treated unit $i$ in the study, i.e. with $A_{i}=1$, we construct a
prediction interval $\widehat{C}^{\mathrm{ITE}}_{i}$ for $\tau_{i}$ such that
$\widehat{C}^{\mathrm{ITE}}_{i}=Y_{i}^{\mathrm{obs}}-\widehat{C}_{0}(X_{i})$,
where $\widehat{C}_{0}(x)$ satisfies

$\mathbb{P}\bigl{(}Y(0)\in\widehat{C}_{0}(X)\mid A=1\bigr{)}\geq 1-\alpha;$
(E.1)

Similarly, for any untreated unit, we construct
$\widehat{C}^{\mathrm{ITE}}_{i}=\widehat{C}_{1}(X_{i})-Y_{i}^{\mathrm{obs}}$,
where $\widehat{C}_{1}(x)$ satisfies

$\mathbb{P}\bigl{(}Y(1)\in\widehat{C}_{1}(X)\mid A=0\bigr{)}\geq 1-\alpha.$
(E.2)

Thus, such construction has a guaranteed coverage for $\tau_{i}$ because

$\displaystyle\mathbb{P}\bigl{(}Y_{i}(1)-Y_{i}(0)\in\widehat{C}^{\mathrm{ITE}}_{i}\bigr{)}$
$\displaystyle=\mathbb{P}(A_{i}=1)\mathbb{P}(Y_{i}(0)\in\widehat{C}_{0}(X_{i})|A_{i}=1)+\mathbb{P}(A_{i}=0)\mathbb{P}(Y_{i}(1)\in\widehat{C}_{1}(X_{i})|A_{i}=0)$
$\displaystyle\geq(1-\alpha)\bigl{(}\mathbb{P}(A_{i}=1)+\mathbb{P}(A_{i}=0)\bigr{)}=1-\alpha,\quad\text{
if }\eqref{eq:treated}\text{ and }\eqref{eq:untreated}\text{ both hold.}$

A more general goal than (E.1) would be to build a prediction set
$\widehat{C}_{0}$ such that

$\mathbb{P}_{(X,Y(0))\sim Q_{X}\times P_{Y(0)\mid
X}}\left(Y(0)\in\widehat{C}_{0}(X)\right)\geq 1-\alpha,$ (E.3)

where $Q_{X}$ is some distribution for $X$. Note that (E.1) can be seen as a
special case of (E.3) with $Q_{X}=P_{X|A=1}$. Based on the untreated samples,
we learn the distribution of $Y(0)\mid X,A=0$. And by the unconfoundedness
assumption, this has the same distribution as $Y(0)\mid X,A=1$, and also
$Y(0)\mid X$.

A more challenging goal potentially of interest in many settings, might be to
obtain prediction sets for future subjects not in the study, for whom neither
$Y_{i}(1)$ nor $Y_{i}(0)$ is observed, but $X_{i}$ is available. Without any
knowledge of the relationship between the distributions of $Y(1)$ and $Y(0)$,
a simple approach would be to obtain a pair of prediction intervals at level
$1-\alpha/2$, namely $(\widehat{Y}^{L}_{1}(x),\widehat{Y}^{R}_{1}(x))$ for
$Y(1)$ and $(\widehat{Y}^{L}_{0}(x),\widehat{Y}^{R}_{0}(x))$ for $Y(0)$. Then
taking the difference of the two sets such that
$\widehat{C}^{\mathrm{ITE}}=\bigl{(}\widehat{Y}^{L}_{1}(x)-\widehat{Y}^{R}_{0}(x),\widehat{Y}^{R}_{1}(x)-\widehat{Y}^{L}_{0}(x)\bigr{)}$
would in principle yield a valid $(1-\alpha)$–coverage for $Y(1)-Y(0)$. See,
e.g. Section 4.1 of Lei and Candès (2020). We also note that Fan and Park
(2007) has provided sharp bounds on the distribution of the treatment effect
which might aid in building a more precise prediction set.

##  S.3 Semiparametric theory and influence functions

In this section, we provide the definition of an influence function in the
literature of semi-parametrics theory and give the derivation of (16).

######  Definition 1.

Given a semiparametric model $\mathcal{F}$, a law ${F}^{*}$ in $\mathcal{F}$,
and a class $\mathcal{A}$ of reg. parametric submodels of $\mathcal{F}$, a
real valued functional

$\theta:\mathcal{F}\rightarrow\mathbb{R}$

is said to be a pathwise differentiable or regular parameter at F* wrt
$\mathcal{A}$ in model $\mathcal{F}$ iff there exists $\psi_{{F}^{*}}({x})$ in
$\mathcal{L}_{2}(\dot{F}^{*})$ such that for each submodel in $\mathcal{A}$,
say indexed by $t$ and with ${F}^{*}={F}_{t^{*}}$, and score, say
$S_{t}\left(t^{*}\right)=s_{t}\left(X;t^{*}\right)$ at $t^{*}$, it holds that

$\left.\frac{\partial}{\partial
t}\theta\left(F_{t}\right)\right|_{t=t^{*}}=\mathbb{E}_{F^{*}}\left[\psi_{F^{*}}(X)S_{t}\left(t^{*}\right)\right]$

$\psi_{{F}^{*}}(.)$ is called a gradient of $\theta$ at ${F}^{*}$ (wrt
$\left.\mathcal{A}\right)$. If, in addition, $\psi_{{F}^{*}}({X})$ has mean
zero under $\mathrm{F}^{*},\psi_{{F}^{*}}({X})$ is most commonly referred to
as an influence function of the functional $\theta$ at $F^{*}$.

In our case of finding the $(1-\alpha)$-quantile for $R|T=1$, let $t$ denote
the index for the parametric submodels, $\theta(t)$ be the desired
$(1-\alpha)$-quantile for $R|T=1$ with $\theta(t^{*}):=r_{\alpha}$. Let
$u(R;\theta):=\mathbbm{1}(R\leq\theta|T=1)-(1-\alpha)$ with $\theta(t^{*})$
satisfies that $\mathbb{E}_{t}[u(R;\theta(t))]=0$. Then,

$\displaystyle 0$ $\displaystyle=\frac{\partial}{\partial
t}\mathbb{E}_{t}\bigl{[}u(R;\theta(t))\bigr{]}\bigr{|}_{t=t^{*}}$ (E.4)
$\displaystyle=\frac{\partial}{\partial
t}\mathbb{E}_{t}\bigl{[}u(R;\theta(t^{*}))\bigr{]}\bigr{|}_{t=t^{*}}+\frac{\partial}{\partial\theta}\mathbb{E}_{t^{*}}[u(R;\theta)]|_{\theta=\theta(t^{*})}\frac{\partial\theta(t)}{\partial
t}\biggr{|}_{t=t^{*}}.$ (E.5)

From this we conclude that

$\displaystyle\frac{\partial\theta(t)}{\partial t}\,\biggr{|}_{t=t^{*}}$
$\displaystyle=-\biggl{\\{}\frac{\partial}{\partial\theta}\mathbb{E}_{t^{*}}\bigl{[}u(R;\theta)\bigr{]}\bigr{|}_{\theta=\theta(t^{*})}\biggr{\\}}^{-1}\frac{\partial}{\partial
t}\mathbb{E}_{t}\bigl{[}u(R;\theta(t^{*}))\bigr{]}\bigr{|}_{t=t^{*}}$ (E.6)
$\displaystyle\propto\frac{\partial}{\partial
t}\mathbb{E}_{t}\biggl{\\{}\mathbb{E}_{t}\Big{(}\mathbbm{1}\\{R\leq\theta\bigr{\\}}\mid
T=1,X\Bigr{)}\Bigr{|}T=1\biggr{\\}}\,\biggr{|}_{t^{*}}$ (E.7)
$\displaystyle=\frac{\partial}{\partial
t}\mathbb{E}_{t}\biggl{\\{}\mathbb{E}_{t}\Big{(}\mathbbm{1}\\{R\leq\theta\bigr{\\}}\mid
T=0,X\Bigr{)}\Bigr{|}T=1\biggr{\\}}\,\biggr{|}_{t^{*}}$ (E.8)
$\displaystyle=\frac{\partial}{\partial
t}\mathbb{E}_{t}\biggl{\\{}\mathbb{E}_{t^{*}}\Bigl{(}\mathbbm{1}\\{R\leq\theta\bigr{\\}}\mid
T=0,X\Bigr{)}\Bigm{|}T=1\biggr{\\}}\,\biggr{|}_{t^{*}}$ (E.9)
$\displaystyle\qquad+\frac{\partial}{\partial
t}\mathbb{E}_{t^{*}}\biggl{\\{}\mathbb{E}_{t}\Bigl{(}\mathbbm{1}\\{R\leq\theta\bigr{\\}}\mid
T=0,X\Bigr{)}\Bigm{|}T=1\biggr{\\}}\,\biggr{|}_{t^{*}}$ (E.10)
$\displaystyle=:\@slowromancap i@+\@slowromancap <EMAIL_ADDRESS>(E.11)

Let $S(Z):=S((1-T)Y,T,X)=\mathbbm{1}\\{T=0\\}S(Y|T=0,X)+S(T|X)+S(X)$ be the
score vector for all the observed data. We analyze the two terms separately.
For the first term, because $R\perp T|X$, we have that

$\displaystyle\@slowromancap i@$
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}E_{t^{*}}\bigl{(}\mathbbm{1}\\{R\leq\theta\bigr{\\}}\mid
T=0,X\bigr{)}S_{X|T=1}\Bigm{|}T=1\biggr{]}$ (E.12)
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}\frac{\mathbbm{1}\\{T=1\\}}{\mathbb{P}(T=1)}E_{t^{*}}\bigl{(}\mathbbm{1}\\{R\leq\theta\bigr{\\}}\mid
T=0,X\bigr{)}S_{X|T=1}\biggr{]}$ (E.13)
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}\frac{\mathbbm{1}\\{T=1\\}}{\mathbb{P}(T=1)}\mathbb{P}(R\leq\theta\mid
X)S_{X|T}\biggr{]}$ (E.14)
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}\frac{\mathbbm{1}\\{T=1\\}}{\mathbb{P}(T=1)}\big{\\{}\mathbb{P}(R\leq\theta\mid
X)-\mathbb{E}(\mathbb{P}(R\leq\theta\mid X)|T)\bigr{\\}}S_{X|T}\biggr{]}$
(E.15)
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}\frac{\mathbbm{1}\\{T=1\\}}{\mathbb{P}(T=1)}\big{\\{}\mathbb{P}(R\leq\theta\mid
X)-\mathbb{E}(\mathbb{P}(R\leq\theta\mid X)|T)\bigr{\\}}S_{X,T}\biggr{]}$
(E.16)
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}\frac{\mathbbm{1}\\{T=1\\}}{\mathbb{P}(T=1)}\big{\\{}\mathbb{P}(R\leq\theta\mid
X)-\mathbb{E}(\mathbb{P}(R\leq\theta\mid X)|T)\bigr{\\}}S_{Z}\biggr{]}$ (E.17)
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}\mathbbm{1}\\{T=1\\}\Bigl{\\{}\mathbb{E}_{t^{*}}\bigl{(}\mathbbm{1}\\{R\leq\theta\bigr{\\}}\mid
X\bigr{)}-(1-\alpha)\Bigr{\\}}S(Z)\biggr{]}/\mathbb{P}(T=1),$ (E.18)

For the second term, for brevity let $A:=\mathbbm{1}\\{R\leq\theta\\}$,

$\displaystyle\@slowromancap ii@$ $\displaystyle=\frac{\partial}{\partial
t}\mathbb{E}_{t^{*}}\biggl{\\{}\frac{\mathbbm{1}\\{T=1\\}}{\mathbb{P}(T=1)}\mathbb{E}_{t}\Bigl{(}\mathbbm{1}\\{R\leq\theta\bigr{\\}}\mid
T=0,X\Bigr{)}\biggr{\\}}\,\biggr{|}_{t^{*}}$
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{\\{}\frac{\mathbbm{1}\\{T=1\\}}{\mathbb{P}(T=1)}\mathbb{E}_{t^{*}}\Bigl{(}AS_{A|T=0,X}\mid
T=0,X\Bigr{)}\biggr{\\}}$
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{\\{}\frac{\mathbbm{1}\\{T=1\\}}{\mathbb{P}(T=1)}\mathbb{E}_{t^{*}}\Bigl{(}AS_{A|T=0,X}(A)\Bigr{|}T=0,X\Bigr{)}\biggr{\\}}\,$
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}\mathbb{E}_{t^{*}}\biggl{\\{}AS_{A|T=0,X}(A,T=0,X)\Bigr{|}T=0,X\biggr{\\}}\biggm{|}T=1\biggr{]}$
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}\mathbb{E}_{t^{*}}\biggl{\\{}\frac{A\mathbbm{1}\\{T=0\\}}{\mathbb{P}(T=0|X)}S_{A|T=0,X}(A,T=0,X)\Bigl{|}X\biggr{\\}}\biggm{|}T=1\biggr{]}$
$\displaystyle=\mathbb{E}_{t^{*}}\biggl{[}\mathbb{E}_{t^{*}}\biggl{\\{}\frac{A\mathbbm{1}\\{T=0\\}}{\mathbb{P}(T=0|X)}S_{A|T,X}(A,T,X)\Bigl{|}X\biggr{\\}}\biggm{|}T=1\biggr{]}$
$\displaystyle=\int\int\int\frac{A\mathbbm{1}\\{T=0\\}}{\mathbb{P}(T=0|X=x)}S_{A|T,x}(r,T,X)f(r|T,x)f(T|x)f(x|T=1)\mathrm{d}r\mathrm{d}x\mathrm{d}T$
$\displaystyle=\int\int\int\frac{A\mathbbm{1}\\{T=0\\}\mathbb{P}(T=1|X=x)}{\mathbb{P}(T=1)\mathbb{P}(T=0|X=x)}S_{A|T,X}(r,T,X)f(r|T,X)f(T|x)f(x)\mathrm{d}r\mathrm{d}x\mathrm{d}T$
$\displaystyle=\mathbb{E}\Bigl{[}A\mathbbm{1}\\{T=0\\}\pi^{\star}(X)S_{A|T,X}(r,T,X)\Bigr{]}/\mathbb{P}(T=1)$
$\displaystyle=\mathbb{E}\Bigl{[}(A-\mathbb{E}(A|T=0,X))\mathbbm{1}\\{T=0\\}\pi^{\star}(X)S_{A|T,X}(R,T,X)\Bigr{]}/\mathbb{P}(T=1)$
$\displaystyle=\mathbb{E}\Bigl{[}\bigl{(}\mathbbm{1}\\{R\leq\theta\\}-m^{\star}(t,X)\bigr{)}\mathbbm{1}\\{T=0\\}\pi^{\star}(X)\bigl{\\{}\mathbbm{1}\\{T=0\\}S_{A|T,X}(R,T,X)+S_{T,X}(T,X)\bigr{\\}}\Bigr{]}/\mathbb{P}(T=1)$
$\displaystyle=\mathbb{E}\Bigl{[}\bigl{(}\mathbbm{1}\\{R\leq\theta\\}-m^{\star}(t,X)\bigr{)}\mathbbm{1}\\{T=0\\}\pi^{\star}(X)S(Z)\Bigr{]}/\mathbb{P}(T=1).$

Combining the two terms together gives us

$\displaystyle\frac{\partial}{\partial t}\theta(t)\,\biggr{|}_{t^{*}}$
$\displaystyle\propto\mathbb{E}_{t^{*}}\biggl{[}\Bigl{\\{}\mathbbm{1}\\{T=1\\}\Bigl{(}m^{\star}(\theta,X)-(1-\alpha)\Bigr{)}$
(E.19)
$\displaystyle\,\,+\mathbbm{1}\\{T=0\\}\pi^{\star}(x)\Bigl{(}\mathbbm{1}\\{R\leq\theta\\}-m^{\star}(\theta,X)\Bigr{)}\Bigr{\\}}S(Z))\biggr{]}.$
(E.20)

We therefore conclude that the following function $\phi(\cdot)$ is
proportional to a mean zero gradient of $\theta(t)$ at $t^{*}$,

$\displaystyle\phi(\theta,X,R,T;\pi^{\star},m^{\star})$
$\displaystyle=\mathbbm{1}\\{T=0\\}\pi^{\star}(X)\Big{[}\mathbbm{1}\\{R\leq\theta\\}-m^{\star}(\theta,X)\Big{]}+\mathbbm{1}\\{T=1\\}\Big{[}m^{\star}(\theta,X)-(1-\alpha)\Big{]}.$
(E.21)

##  S.4 Proof of Lemma 2

###### Proof of Lemma 2.

Throughout the proof, we will write $R$ instead of $R(X,Y)$ for convenience.
Firstly, note that

$\begin{split}\mathbb{P}_{(X,Y)\sim Q_{X}\times P_{Y\mid
X}}(R(X,Y)\leq\theta\big{|}\theta)&=\mathbb{E}_{X\sim
Q_{X}}\bigl{[}\mathbb{P}(R\leq\theta\mid X,\theta)\big{|}\theta\bigr{]}\\\
&=\mathbb{E}_{X\sim Q_{X}}\bigl{[}\mathbb{P}(R\leq\theta\mid
X,T=1,\theta)\big{|}\theta\bigr{]}\\\
&=\int_{\chi}\int_{-\infty}^{\theta}p(x|T=1)p(r|x,T=1)\mathrm{d}r\mathrm{d}x\\\
&=\int_{\chi}\int_{-\infty}^{\theta}p(x|T=1)p(r|x,T=0)\mathrm{d}r\mathrm{d}x\text{
(Using~{}\eqref{eq:conditional-dist-same})}\\\
&=\int_{\chi}\int_{-\infty}^{\theta}\frac{p(x|T=1)}{p(x|T=0)}p(r|x,T=0)p(x|T=0)\mathrm{d}r\mathrm{d}x\\\
&=\mathbb{E}\biggl{[}\frac{p(X|T=1)}{p(X|T=0)}\mathbbm{1}\\{R\leq\theta\\}\,\Bigr{|}\,T=0,\theta\biggr{]}.\end{split}$
(E.22)

On the other hand, we can prove that

$P[\mathrm{IF}(\theta,X,R,T;\pi,m)]=\mathbb{P}(T=1)\biggl{\\{}\mathbb{E}\biggl{[}\frac{p(X|T=1)}{p(X|T=0)}\mathbbm{1}\\{R\leq\theta\\}\,\Bigr{|}\,T=0,\theta\biggr{]}-(1-\alpha)\biggr{\\}}.$
(E.23)

This combined with (E.22) completes the proof of (19). We will prove (E.23) in
two steps.

$P[\mathrm{IF}(\theta,X,R,T;\pi,m)]=\mathbb{E}\Bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}\mathbb{P}(R\leq\theta|X,\theta)-(1-\alpha)\bigr{\\}}\big{|}\theta\Bigr{]},$
(Step 1)

and

$\begin{split}&\mathbb{E}\bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}\mathbb{P}(R\leq\theta|X,\theta)-(1-\alpha)\bigr{\\}}\big{|}\theta\bigr{]}\\\
&\quad=\mathbb{P}(T=1)\biggl{\\{}\mathbb{E}\biggl{[}\frac{p(X|T=1)}{p(X|T=0)}\mathbbm{1}\\{R\leq\theta\\}\,\Bigr{|}\,T=0,\theta\biggr{]}-(1-\alpha)\biggr{\\}}.\end{split}$
(Step 2)

In the proof of (Step 1), we will use the fact that either $\pi(\cdot)$ or
$m(\cdot,\cdot)$ represents the correct density ratio or the correct
conditional distribution function. The proof of (Step 2) follows essentially
from Bayes rule.

#### Proof of (Step 1).

If $\pi(x)=\mathbb{P}(T=1|X=x)/\mathbb{P}(T=0|X=x)$ for all $x$ (i.e., density
ratio is correct), then we have

$\mathbb{E}\left[\mathbbm{1}\\{T=0\\}\pi(X)\big{|}X=x,R\right]=\mathbb{P}(T=0|X=x,R)\frac{\mathbb{P}(T=1|X=x)}{\mathbb{P}(T=0|X=x)}=\mathbb{P}(T=1|X=x).$

The second equality here follows because $T$ is independent of $R$ given $X$.
This implies that

$\displaystyle\mathbb{E}\left[\mathbbm{1}\\{T=0\\}\pi(X)\big{\\{}\mathbbm{1}\\{R\leq\theta\\}-m(\theta,X)\big{\\}}\big{|}\theta\right]$
$\displaystyle\quad=\mathbb{E}\Bigl{[}\mathbb{P}(T=1|X)\big{\\{}\mathbbm{1}\\{R\leq\theta\\}-m(\theta,X)\big{\\}}\big{|}\theta\Bigr{]}$
$\displaystyle\quad=\mathbb{E}\Bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}\mathbb{P}(R\leq\theta|X,\theta)-m(\theta,X)\bigr{\\}}\big{|}\theta\Bigr{]}.$

Similarly,

$\mathbb{E}\Bigl{[}\mathbbm{1}\\{T=1\\}\left\\{m(\theta,X)-(1-\alpha)\right\\}\big{|}\theta\Bigr{]}=\mathbb{E}\Bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}m(\theta,X)-(1-\alpha)\bigr{\\}}\big{|}\theta\Bigr{]}.$

Hence, if $\pi(\cdot)$ is the true density ratio, then

$\displaystyle P[\mathrm{IF}(\theta,X,R,T;\pi,m)]$
$\displaystyle=\mathbb{E}\Bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}\mathbb{P}(R\leq\theta|X,\theta)-(1-\alpha)\bigr{\\}}\big{|}\theta\Bigr{]}.$

This completes the proof of (Step 1) when $\pi(\cdot)$ is the true density
ratio.

If $m(\gamma,x)=\mathbb{E}[\mathbbm{1}\\{R\leq\gamma\\}|X=x]$ for all
$\gamma\in\mathbb{R}$,
$x\in{\mathchoice{\raisebox{0.0pt}{$\displaystyle\chi$}}{\raisebox{0.0pt}{$\textstyle\chi$}}{\raisebox{0.0pt}{$\scriptstyle\chi$}}{\raisebox{0.0pt}{$\scriptscriptstyle\chi$}}}$
(i.e., the conditional mean is correct), then using the conditional
independence of $R$ and $T$ given $X$, we have

$\mathbb{E}\Bigl{[}\mathbbm{1}\\{T=0\\}\pi(X)\bigl{\\{}\mathbbm{1}\\{R\leq\theta\\}-m(\theta,X)\bigr{\\}}\big{|}\theta\Bigr{]}=0.$

Hence,

$\displaystyle P[\mathrm{IF}(\theta,X,R,T;\pi,m)]$
$\displaystyle=\mathbb{E}\Bigl{[}\mathbbm{1}\\{T=1\\}\bigl{\\{}m(\theta,X)-(1-\alpha)\bigr{\\}}\big{|}\theta\Bigr{]}$
$\displaystyle=\mathbb{E}\Bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}m(\theta,X)-(1-\alpha)\bigr{\\}}\big{|}\theta\Bigr{]}$
$\displaystyle=\mathbb{E}\Bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}\mathbb{P}(R\leq\theta|X,\theta)-(1-\alpha)\bigr{\\}}\big{|}\theta\Bigr{]}.$

This completes the proof of (Step 1) if $m(\cdot,\cdot)$ is the true
conditional mean function.

#### Proof of (Step 2).

$\begin{split}&\mathbb{E}\bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}\mathbb{P}(R\leq\theta|X,\theta)-(1-\alpha)\bigr{\\}}\big{|}\theta\bigr{]}\\\
&\quad=\mathbb{E}\bigl{[}\mathbb{P}(T=1|X)\mathbb{P}(R\leq\theta|X,\theta)\big{|}\theta\bigr{]}-\mathbb{P}(T=1)(1-\alpha)\\\
&\quad=\mathbb{E}\biggl{[}\mathbbm{1}\\{T=0\\}\frac{\mathbb{P}(T=1|X)}{\mathbb{P}(T=0|X)}\mathbb{P}(R\leq\theta|X,\theta)\big{|}\theta\biggr{]}-\mathbb{P}(T=1)(1-\alpha)\\\
&\quad=\mathbb{E}\biggl{[}\mathbbm{1}\\{T=0\\}\frac{\mathbb{P}(T=1|X)}{\mathbb{P}(T=0|X)}\mathbbm{1}\\{R\leq\theta\\}\big{|}\theta\biggr{]}-\mathbb{P}(T=1)(1-\alpha)\\\
&\quad\stackrel{{\scriptstyle(b)}}{{=}}\frac{\mathbb{P}(T=1)}{\mathbb{P}(T=0)}\mathbb{E}\biggl{[}\mathbbm{1}\\{T=0\\}\frac{p(X|T=1)}{p(X|T=0)}\mathbbm{1}\\{R\leq\theta\\}\bigl{|}\theta\biggr{]}-\mathbb{P}(T=1)(1-\alpha)\\\
&\quad=\frac{\mathbb{P}(T=1)}{\mathbb{P}(T=0)}\mathbb{E}_{T}\biggl{\\{}\mathbb{E}\biggl{[}\mathbbm{1}\\{T=0\\}\frac{p(X|T=1)}{p(X|T=0)}\mathbbm{1}\\{R\leq\theta\\}\Big{|}T,\theta\biggr{]}\biggr{\\}}-\mathbb{P}(T=1)(1-\alpha)\\\
&\quad=\frac{\mathbb{P}(T=1)}{\mathbb{P}(T=0)}\mathbb{P}(T=0)\mathbb{E}\biggl{[}\frac{p(X|T=1)}{p(X|T=0)}\mathbbm{1}\\{R\leq\theta\\}\Big{|}T=0,\theta\biggr{]}-\mathbb{P}(T=1)(1-\alpha)\\\
&\quad=\mathbb{P}(T=1)\biggl{\\{}\mathbb{E}\biggl{[}\frac{p(X|T=1)}{p(X|T=0)}\mathbbm{1}\\{R\leq\theta\\}\,\Bigr{|}\,T=0,\theta\biggr{]}-(1-\alpha)\biggr{\\}},\end{split}$
(E.24)

where equality (b) comes from Bayes rule. This completes the proof of (Step
2). ∎

##  S.5 Proof of Theorem 2

###### Proof.

By definition of the $\mathrm{IF}$ function, it holds
$\forall\gamma\in\mathbb{R}$,

$\displaystyle
P\bigl{[}\mathrm{IF}(\gamma,X,R,{\widehat{\pi}},\widehat{m})\bigr{]}$
$\displaystyle=P\Bigl{[}\mathbb{P}(T=0|X){\widehat{\pi}}(X)\bigr{\\{}\mathbb{P}(R\leq\gamma|X)-\widehat{m}(\gamma,X)\bigr{\\}}$
(E.25)
$\displaystyle\qquad+\mathbb{P}(T=1|X)\bigl{\\{}\widehat{m}(\gamma,X)-(1-\alpha)\bigr{\\}}\Bigr{]}$
$\displaystyle=P\Bigl{[}\mathbb{P}(T=0|X)\bigl{\\{}\widehat{\pi}(X)-\pi^{\star}(X)\bigr{\\}}\bigr{\\{}m^{\star}(\gamma,X)-\widehat{m}(\gamma,X)\bigr{\\}}\Bigr{]}$
$\displaystyle\qquad+P\Bigl{[}\mathbb{P}(T=0|X)\pi^{\star}(X)\bigl{\\{}m^{\star}(\gamma,X)-\widehat{m}(\gamma,X))\bigr{\\}}$
$\displaystyle\qquad+\mathbb{P}(T=1|X)\bigl{\\{}\widehat{m}(\gamma,X)-(1-\alpha)\bigr{\\}}\Bigr{]}$
$\displaystyle=P\Bigl{[}\mathbb{P}(T=0|X)\bigl{\\{}\widehat{\pi}(X)-\pi^{\star}(X)\bigr{\\}}\bigr{\\{}m^{\star}(\gamma,X)-\widehat{m}(\gamma,X)\bigr{\\}}\Bigr{]}$
$\displaystyle\qquad+P\Bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}m^{\star}(\gamma,X)-\widehat{m}(\gamma,X))\bigr{\\}}$
$\displaystyle\qquad+\mathbb{P}(T=1|X)\bigl{\\{}\widehat{m}(\gamma,X)-(1-\alpha)\bigr{\\}}\Bigr{]}$
$\displaystyle=P\Bigl{[}\mathbb{P}(T=0|X)\bigl{\\{}\widehat{\pi}(X)-\pi^{\star}(X)\bigr{\\}}\bigr{\\{}m^{\star}(\gamma,X)-\widehat{m}(\gamma,X)\bigr{\\}}\Bigr{]}$
$\displaystyle\qquad+P\Bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}m^{\star}(\gamma,X)-(1-\alpha)\bigr{\\}}\Bigr{]}.$

Repeating the same calculation in (E.25) with either $\widehat{\pi}$ replaced
by $\pi^{\star}$ and $\widehat{m}$ replaced by $m^{\star}$ yields
$\forall\gamma\in\mathbb{R}$,

$\displaystyle P\bigl{[}\mathrm{IF}(\gamma,X,R,\pi,m)\bigr{]}$
$\displaystyle=P\Bigl{[}\mathbb{P}(T=1|X)\bigl{\\{}m^{\star}(\gamma,X)-(1-\alpha)\bigr{\\}}\Bigr{]}.$
(E.26)

Therefore,

$\displaystyle\sup_{\gamma\in\mathbb{R}}\Bigl{|}P\bigl{[}\mathrm{IF}(\gamma,X,R,\widehat{\pi},\widehat{m})-\mathrm{IF}(\gamma,X,R,{\pi^{\star}},{m^{\star}})\bigr{]}\Bigr{|}$
(E.27)
$\displaystyle\quad=\sup_{\gamma\in\mathbb{R}}\Bigl{|}P\bigl{[}\mathbb{P}(T=0|X)\bigl{\\{}\pi^{\star}(X)-\widehat{\pi}(X)\bigr{\\}}\bigr{\\{}m^{\star}(\gamma,X)-\widehat{m}(\gamma,X)\bigr{\\}}\bigr{]}\Bigr{|}$
$\displaystyle\quad\leq\|\widehat{\pi}-\pi^{\star}\|_{2}\sup_{\gamma}\|\widehat{m}(\gamma,\cdot)-m^{\star}(\gamma,\cdot)\|_{2}.$

The last inequality here follows from Cauchy–Schwarz inequality. ∎

##  S.6 Proof of Theorem 3

###### Proof.

Without loss of generality assume the indexes in $\mathcal{I}_{2}$ is
$1,\dots,n$ with $n:=|\mathcal{I}_{2}|$, and we expand
$\mathbb{P}_{\mathcal{I}_{2}}\bigl{[}\mathrm{IF}(\theta,\mathcal{D}_{2};\widehat{\pi},\widehat{m})\bigr{]}-P\bigl{[}\mathrm{IF}(\theta,X,R,T;\widehat{\pi},\widehat{m})\bigr{]}$
into three parts,

$\displaystyle\mathbb{P}_{\mathcal{I}_{2}}\bigl{[}\mathrm{IF}(\theta,\mathcal{D}_{2};\widehat{\pi},\widehat{m})\bigr{]}-P\bigl{[}\mathrm{IF}(\theta,X,R,T;\widehat{\pi},\widehat{m})\bigr{]}$
(E.28) $\displaystyle=\;$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}\\{t_{i}=0\\}\widehat{\pi}(x_{i})\mathbbm{1}\\{r_{i}\leq\theta\\}-P\bigl{[}\mathbbm{1}\\{t_{i}=0\\}\widehat{\pi}(x_{i})\mathbbm{1}\\{r_{i}\leq\theta\\}\bigr{]}$
(E.29)
$\displaystyle\quad+\frac{1}{n}\sum_{i=1}^{n}\widehat{m}(\theta,x_{i})\bigl{(}\mathbbm{1}\\{t_{i}=1\\}-\mathbbm{1}\\{t_{i}=0\\}\bigr{)}-P\bigl{[}\widehat{m}(\theta,x_{i})(\mathbbm{1}\\{t_{i}=1\\}-\mathbbm{1}\\{t_{i}=0\\})\bigr{]}$
(E.30)
$\displaystyle\quad-(1-\alpha)\biggl{[}\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}\\{t_{i}=1\\}-\mathbb{P}(T=1)\biggr{]}$
(E.31) $\displaystyle=$
$\displaystyle:\;\mathcal{R}_{1}(\theta)+\mathcal{R}_{2}(\theta)+\mathcal{R}_{3},$
(E.32)

where the three terms will be controlled separately. In particular,
$\sup_{\theta}|\mathcal{R}_{1}(\theta)|$ and
$\sup_{\theta}|\mathcal{R}_{2}(\theta)|$ will be bounded using tools from the
empirical processes theory. Also notice that conditional on training data
$\mathcal{D}_{1}$, $\widehat{\pi}$ and $\widehat{m}$ are non-random functions
and for ease of notation, we treat these functions as non-random and omit the
conditioning part from this point onwards.

For $W_{i}=(X_{i},Y_{i},T_{i}),i\in\mathcal{I}_{2}$ and any function
$f:\mathbb{R}^{d+2}\to\mathbb{R}$, for notation simplicity we define

$\mathbb{G}_{n}f~{}:=~{}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\bigl{\\{}f(W_{i})-\mathbb{E}[f(W_{i})]\bigr{\\}}.$

#### Bound on $\sup_{\theta}|\mathcal{R}_{1}(\theta)|$:

We have a class of functions

$\mathcal{F}=\bigl{\\{}f:f_{\widehat{\pi},\theta}(w)=\mathbbm{1}\\{t=1\\}\widehat{\pi}(x)\mathbbm{1}{\\{R(x,y)\leq\theta\\}},\forall\widehat{\pi}(\cdot)\in\mathcal{F}_{\pi}\bigr{\\}}.$

Notice that $\forall\theta\in\mathbb{R}$ and
$\widehat{\pi}(\cdot)\in\mathcal{F}_{\pi}$, we have
$|f_{\theta}(w)|\leq\pi_{0}\mathbbm{1}\\{t=1\\}$. Therefore,
$F(w):=\pi_{0}\mathbbm{1}\\{t=1\\}$ is an envelope function of
$\\{f_{\theta}(\cdot):\theta\in\mathbb{R}\\}$. Let $\|\cdot\|_{\mathcal{F}}$
denote the supremum norm $\|z\|_{\mathcal{F}}=\sup_{f\in\mathcal{F}}|z(f)|$.

Applying Lemma 4 with $s(t,x)=\mathbbm{1}\\{t=0\\}\widehat{\pi}(x)$ and
$h(x,y)=R(x,y)$ gives us