# Model-free selective inference under covariate shift
via weighted conformal p-values

Ying Jin  Department of Statistics, Stanford University  Emmanuel J. Candès
Department of Statistics, Stanford University  Department of Mathematics,
Stanford University

###### Abstract

This paper introduces novel weighted conformal p-values and methods for model-
free selective inference. The problem is as follows: given test units with
covariates $X$ and missing responses $Y$, how do we select units for which the
responses $Y$ are larger than user-specified values while controlling the
proportion of false positives? Can we achieve this without any modeling
assumptions on the data and without any restriction on the model for
predicting the responses? Last, methods should be applicable when there is a
covariate shift between training and test data, which commonly occurs in
practice.

We answer these questions by first leveraging any prediction model to produce
a class of well-calibrated weighted conformal p-values, which control the
type-I error in detecting a large response. These p-values cannot be passed
onto classical multiple testing procedures since they may not obey a well-
known positive dependence property. Hence, we introduce weighted conformalized
selection (WCS), a new procedure which controls false discovery rate (FDR) in
finite samples. Besides prediction-assisted candidate selection, WCS (1)
allows to infer multiple individual treatment effects, and (2) extends to
outlier detection with inlier distributions shifts. We demonstrate performance
via simulations and applications to causal inference, drug discovery, and
outlier detection datasets.

##  1 Introduction

Many scientific discovery and decision making tasks concern the selection of
promising candidates from a potentially large pool. In drug discovery,
scientists aim to find drug candidates—from a huge chemical space of molecules
or compounds—with sufficiently high binding affinity to a target. College
admission officers or hiring managers seek applicants with a high potential of
excellence. In healthcare, one would like to allocate resources to individuals
who would benefit the most. The list of such examples is endless. In all these
problems, a decision maker would like to identify those units for which an
unknown outcome of interest (affinity for a target molecule/job
performance/reduction in stress level) takes on a large value.

Machine learning has shown tremendous potential for accelerating these
resource- and time-consuming processes whereby predictions for the unknown
outcomes are used to shortlist promising candidates before any further in-
depth investigation [LDG13, CH18, SSKMM19]. This is a sharp departure from
traditional strategies that either directly evaluate the outcomes via a
physical screening to determine drug binding affinities, or “predict” the
unknown outcomes via human judgement in healthcare or job hiring applications.
In the new paradigm, a prediction model draws evidence from training data to
inform values the unknown outcome may take on given the features of a unit.

If the intent is thus to use machine learning to screen candidates, we would
need to understand and control the uncertainty in black-box prediction models
as a selection tool. First, gathering evidence for new discoveries from past
instances is challenging as they are not necessarily similar; i.e. they may be
drawn from distinct distributions. Second, since a subsequent investigation
applied to the shortlisted candidates is often resource intensive, it is
meaningful to limit the _error rate on the selected_. It is not so much the
accuracy of the prediction that matters, rather the proportion of selected
units not passing the test. Keeping such a proportion at sufficiently low
levels is a non-traditional criterion in typical prediction problems.

This paper studies the reliable selection of promising units/individuals whose
unknown outcomes are sufficiently large. We develop methods that apply in a
model-free fashion—without making any modeling assumptions on the data
generating process while controlling the expected proportion of false
discoveries among the selected.

###  1.1 Reliable selection under distribution shift

Formally, suppose test samples
$\\{(X_{j},Y_{j})\\}_{j\in\mathcal{D}_{\textnormal{test}}}$ are drawn i.i.d.
from an unknown distribution $Q$, whose outcomes
$\\{Y_{j}\\}_{j\in\mathcal{D}_{\textnormal{test}}}$ are unobserved. Our goal
is to identify a subset $\mathcal{R}\subseteq\mathcal{D}_{\textnormal{test}}$
with outcomes above some known, potentially random, thresholds
$\\{c_{j}\\}_{j\in\mathcal{D}_{\textnormal{test}}}$. To draw statistical
evidence for large outcomes in $\mathcal{D}_{\textnormal{test}}$, we assume we
hold an independent set of i.i.d. calibration data
$\\{(X_{i},Y_{i})\\}_{i\in\mathcal{D}_{\textnormal{calib}}}$ whose outcomes
are, this time, observed.

Before we begin our discussion, we note that a challenge that underlies almost
all applications is that samples in $\mathcal{D}_{\textnormal{calib}}$ are
often different from those in $\mathcal{D}_{\textnormal{test}}$. That is, in
many cases, a point $(X_{j},Y_{j})$ in $\mathcal{D}_{\textnormal{test}}$ is
sampled from $Q$ while a point $(X_{i},Y_{i})$ in
$\mathcal{D}_{\textnormal{calib}}$ is sampled from a different distribution
$P$. We elaborate on the distribution shift issue by considering two
motivating examples.

##### Drug discovery.

Virtual screening, which identifies viable drug candidates using predicted
affinity scores from supervised machine learning models, is increasingly
popular in early stages of drug discovery [Hua07, KMLH17, VCC+19, DDJ+21]. In
this application, $\mathcal{D}_{\textnormal{calib}}$ is the set of drugs with
known affinities, while $\mathcal{D}_{\textnormal{test}}$ is the set of new
drugs whose binding affinities are of interest. In practice, there usually is
a distribution shift between the two batches: for instance, an experimenter
may favor a specific structure when choosing candidates to physically screen
their true binding affinities [PPS+15]. Since physically screened drugs are
subsequently added to the calibration dataset
$\mathcal{D}_{\textnormal{calib}}$, they may not be representative of the new
drugs in $\mathcal{D}_{\textnormal{test}}$. Accounting for such distribution
shifts when applying virtual screening is crucial for the reliability of the
whole drug discovery process, and has attracted recent attention [Krs21].

##### College admission.

A college may be interested in prospective students who are likely to graduate
after four years of undergraduate education (a binary outcome of interest)
among all applicants in $\mathcal{D}_{\textnormal{test}}$ as to maintain a
reasonable graduation rate. The issue is that
$\mathcal{D}_{\textnormal{calib}}$ usually consists of students who _were_
admitted in the past as the outcome is only known for these students. Since
previous cohorts may have been selected by applying various admission
criteria, the distribution of the training data is different from that of the
new applicants even if the distribution of applicants does not much vary over
the years. Similar concerns apply to job hiring when using prediction models
to select suitable applicants [FRTT12, SS16]; the way an institution chooses
to record the outcomes of former candidates can lead to a discrepancy between
the documented candidates and those under consideration.

In this paper, we propose to address the distribution shift issue by
considering a scenario where the shift can be entirely attributed to observed
covariates; this is commonly referred to as a covariate shift in the
literature [SKM07, TFBCR19]. Concretely, we assume that there is some function
$w\colon\mathcal{X}\to\mathbb{R}^{+}$ such that

$\displaystyle\frac{{\,\mathrm{d}}Q}{{\,\mathrm{d}}P}(x,y)=w(x),\quad
P\textrm{-almost surely},$ (1)

where ${\,\mathrm{d}}Q/{\,\mathrm{d}}P$ denotes the Radon-Nikodym derivative.
While all kinds of distribution shift could happen in practice, covariate
shift—widely used to characterize distribution shifts caused by selection on
known features—covers many important cases. In drug discovery, if the
preference for choosing which samples to screen only depends upon the
covariates, e.g. the sequence of amino acids of the compound or the physical
properties of the compound, we are dealing with the covariate shift (1). In
college admission or job hiring, if a specific type of applicant—characterized
by certain features—were favored in previous admission/hiring cycles, such a
preference would yield the shift (1) between previous students/employees and
current applicants.

Returning to our decision problem, we are interested in finding as many units
$j\in\mathcal{D}_{\textnormal{test}}$ as possible with $Y_{j}>c_{j}$, while
controlling the false discovery rate (FDR) below some pre-specified level
$q\in(0,1)$. The FDR is the expected fraction of false discoveries, where a
false discovery is a selected unit for which $Y_{j}\leq c_{j}$, i.e., the
outcome turns out to fall below the threshold. Formally, we wish to find
$\mathcal{R}$ obeying

$\displaystyle\textnormal{FDR}:=\mathbb{E}\Bigg{[}\frac{\sum_{j\in\mathcal{D}_{\textnormal{test}}}\operatorname{\mathds{1}}\\{j\in\mathcal{R},Y_{j}\leq
c_{j}\\}}{\max\\{1,|\mathcal{R}|\\}}\Bigg{]}\leq q,$ (2)

where $q$ is the nominal target FDR level. The expectation in (2) is taken
over both the new test samples and the calibration data/screening procedure.
The FDR quantifies the tradeoff between resources dedicated to the shortlisted
candidates and the benefits from true positives [BH95]. By controlling the
FDR, we ensure that a sufficiently large proportion of follow-up
resources—such as clinical trials in drug discovery, human evaluation of
student profiles and interviews—to confirm predictions from machine learning
models are devoted to interesting candidates.

The problem of FDR-controlled selection with machine learning predictions has
been studied in an earlier paper [JC22] under the assumption of no
distribution shift. This means the new instances are exchangeable with the
known calibration data, which, as we discussed above, is rarely the case in
practice. Therefore, there is a pressing need for novel techniques to address
the distribution shift issue.

###  1.2 Preview of theoretical contributions

Our strategy consists in (1) constructing a calibrated confidence measure
(p-value), which applies to any predictive model, for detecting a large
unknown outcome, and (2) in employing multiple testing ideas to build the
selection set from test p-values.

##### Weighted conformal p-value: a calibrated confidence measure.

We introduce random variables
$\\{p_{j}\\}_{j\in\mathcal{D}_{\textnormal{test}}}$, which resemble p-values
in the sense that for each $j\in\mathcal{D}_{\textnormal{test}}$, they obey

$\displaystyle\mathbb{P}(p_{j}\leq t,~{}Y_{j}\leq c_{j})\leq
t,\quad\textrm{for all}~{}t\in[0,1],$ (3)

where the probability is over both the calibration data and the test sample
$(X_{j},Y_{j},c_{j})$. This is similar to the definition of a classical
p-value, hence, we call $p_{j}$ a“p-value”. With this, for any
$\alpha\in(0,1)$, selecting $j$ if and only if $p_{j}\leq\alpha$ controls the
type-I error below $\alpha$ (the chance that $Y_{j}$ is below threshold is at
most $\alpha$). Note however that (3) is perhaps an unconventional notion of
type-I error, since it accounts for the randomness in the “null hypothesis”
$H_{j}\colon Y_{j}\leq c_{j}$.

We construct $p_{j}$ using ideas from conformal prediction [VGS05, TFBCR19].
Take any monotone function $V\colon\mathcal{X}\times\mathcal{Y}\to\mathbb{R}$
(see Definition 2.1) built upon any prediction model that is independent of
$\mathcal{D}_{\textnormal{calib}}$ and $\mathcal{D}_{\textnormal{test}}$. A
concrete instance is $V(x,y)=y-\widehat{\mu}(x)$, where
$\widehat{\mu}\colon\mathcal{X}\to\mathcal{Y}$ is any predictor trained on a
held-out dataset. With scores $V_{i}=V(X_{i},Y_{i})$,
$i\in\mathcal{D}_{\textnormal{calib}}$, and weight function $w(\cdot)$ given
by (1), we construct $p_{j}$ to be approximately equal to
$\widehat{F}(\widehat{V}_{j})$, where $\widehat{F}(\cdot)$ is the cumulative
distribution function (c.d.f.) of the empirical distribution

$\frac{\sum_{i\in\mathcal{D}_{\textnormal{calib}}}\,w(X_{i})\,\delta_{V_{i}}}{\sum_{i\in\mathcal{D}_{\textnormal{calib}}}w(X_{i})}.$

While the precise definition of $p_{j}$ is in Section 2, Figure 1 visualizes
the idea.

Intuitively, $p_{j}$ contrasts the value of $\widehat{V}_{j}=V(X_{j},c_{j})$
against the distribution of the unknown ‘oracle’ score $V_{j}=V(X_{j},Y_{j})$.
The latter is approximated by $\widehat{F}(\cdot)$, which reweights the
training data with weights reflecting the covariate shift (1). As such,
$\widehat{F}(V_{j})$ is approximately uniformly distributed in $[0,1]$. Thus,
a notably small value of $p_{j}\approx\widehat{F}(\widehat{V}_{j})$ relative
to the uniform distribution informs that $\widehat{V}_{j}$ is smaller than
what is typical of $V_{j}$, which further suggests evidence for $Y_{j}\geq
c_{j}$. Leveraging conformal prediction techniques, such high-level ideas are
made exact to achieve (3).

Figure 1: Visualization of weighted conformal p-values. Left: Distribution of
scores in the calibration set (gray) and individual values of $w(X_{i})$
(orange) for $i\in\mathcal{D}_{\textnormal{calib}}$. Right: Weighted empirical
c.d.f. of $\\{V_{i}\\}_{i\in\mathcal{D}_{\textnormal{calib}}}$; the red dashed
lines illustrate how $p_{j}$ is computed from $\widehat{V}_{j}$, in which
$j\in\mathcal{D}_{\textnormal{test}}$.

##### Intricate dependence among weighted conformal p-values.

Going beyond per-unit type-I error control, selecting from multiple test
samples requires dealing with multiple p-values
$\\{p_{j}\\}_{j\in\mathcal{D}_{\textnormal{test}}}$. A natural idea is then to
use multiple testing ideas as if $\\{p_{j}\\}$ were classical p-values.
However, the situation is complicated here because
$\\{p_{j}\\}_{j\in\mathcal{D}_{\textnormal{test}}}$ all depend on the same set
of calibration data, and, in sharp contrast to classical multiple testing,
there is another source of randomness in the unknown outcomes.

Our second theoretical result shows the mutual dependence among multiple
p-values is particularly challenging. It states that the favorable positive
dependence structure among p-values, which was the key to FDR control with
other p-values derived with conformal prediction ideas [BCL+21, JC22], can be
violated in the presence of data-dependent weights in (4). As a result, it
remains unclear whether applying existing multiple testing procedures controls
the FDR in finite samples.

##### Weighted Conformalized Selection (WCS).

Last, we introduce WCS, a new procedure that achieves finite-sample FDR
control under distribution shift. The idea is to calibrate the selection
threshold for each $j\in\mathcal{D}_{\textnormal{test}}$ using a set of
“auxiliary” p-values
$\\{\bar{p}_{\ell}^{(j)}\\}_{\ell\in\mathcal{D}_{\textnormal{test}}}$. These
p-values happen to obey a special conditional independence property, namely,

$\displaystyle\bar{p}_{j}^{(j)}\perp\\!\\!\\!\perp\\{\bar{p}_{\ell}^{(j)}\\}_{\ell\neq
j}\,\big{|}\,\mathcal{Z}_{j},$

where $\mathcal{Z}_{j}$ is the unordered set of the calibration data and the
$j$-th test sample. We show that FDR control is valid regardless of the
machine learning model used in constructing $V(\cdot,\cdot)$, and applies to a
wide range of scenarios where the cutoffs $c_{j}$ can themselves be random
variables.

###  1.3 Broader scope

Our methods happen to be applicable to additional statistical inference
problems.

##### Multiple counterfactual inference.

The individual treatment effect (ITE) is a random variable that characterizes
the difference between an individual’s outcome when receiving a treatment
versus not, whose inference relies on comparing the observed outcome with the
counterfactual [IR15, LC21, JRC23]. To infer the counterfactual of a unit
under treatment, the calibration data are units who do not receive the
treatment. If the treatment assignment depends on the covariates, then there
may be, and usually will be, a covariate shift between treated and control
units. We discuss this problem in Section 5.

##### Outlier detection.

Imagine we have access to i.i.d. inliers
$\\{Z_{i}\\}_{i\in\mathcal{D}_{\textnormal{calib}}}$ from an unknown
distribution $P$. Then outlier detection [BCL+21] seeks to find outliers among
new samples $\\{Z_{j}\\}_{j\in\mathcal{D}_{\textnormal{test}}}$, i.e., those
which do not follow the inlier distribution $P$. A potential application is
fraud detection in financial transactions, where controlling the FDR ensures
reasonable resource allocation in follow-up inquiries. However, financial
behavior varies with users and it is possible that normal transactions of
interest follow different distributions across different populations. A
related problem is to detect concept drifts; that is, allowing a subset of
features $X\subseteq Z$ to follow a distinct distribution, and only detecting
test samples whose $Z{\,|\,}X$ distribution differs from $P_{Z{\,|\,}X}$. We
discuss this problem in Section 6.

As in (3), the FDR (2) is marginalized over the randomness in the hypotheses.
To this end, WCS requires all $\\{X_{n+j},Y_{n+j}\\}_{j=1}^{m}$ to be i.i.d.
draws from a super-population. To relax this, we generalize our methods to
control a notion of FDR that is conditional on the random hypotheses. This is
useful for handling imbalanced data in classification, and the aforementioned
outlier detection problem.

###  1.4 Data and code

An R package, ConfSelect, that implements Weighted Conformalized Selection (as
well as the version without weights) together with code to reproduce all
experiments in this paper, can be found in the GitHub repository
https://github.com/ying531/conformal-selection.

###  1.5 Related work

Our methods build upon the conformal inference framework [VGS05] to infer
unknown outcomes. There, the theoretical guarantee for conformal prediction
sets is usually marginally valid over the randomness in a single test point.
As discussed in [JC22], one might be interested in multiple test samples
simultaneously; in such situations, standard conformal prediction tools are
insufficient.

This work is connected to a recent line of research on selective inference
issues arising in predictive inference. Some works use conformal p-values
(whose definition differs from ours) with ‘full’ observations
$\\{(X_{j},Y_{j})\\}_{j\in\mathcal{D}_{\textnormal{test}}}$—that is, the
response is observed—for outlier detection, i.e., for identifying whether
$(X_{j},Y_{j})$ follows the same distribution as the calibration data
$\\{(X_{i},Y_{i})\\}_{i\in\mathcal{D}_{\textnormal{calib}}}$ [BCL+21, LSS22,
MLMR22]. Whereas these methods do not consider distribution shift between
calibration and test inliers, Section 6 extends all this. Here, the response
$Y_{n+j}$—the object of inference—is not observed and, therefore, the problem
is completely different. Consequently, methods and techniques to deal with
this are new. Our focus on FDR control is similar in spirit to other recent
papers; for instance, rather than marginal coverage of prediction sets,
[FSJB22] studies the number of false positives, and [BHRZ23] the miscoverage
rate on selected units.

Our methodology draws ideas from the multiple hypothesis testing literature,
where the FDR is a popular type-I error criterion in exploratory analysis for
controlling the proportion of ‘false leads’ in follow-up studies [BH95, BY01].
This paper is distinct from the conventional setting. First, testing a random
outcome instead of a model parameter leads to a random null set and a
complicated dependence structure. Second, our inference relies on (weighted)
exchangeability of the data while imposing no assumptions on their
distributions. In contrast, a null hypothesis often specifies the distribution
of test statistics or p-values. Later on, we shall also draw connections to a
recent literature dealing with dependence in multiple testing, such as that
concerned with e-values [VW21, WR20] and conditional calibration of the BH
procedure [FL22].

The application of our method is related to a substantial literature that
emphasizes the importance of uncertainty quantification in drug discovery
[NCBE14, SNB17, AHBC17, SAN+18, CCB19, WJR22]. Although existing works aim to
employ conformal inference to assign confidence to machine prediction, this is
often done in heuristic ways, which potentially invalidates the guarantees of
conformal prediction due to selection issues. Instead, we provide here a
principled approach to prioritizing drug candidates with clear and
interpretable error control (i.e., guaranteeing a sufficiently high ‘hit rate’
[WJR22]), and offer a potential solution to the widely recognized distribution
shift problem [Krs21, FL23].

The application to individual treatment effects in Section 5 is connected to
[CDLM21]; the distinction is that they focus on a summary statistic such as a
population quantile of the ITEs, while our method detects individuals with
positive ITEs. Also related is [DWR21], which tests individuals for positive
treatment effects with FDR control; however, in their setting a positive ITE
means the whole distribution of the control outcome is dominated by that of
the treated outcome, whereas we directly compare two random variables.
Accordingly, our techniques are very different from these two references.

##  2 Weighted conformal p-values

###  2.1 Construction of weighted conformal p-values

Our weighted conformal p-values build upon any _monotone_ score function,
defined as follows.

######  Definition 2.1 (Monotonicity).

A score function
$V(\cdot,\cdot)\colon\mathcal{X}\times\mathcal{Y}\to\mathbb{R}$ is monotone if
$V(x,y)\leq V(x,y^{\prime})$ holds for any $x\in\mathcal{X}$ and any
$y,y^{\prime}\in\mathcal{Y}$ obeying $y\leq y^{\prime}$.

Intuitively, the score function (often referred to as nonconformity score in
conformal prediction [VGS05]) describes how well a hypothesized value
$y\in\mathcal{Y}$ _conforms_ to the machine prediction. A popular and monotone
nonconformity score is $V(x,y)=y-\widehat{\mu}(x)$, where
$\widehat{\mu}\colon\mathcal{X}\to\mathbb{R}$ is a point prediction; other
choices include those based on quantile regression [RPC19] or estimates of the
conditional c.d.f. [CWZ21].

From now on, write $\mathcal{D}_{\textnormal{calib}}=\\{1,\dots,n\\}$ and
$\mathcal{D}_{\textnormal{test}}=\\{n+1,\dots,n+m\\}$. Compute
$V_{i}=V(X_{i},Y_{i})$ for $i=1,\dots,n$ and
$\widehat{V}_{n+j}=V(X_{n+j},c_{n+j})$ for $j=1,\dots,m$. Our weighted
conformal p-values are defined as

$\displaystyle
p_{j}=\frac{\sum_{i=1}^{n}w(X_{i})\operatorname{\mathds{1}}{\\{V_{i}<\widehat{V}_{n+j}\\}}+(w(X_{n+j})+\sum_{i=1}^{n}w(X_{i})\operatorname{\mathds{1}}{\\{V_{i}=\widehat{V}_{n+j}\\}})\cdot
U_{j}}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})},$ (4)

where $w(\cdot)$ is the covariate shift function in (1), and
$U_{j}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}\textrm{Unif}([0,1])$
are tie-breaking random variables. When $w(\cdot)\equiv 1$, our p-values
reduce to the conformal p-values in [JC22], see Appendix A.1. We also refer
readers to Appendix A.2 for the connection between our p-values and conformal
prediction intervals.

Roughly speaking, with a monotone score, the p-value (4) measures how small
$c_{n+j}$ is compared with typical values of $Y_{n+j}$, and provides
calibrated evidence for a large outcome as expressed in (3). This can be
derived using conformal inference theory [TFBCR19] and the monotocity of $V$.
We include a formal result here for completeness, whose proof is in Appendix
B.1.

######  Lemma 2.2.

The tail inequality (3) holds if the covariate shift (1) holds and the score
function is monotone.

Consequently, selecting a unit if $p_{j}\leq\alpha$ controls the type-I error
at level $\alpha$ for each $j\in\mathcal{D}_{\textnormal{test}}$. We shall
however see that dealing with multiple test samples is far more challenging.

###  2.2 Weighted conformal p-values are not PRDS

Applying multiple testing procedures to our p-values to obtain a selection set
naturally comes to mind. Previous works applying the Benjamini-Hochberg (BH)
procedure [BH95] to p-values obtained via conformal inference ideas indeed
control the FDR [BCL+21, JC22, MLMR22]. This is because the p-values obey a
favorable dependence structure called _positive dependence on a subset_
(PRDS).

######  Definition 2.3 (PRDS).

A random vector $X=(X_{1},\dots,X_{m})$ is PRDS on a subset $\mathcal{I}$ if
for any $i\in\mathcal{I}$ and any increasing set $D$, the probability
$\mathbb{P}(X\in D{\,|\,}X_{i}=x)$ is increasing in $x$.

Above, a set $D\subseteq\mathbb{R}^{m}$ is _increasing_ if $a\in D$ and
$b\succeq a$ implies $b\in D$ ($b\succeq a$ means that all the components of
$b$ are larger than or equal to those of $a$). It is well-known that the BH
procedure controls the FDR when applied to PRDS p-values [BY01]. The PRDS
property among (unweighted) conformal p-values was first studied in
[BCL+21],111They also show the non-randomized unweighted conformal p-values
are PRDS when the scores are almost surely distinct. We do not elaborate on
this distinction as it is not the focus here. and generalized in [JC22] to
plug-in values $c_{n+j}$ in lieu of $Y_{n+j}$, thereby forming the basis for
FDR control. If the PRDS property were to hold for the weighted conformal
p-values, then applying the BH procedure would also control the FDR. This is
however not always the case; a constructive proof of Proposition 2.4 is in
Appendix B.3.

######  Proposition 2.4.

There exist a sample size $n$, a weight function
$w\colon\mathcal{X}\to\mathbb{R}^{+}$, a nonconformity score
$V\colon\mathcal{X}\times\mathcal{Y}\to\mathbb{R}$, and distributions $P$ and
$Q$ obeying (1), such that the weighted conformal p-values (4) are not PRDS
for
$\\{X_{i},Y_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}P$
and
$\\{X_{n+j},Y_{n+j}\\}_{j=1}^{m}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}Q$.

A little more concretely, we find that the PRDS property is likely to fail
when the nonconformity scores are negatively correlated with the weights. Note
that the dependence among conformal p-values arises from sharing the
calibration data. Under exchangeability (i.e. taking equal weightes
$w(x)\equiv 1$), a smaller conformal p-value is (only) associated with larger
calibration scores, and hence smaller p-values for other test points. However,
with negatively associated weights and scores, a smaller p-value may also be
due to smaller calibration weights—instead of larger calibration scores—which
implies that other p-values may take on larger values.

###  2.3 BH procedure controls FDR asymptotically

Before introducing our new selection methods, we nevertheless show asymptotic
FDR control using the BH procedure.

######  Theorem 2.5.

Suppose $w(\cdot)$ is uniformly bounded by a fixed constant, the covariate
shift (1) holds, and $\\{c_{n+j}\\}_{j=1}^{m}$ are i.i.d. random variables.
Fix any $q\in(0,1)$, and let $\mathcal{R}$ be the output of the BH procedure
applied to $\\{p_{j}\\}_{j=1}^{m}$ in (4). Then the following results hold.

  1. (i)

For any fixed $m$, it holds that $\limsup_{n\to\infty}\textnormal{FDR}\leq q$.

  2. (ii)

Suppose $m,n\to\infty$. Under a mild technical condition (see Appendix A.4),
$\limsup_{m,n\to\infty}\textnormal{FDR}\leq q$. Furthermore, the asymptotic
FDR and power of BH can be exactly characterized in this case.

We provide the complete technical version of Theorem 2.5 in Theorem A.3, whose
proof is in Appendix B.4. In fact, we prove that as $n\to\infty$, the weighted
conformal p-values converge to i.i.d. random variables that are dominated by
the uniform distribution, which ensures asymptotic FDR control (see Appendix
B.4 for details). The technical condition we impose resembles that in [STS04],
which ensures the existence of a limit point of the rejection threshold and
enables the asymptotic analysis.

We also remark that the PRDS property is a sufficient, but not necessary,
condition for the BH procedure to control the FDR. In fact, we will see that
the BH procedure with weighted conformal p-values empirically controls the FDR
in all of our numerical experiments, hence it remains a reasonable option in
practice.

##  3 Finite-sample FDR control

We now introduce a new multiple testing procedure, WCS, that controls the FDR
in finite samples with weighted conformal p-values. Our method builds on the
following p-values:

$\displaystyle{p}_{j}=\frac{\sum_{i=1}^{n}w(X_{i})\operatorname{\mathds{1}}{\\{V_{i}<\widehat{V}_{n+j}\\}}+w(X_{n+j})}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})}.$
(5)

It is a slight modification of (4) up to random tie-breaking. The asymptotic
results in Theorem 2.5 also hold with these p-values as long as the
distributions of the $V_{i}$’s and $\widehat{V}_{n+j}$’s do not have point
masses.

###  3.1 Weighted Conformalized Selection

As before, we compute $V_{i}=V(X_{i},Y_{i})$ for $i=1,\dots,n$ and
$\widehat{V}_{n+j}=V(X_{n+j},c_{n+j})$ for $j=1,\dots,m$. For each $j$, we
compute a set of auxiliary p-values

$\displaystyle{p}_{\ell}^{(j)}:=\frac{\sum_{i=1}^{n}w(X_{i})\operatorname{\mathds{1}}{\\{V_{i}<\widehat{V}_{n+\ell}\\}}+w(X_{n+j})\operatorname{\mathds{1}}{\\{\widehat{V}_{n+j}<\widehat{V}_{n+\ell}\\}}}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})},\quad\ell\neq
j.$ (6)

Then, we let $\widehat{\mathcal{R}}_{j\to 0}$ be the rejection set of BH
applied to
$\\{p_{1}^{(j)},\dots,p_{j-1}^{(j)},0,p_{j+1}^{(j)},\dots,p_{m}^{(j)}\\}$ at
the nominal level $q$. Note that $j\in\widehat{\mathcal{R}}_{j\to 0}$ by
default, and set $s_{j}=\frac{q|\widehat{\mathcal{R}}_{j\to 0}|}{m}$. We then
compute the first-step rejection set $\mathcal{R}^{(1)}:=\\{j\colon p_{j}\leq
s_{j}\\}$. Finally, we prune $\mathcal{R}^{(1)}$ to obtain the final selection
set $\mathcal{R}$ using any of the following three methods.

  1. (a)

_Heterogeneous pruning_ : generate i.i.d. random variables
$\xi_{j}\sim\textrm{Unif}[0,1]$, and set

$\displaystyle\mathcal{R}:=\mathcal{R}_{\textnormal{hete}}=\Big{\\{}j\colon
p_{j}\leq s_{j},~{}\xi_{j}|\widehat{\mathcal{R}}_{j\to 0}|\leq
r_{\textnormal{hete}}^{*}\Big{\\}},$ (7)

where
$r_{\textnormal{hete}}^{*}:=\max\big{\\{}r\colon\sum_{j=1}^{m}\operatorname{\mathds{1}}{\\{p_{j}\leq
s_{j},\,\xi_{j}|\widehat{\mathcal{R}}_{j\to 0}|\leq r\\}}\geq r\big{\\}}.$

  2. (b)

_Homogeneous pruning_ : generate an independent $\xi\sim\textrm{Unif}[0,1]$,
and set

$\displaystyle\mathcal{R}:=\mathcal{R}_{\textnormal{homo}}=\Big{\\{}j\colon
p_{j}\leq s_{j},~{}\xi|\widehat{\mathcal{R}}_{j\to 0}|\leq
r_{\textnormal{homo}}^{*}\Big{\\}},$ (8)

where
$r_{\textnormal{homo}}^{*}:=\max\big{\\{}r\colon\sum_{j=1}^{m}\operatorname{\mathds{1}}{\\{p_{j}\leq
s_{j},\,\xi|\widehat{\mathcal{R}}_{j\to 0}|\leq r\\}}\geq r\big{\\}}.$

  3. (c)

_Deterministic pruning_ : define the rejection set

$\displaystyle\mathcal{R}:=\mathcal{R}_{\textnormal{dtm}}=\Big{\\{}j\colon
p_{j}\leq s_{j},~{}|\widehat{\mathcal{R}}_{j\to 0}|\leq
r_{\textnormal{dtm}}^{*}\Big{\\}},$ (9)

where
$r_{\textnormal{dtm}}^{*}:=\max\big{\\{}r\colon\sum_{j=1}^{m}\operatorname{\mathds{1}}{\\{p_{j}\leq
s_{j},\,|\widehat{\mathcal{R}}_{j\to 0}|\leq r\\}}\geq r\big{\\}}.$ In all
options, we use $p_{j}$ defined in (5).

It is straightforward to see that
$\mathcal{R}_{\textnormal{dtm}}\subseteq\mathcal{R}_{\textnormal{hete}}$ and
$\mathcal{R}_{\textnormal{dtm}}\subseteq\mathcal{R}_{\textnormal{homo}}$,
i.e., random pruning leads to larger rejection sets. The selection procedure
is summarized in Algorithm 1.

Algorithm 1 Weighted Conformalized Selection

0: Calibration data $\\{(X_{i},Y_{i})\\}_{i=1}^{n}$, test data
$\\{X_{n+j}\\}_{j=1}^{m}$, thresholds $\\{c_{n+j}\\}_{j=1}^{m}$, weight
$w(\cdot)$, FDR target $q\in(0,1)$, monotone nonconformity score
$V\colon\mathcal{X}\times\mathcal{Y}\to\mathbb{R}$, pruning method
$\in\\{\texttt{hete},\texttt{homo},\texttt{dtm}\\}$.

1: Compute $V_{i}=V(X_{i},Y_{i})$ for $i=1,\dots,n$, and
$\widehat{V}_{n+j}=V(X_{n+j},c_{n+j})$ for $j=1,\dots,m$.

2: Construct weighted conformal p-values $\\{p_{j}\\}_{j=1}^{m}$ as in (5). \-
First-step selection -

3: for $j=1,\dots,m$ do

4: Compute p-values $\\{{p}_{\ell}^{(j)}\\}$ as in (6).

5: (BH procedure) Compute $k^{*}_{j}=\max\big{\\{}k\colon 1+\sum_{\ell\neq
j}\operatorname{\mathds{1}}\\{{p}_{\ell}^{(j)}\leq qk/m\\}\geq k\big{\\}}$.

6: Compute $\widehat{\mathcal{R}}_{j\to 0}=\\{j\\}\cup\\{\ell\neq
j\colon{p}_{\ell}^{(j)}\leq qk^{*}_{j}/m\\}$.

7: end for

8: Compute the first-step selection set
$\mathcal{R}^{(1)}=\\{j\colon{p}_{j}\leq q|\widehat{\mathcal{R}}_{j\to
0}|/m\\}$. \- Second-step pruning -

9: Compute $\mathcal{R}=\mathcal{R}_{\textrm{hete}}$ as in (7) or
$\mathcal{R}=\mathcal{R}_{\textrm{homo}}$ as in (8) or
$\mathcal{R}=\mathcal{R}_{\textrm{dtm}}$ as in (9).

9: Selection set $\mathcal{R}$.

###  3.2 Theoretical guarantee

The following theorem, whose proof is in Appendix C.1, shows exact finite-
sample FDR control.

######  Theorem 3.1.

Write $Z_{i}=(X_{i},Y_{i})$ for $i=1,\dots,m+n$, and
$\widetilde{Z}_{n+j}=(X_{n+j},c_{n+j})$ for $j=1,\dots,m$. Suppose
$\\{Z_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}P$
and
$\\{Z_{n+j}\\}_{j=1}^{m}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}Q$,
and (1) holds for the input weight function $w(\cdot)$ of Algorithm 1. Assume
that for each $j=1,\dots,m$, the samples in
$\\{Z_{1},\dots,Z_{n},Z_{n+j}\\}\cup\\{\widetilde{Z}_{n+\ell}\\}_{\ell\neq j}$
are mutually independent. Then with either
$\mathcal{R}\in\\{\mathcal{R}_{\textnormal{hete}},\mathcal{R}_{\textnormal{homo}},\mathcal{R}_{\textnormal{dtm}}\\}$,
it holds that

$\displaystyle\mathbb{E}\Bigg{[}\frac{\sum_{j=1}^{m}\operatorname{\mathds{1}}{\\{j\in\mathcal{R},Y_{n+j}\leq
c_{n+j}\\}}}{1\vee|\mathcal{R}|}\Bigg{]}\leq q$

(the expectation is taken over both calibration and test data).

As we allow for random thresholds $c_{n+j}$, the independence assumption rules
out the cases where the thresholds $c_{n+j}$ are adversarially chosen; the
reason is that they would break certain exchangeability properties between
scores $\\{\widehat{V}_{n+\ell}\\}_{\ell=1}^{m}$ and calibration scores
$\\{V_{i}\\}_{i=1}^{n}$. Independence holds if $c_{n+j}$ is pre-determined, or
more generally, if $c_{n+j}$ is a random variable associated with independent
test samples (such as when $\\{X_{n+j},c_{n+j},Y_{n+j}\\}_{j=1}^{m}$ are
i.i.d. tuples that are independent of the calibration data). Examples include
individual treatment effects studied in Section 5, and the drug-target
interaction prediction task studied in Section 4.2. Other examples with random
thresholds related to healthcare are in [JC22, Section 2.4].

While the $p_{j}$’s are still marginally stochastically larger than uniforms,
their complicated dependence requires additional care, and this is why we use
the auxiliary p-values $\\{p_{\ell}^{(j)}\\}$ as ‘calibrators’ to determine a
potentially different rejection rule than naively applying the BH procedure.
The proof of Theorem 3.1 relies on comparing our p-values (5) to oracle
p-values

$\displaystyle\bar{p}_{j}=\frac{\sum_{i=1}^{n}w(X_{i})\operatorname{\mathds{1}}{\\{V_{i}<{V}_{n+j}\\}}+w(X_{n+j})}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})}.$
(10)

As we shall see, replacing our conformal p-values with their oracle
counterparts does not change the first-step rejection set
$\widehat{\mathcal{R}}_{j\to 0}$. The advantage of working with the oracle
p-values is a friendly dependence structure, which ultimately ensures FDR
control. The main ideas of the argument are:

  * •

Randomization reduction: for each $j$, $\widehat{\mathcal{R}}_{j\to 0}$ can be
expressed as $\widehat{\mathcal{R}}_{j\to 0}=f(p_{j},\bm{p}_{-j}^{(j)})$,
where $\bm{p}_{-j}^{(j)}:=\\{p_{\ell}^{(j)}\\}_{\ell\neq j}$. For all three
pruning options, it holds that

$\textnormal{FDR}\leq\sum_{j=1}^{m}\mathbb{E}\bigg{[}\frac{\operatorname{\mathds{1}}\\{p_{j}\leq
q|f(p_{j},\bm{p}_{-j}^{(j)})|/m\\}}{|f(p_{j},{\bm{p}}_{-j}^{(j)})|}\bigg{]}.$

  * •

Leave-one-out: $\widehat{\mathcal{R}}_{j\to 0}$ remains the same if we replace
all $\\{p_{\ell}\\}$ with $\\{\bar{p}_{\ell}\\}$ defined above, i.e.,
$f(p_{j},\bm{p}_{-j}^{(j)})=f(\bar{p}_{j},\bar{\bm{p}}_{-j}^{(j)})$. Since
$\widehat{\mathcal{R}}_{j\to 0}$ is obtained by BH, it follows that
$f(\bar{p}_{j},\bar{\bm{p}}_{-j}^{(j)})=f(0,\bar{\bm{p}}_{-j}^{(j)})$.

  * •

Conditional independence: under the covariate shift (1), the oracle p-values
possess a nice conditional independence structure, namely,
$\bar{p}_{j}\perp\\!\\!\\!\perp\bar{\bm{p}}_{-j}^{(j)}{\,|\,}\mathcal{Z}_{j}$,
where $\mathcal{Z}_{j}=[Z_{1},\dots,Z_{n},Z_{n+j}]$ is an unordered set of
$Z_{i}=(X_{i},Y_{i})$, $i=1,\dots,n,n+j$. This, together with the first two
steps and the fact that $\bar{p}_{j}$ is stochastically larger than a uniform
random variable conditional on $\mathcal{Z}_{j}$, gives $\textnormal{FDR}\leq
q$.

Along the way, we shall explore interesting connections between Algorithm 1
and existing ideas in the multiple testing literature.

######  Remark 3.2.

Algorithm 1 is related to the conditional calibration method [FL22], which
utilizes sufficient statistics to calibrate a rejection threshold for each
individual hypothesis to achieve finite-sample FDR control. Indeed, for each
$j$, $s_{j}:={q|\widehat{\mathcal{R}}_{j\to 0}|}/{m}$ can be viewed as the
‘calibrated threshold’ for p-values in their framework; the unordered set
$\mathcal{Z}_{j}$ (after a careful leave-one-out analysis in our proof) plays
a similar role as a ‘sufficient statistic’. Our heterogeneous pruning is
similar to their random pruning step, while the other two options generalize
their approach. In addition, the procedure and its theoretical analysis are
specific to our problem, and they are significantly different.

######  Remark 3.3.

Algorithm 1 is also connected to the eBH procedure [WR20], a generalization of
the BH procedure to e-values. In the conventional setting with a set of
(deterministic) hypotheses $\\{H_{j}\\}_{j=1}^{m}$, e-values are nonnegative
random variables $\\{e_{j}\\}_{j=1}^{m}$ such that $\mathbb{E}[e_{j}]\leq 1$
if $H_{j}$ is null. For a target level $q\in(0,1)$, the eBH procedure outputs
$\mathcal{R}_{\textnormal{eBH}}:=\\{j\colon e_{j}\geq m/(q\widehat{k})\\}$,
where
$\widehat{k}=\max\big{\\{}k\colon\sum_{j=1}^{m}\operatorname{\mathds{1}}\\{e_{j}\geq
m/(qk)\\}\geq k\big{\\}}$. One can check that $\mathcal{R}_{\textnormal{dtm}}$
is equivalent to $\mathcal{R}_{\textnormal{eBH}}$ applied to

$\displaystyle e_{j}:=\frac{\operatorname{\mathds{1}}\\{p_{j}\leq
q|\widehat{\mathcal{R}}_{j\to 0}|/m\\}}{q|\widehat{\mathcal{R}}_{j\to
0}|/m},\quad j=1,\dots,m.$ (11)

Similar to (3), the $e_{j}$’s above obey a generalized notion of “null”
e-values, in the sense that

$\displaystyle\mathbb{E}\big{[}e_{j}\operatorname{\mathds{1}}\\{j\in\mathcal{H}_{0}\\}\big{]}\leq
1,\quad\textrm{for all}~{}j=1,\dots,m,$

see Appendix C.1 for details. Furthermore, using the generalized e-values
(11), $\mathcal{R}_{\textnormal{hete}}$ and $\mathcal{R}_{\textnormal{homo}}$
are equivalent to running eBH using $\\{e_{j}/\xi_{j}\\}$ and
$\\{e_{j}/\xi\\}$, respectively. Note that $\\{e_{j}/\xi_{j}\\}$ and
$\\{e_{j}/\xi_{j}\\}$ are no longer e-values, yet our procedures still control
the FDR while achieving higher power.

Our next result shows that the extra randomness in the second step does not
incur too much additional variation for $\mathcal{R}_{\textnormal{hete}}$ and
$\mathcal{R}_{\textnormal{homo}}$. The proof of Proposition 3.4 is in Appendix
C.3.

######  Proposition 3.4.

Suppose $\|w\|_{\infty}\leq M$ for some constant $M>0$. Suppose the
distributions of $\\{V(X_{i},Y_{i})\\}_{i=1}^{n}$ and
$\\{V(X_{n+j},c_{n+j})\\}_{j=1}^{m}$ have no point mass. Let
$\mathcal{R}_{\textnormal{BH}}$ be the rejection set of BH with weighted
conformal $p$-values (4), and $\mathcal{R}$ be any of the three selections
from Algorithm 1. Let $\mathcal{R}^{(1)}=\\{j\colon p_{j}\leq
q|\widehat{\mathcal{R}}_{j\to 0}|/m\\}$ be the first-step selection set. Then
the following holds:

  1. (i)

If $m$ is fixed, then
$\lim_{n\to\infty}\mathbb{P}(\mathcal{R}_{\textnormal{BH}}=\mathcal{R}=\mathcal{R}^{(1)})=1$
for each
$\mathcal{R}\in\\{\mathcal{R}_{\textnormal{homo}},\mathcal{R}_{\textnormal{hete}},\mathcal{R}_{\textnormal{dtm}}\\}$.

  2. (ii)

If $m,n\to\infty$ and the regularity conditions in case (ii) of Theorem 2.5
hold, then
$\frac{|\mathcal{R}^{(1)}\Delta\mathcal{R}_{\textnormal{BH}}|}{|\mathcal{R}^{(1)}|}~{}{\stackrel{{\scriptstyle\textnormal{a.s.}}}{{\to}}}~{}0$,
$\frac{|\mathcal{R}^{(1)}\Delta\mathcal{R}_{\textnormal{BH}}|}{|\mathcal{R}_{\textnormal{BH}}|}~{}{\stackrel{{\scriptstyle\textnormal{a.s.}}}{{\to}}}~{}0$,
$\frac{|\mathcal{R}^{(1)}\Delta\mathcal{R}_{\textnormal{homo}}|}{|\mathcal{R}^{(1)}|}~{}{\stackrel{{\scriptstyle\textrm{P}}}{{\to}}}~{}0$,
and
$\frac{|\mathcal{R}_{\textnormal{BH}}\Delta\mathcal{R}_{\textnormal{homo}}|}{|\mathcal{R}_{\textnormal{BH}}|}~{}{\stackrel{{\scriptstyle\textrm{P}}}{{\to}}}~{}0$.

Proposition 3.4 shows that when the size of the calibration set is
sufficiently large, our procedure is close to the BH procedure applied to the
weighted conformal p-values, the latter being a deterministic function of the
data. In particular, our first-step rejection set is asymptotically equivalent
to BH under mild conditions, and the same applies to
$\mathcal{R}_{\textnormal{homo}}$. The asymptotic analysis of
$\mathcal{R}_{\textnormal{hete}}$ and $\mathcal{R}_{\textnormal{dtm}}$ is
challenging. In our numerical experiments, we find that
$\mathcal{R}_{\textnormal{hete}}$ is often very close to
$\mathcal{R}_{\textnormal{homo}}$, while $\mathcal{R}_{\textnormal{dtm}}$
usually makes too few rejections.

###  3.3 FDR bounds with estimated weights

In practice, the covariate shift $w(\cdot)$ may be unknown. When the conformal
p-values are computed with some fitted weights $\widehat{w}(X_{i})$, Theorem
3.5 establishes upper bounds on the FDR. The proof is in Appendix C.4.

######  Theorem 3.5.

Take $w(\cdot):=\widehat{w}(\cdot)$ as the input weight function in Algorithm
1 and assume $\widehat{w}(\cdot)$ is estimated in a training process with data
independent from
$\\{(X_{i},Y_{i})\\}_{i=1}^{n}\cup\\{(X_{n+j},c_{n+j})\\}_{j=1}^{m}$. Under
the conditions of Theorem 3.1, for each
$\mathcal{R}\in\\{\mathcal{R}_{\textnormal{homo}},\mathcal{R}_{\textnormal{hete}},\mathcal{R}_{\textnormal{dtm}}\\}$,
we have

$\displaystyle\textnormal{FDR}\leq
q\cdot\mathbb{E}\bigg{[}\frac{\widehat{\gamma}^{2}}{1+q(\widehat{\gamma}^{2}-1)/m}\bigg{]},$

where
$\widehat{\gamma}:=\sup_{x\in\mathcal{X}}\max\\{\widehat{w}(x)/w(x),\,w(x)/\widehat{w}(x)\\}$.

Above, $\widehat{\gamma}\geq 1$ measures the estimation error in
$\widehat{w}(\cdot)$ relative to the true weights. When both $w(\cdot)$ and
$\widehat{w}(\cdot)$ are bounded away from zero and infinity,
$\widehat{\gamma}-1$ is of the same order as $\sup_{x}|\widehat{w}(x)-w(x)|$,
and hence the FDR inflation in Theorem 3.1 converges to zero if $\widehat{w}$
is consistent.

##  4 Application to drug discovery datasets

As a direct application of WCS, we consider the goal of prioritizing drug
candidates. We focus on two tasks: (i) drug property prediction, i.e.,
selecting molecules that bind to a target protein, and (ii) drug-target
interaction prediction, i.e., selecting drug-target pairs with high affinity
scores. We use the DeepPurpose library [HFG+20] for data pre-processing and
model training.

###  4.1 Drug property prediction

Our first goal is to find molecules that may bind to a target protein for HIV.
Machine learning models are trained on a subset of screened molecules from a
drug library, and then used to predict the remaining ones.

We use the HIV screening dataset in the DeepPurpose library with a total size
of $n_{\textrm{tot}}=41127$. In this dataset, the covariate $X\in\mathcal{X}$
is a string that represents the chemical structure of a molecule (encoded by
Extended-Connectivity FingerPrints, ECFP), and the response $Y\in\\{0,1\\}$ is
binary, indicating whether the molecule binds to the target protein. Our goal
is to select as many new drugs with $Y=1$ as possible while controlling the
FDR below a specified level. This can be viewed as the goal (2) with
$c_{n+j}\equiv 0.5$.

Oftentimes, experimenters introduce a bias by selecting the first batch of
molecules to screen (the training data), which results in a covariate shift
between training (calibration) and test data. Here, we mimic an
experimentation procedure that builds upon a pre-trained prediction model
$\widehat{\mu}\colon\mathcal{X}\to\mathbb{R}$ for binding activity, so that
those with higher predicted values are more likely to be included in the
training (calibration) data. In our experiment, to reduce computational time,
we train a single model for both predicting test samples and for selecting the
calibration fold. We take a subset of size $0.4\times n_{\textrm{tot}}$ as the
training set $\mathcal{D}_{\textnormal{train}}$, on which we train a small
neural network $\widehat{\mu}(\cdot)$ with three layers trained over three
epochs. The remaining $0.6\times n_{\textrm{tot}}$ samples are randomly
selected as the calibration set $\mathcal{D}_{\textnormal{calib}}$ with
probability $p(x)=\min\\{0.8,\sigma(\widehat{\mu}(x)-\bar{\mu})\\}$; here,
$x\in\mathcal{X}$,
$\bar{\mu}=\frac{1}{|\mathcal{D}_{\textnormal{train}}|}\sum_{i\in\mathcal{D}_{\textnormal{train}}}\widehat{\mu}(X_{i})$
is the average prediction on the training fold, and
$\sigma(t)=e^{t}/(1+e^{t})$ is the sigmoid function. The covariate shift (1)
is thus of the form $w(x)\propto 1/p(x)$, which we assume is known.

We compare the BH procedure with weighted conformal p-values (4), as well as
our Algorithms 1 and 2. The last one is applicable because we here take
$c=0.5$ for binary classification. We consider two scores used in BH and
Algorithm 1:

  * •

res: $V(x,y)=y-\widehat{\mu}(x)$.

  * •

clip: the score
$V(x,y)=M\cdot\operatorname{\mathds{1}}\\{y>0\\}-\widehat{\mu}(x)$, with
$M>2\sup_{x}|\widehat{\mu}(x)|$.

The empirical FDR over $N=200$ independent runs for FDR target levels
$q\in\\{0.1,0.2,0.5\\}$ is summarized in Figure 2. All algorithms empirically
control the FDR below the nominal levels (up to random fluctuation), showing
the reliability of WCS in prioritizing drug discovery under covariate shifts.
The two scores yield similar power as in Figure 3.

Figure 2: Empirical FDR for drug property prediction. The label wBH is a
shorthand for BH, and wCC.* for Algorithm 1 (for clip and res) or Algorithm 2
(for sub) with three pruning options
$*\in\\{\texttt{hete},\texttt{homo},\texttt{dtm}\\}$. The red dashed lines are
the nominal FDR levels. Figure 3: Empirical power for drug property
prediction. Everything else is as in Figure 2.

###  4.2 Drug-target interaction prediction

We then consider drug-target interaction (DTI) prediction, where the goal is
to select drug-target pairs with a high binding score. This task is relevant
if a therapeutic company is interested in prioritizing resources for drug
candidates that may be effective for any target they are interested in. We use
the DAVIS dataset [DHH+11], which records real-valued binding affinities for
$n_{\textrm{tot}}=30060$ drug-target pairs. The drugs and the targets are
encoded into numeric features using ECFP and Conjoint triad feature (CTF)
[SZL+07, STW+09].

We essentially follow the same procedure as in Section 4.1 and only rehearse
the main components. First, we randomly sample a subset of size
$n_{\textnormal{train}}=0.2\times n_{\textrm{tot}}$ as
$\mathcal{D}_{\textnormal{train}}$, on which we train a regression model
$\widehat{\mu}(\cdot)\colon\mathcal{X}\to\mathbb{R}$ using a $3$-layer neural
network with $10$ training epochs. This relatively simple model is suitable
for numerical experiments on CPUs (one can surely use more complicated
prediction models in practice). As before, we use the same model
$\widehat{\mu}(\cdot)$ for both predicting test samples and selecting
calibration data into first-batch screening (in practice, they can of course
be different). We sample a drug-target pair for inclusion in the calibration
set with probability $p(x)=\sigma(2\widehat{\mu}(x)-\bar{\mu})$, where
$\bar{\mu}$ is the average training prediction as before. Finally, among those
not selected in the calibration data, we randomly sample a subset of size
$m=5000$ as test data.

We choose a more complicated threshold $c_{j}$ for the continuous response.
For a drug-target pair $X_{n+j}$, we set $c_{j}$ to be the
$q_{\textrm{pop}}$-th quantile of the binding scores of all drug-target pairs
in $\mathcal{D}_{\textnormal{train}}$ with the same target. We evaluate
$q_{\textrm{pop}}\in\\{0.7,0.8\\}$. Thus,
$c_{j}:=c(X_{n+j},\mathcal{D}_{\textnormal{train}},q_{\textrm{pop}})$ where
$c$ is a mapping that takes both $\mathcal{D}_{\textnormal{train}}$ and the
target information in $X_{n+j}$ as inputs. Lastly, we evaluate the BH
procedure and Algorithm 1 with scores:

  * •

res: $V(x,y)=y-\widehat{\mu}(x)$.

  * •

clip:
$V(x,y)=M\cdot\operatorname{\mathds{1}}\\{y>c(x,\mathcal{D}_{\textnormal{train}},q_{\textrm{pop}})\\}+c(x,\mathcal{D}_{\textnormal{train}},q_{\textrm{pop}})\operatorname{\mathds{1}}\\{y\leq
c(x,\mathcal{D}_{\textnormal{train}},q_{\textrm{pop}})\\}-\widehat{\mu}(x)$ in
which $M=100$.

We always use $\mathcal{D}_{\textnormal{calib}}$ as the calibration set, and
set the FDR target as $q\in\\{0.1,0.2,0.5\\}$.

Figure 4 shows false discovery proportions (FDPs) for $q_{\textrm{pop}}=0.8$
in $N=200$ independent. Similar results for $q_{\textrm{pop}}=0.7$ are
presented in Appendix A.5. We see that the FDR is controlled at the desired
level for all algorithms and nonconformity scores. This shows the validity of
our algorithms and the plausibility of the independence assumptions we make on
the drug-target pairs. In this task, we do not observe much difference between
deterministic pruning (WCS.dtm) and the other two pruning options.
Furthermore, we observe that the FDPs across replications tightly concentrate
especially for clip and $q\in\\{0.2,0.5\\}$, showing that our algorithms are
stable with respect to data splitting (i.e., the randomness in choosing the
initial screening sets and in the training process). Comparing the two
nonconformity scores, we see that clip exploits the error budget and obtains a
realized FDR, which is very close to the nominal level, while res has a much
lower FDR in all cases. Not surprisingly, Figure 5 shows that clip has much
higher power.

Figure 4: Empirical FDR for drug-target interaction prediction with
$q_{\textrm{pop}}=0.8$. The shorthand WBH stands for BH procedure, and WCS.*
for Algorithm 1 with three pruning options
$*\in\\{\texttt{hete},\texttt{homo},\texttt{dtm}\\}$. The red dashed lines are
the nominal FDR levels. Solid lines are empirical averages (here empirical
FDR). Figure 5: Empirical power for drug-target interaction prediction with
$q_{\textrm{pop}}=0.8$. Everything else is as in Figure 4.

##  5 Application to multiple individual treatment effects

We now apply our method to infer individual treatment effects (ITEs) [LC21,
JRC23] under the potential outcomes framework [IR15]. The ITE describes the
difference between an individual’s outcomes when receiving a treatment versus
not; its variation comes from both individual characteristics and intrinsic
uncertainty in the outcomes. Inference on multiple ITEs is useful in assisting
reliable personalized decision making.

We consider a super-population setting, where
$(X_{i},O_{i}(1),O_{i}(0),T_{i})$ are drawn i.i.d. from a joint distribution
$\mathcal{P}$ (distinct from the distributions $P$ and $Q$ we use for the
data). Here, $X_{i}\in\mathcal{X}$ is the features, $T_{i}\in\\{0,1\\}$
indicates whether unit $i$ receives treatment, and $O_{i}(1),O_{i}(0)$ are the
potential outcomes under treatment and not, respectively. Under the standard
SUTVA [IR15], we observe $(X_{i},O_{i},T_{i})$, where
$O_{i}=O_{i}(T_{i})=T_{i}O_{i}(1)+(1-T_{i})O_{i}(0)$. We will specify how the
treatment $T_{i}$ is allocated later on. We focus on units _in the study_ ,
i.e., those who have received a certain treatment option. The ITE of unit $i$
is the random variable $\Delta_{i}=O_{i}(1)-O_{i}(0)$. As only one potential
outcome is observed, a crucial part of inferring the ITE is to predict the
counterfactual, i.e., $O_{i}(1)$ for those units with $T_{i}=0$, and
$O_{i}(0)$ for $T_{i}=1$. In the following, we consider simultaneous inference
for the ITEs of a set of units in the control group (inference for treated
units is similar).

Formally, the test samples are $\\{X_{n+j}\\}_{j=1}^{m}$ are units in the
control group (the outcome $O_{n+j}$ is observed). Our goal is to screen for
positive ITEs with FDR control, i.e., finding a subset
$\mathcal{R}\subseteq\\{1,\dots,m\\}$, such that

$\displaystyle\textnormal{FDR}:=\mathbb{E}\Bigg{[}\frac{\sum_{j=1}^{m}\operatorname{\mathds{1}}\\{O_{n+j}(1)\leq
O_{n+j}(0),~{}j\in\mathcal{R}\\}}{1\vee|\mathcal{R}|}\Bigg{]}\leq q.$ (12)

To cast the counterfactual problem in our framework, set $Y_{n+j}:=O_{n+j}(1)$
to be the unobserved outcome, and the thresholds $c_{n+j}:=O_{n+j}(0)$ to be
the observed outcomes. To infer $Y_{n+j}$, the training data is
$\\{(X_{i},Y_{i})\\}_{i=1}^{n}=\\{(X_{i},O_{i}(1))\\}_{i=1}^{n}$ for which
$T_{i}=1$. Conditional on the treatment status, the training samples are
i.i.d. from $P:=\mathcal{P}_{X,O(1){\,|\,}T=1}$, while the test samples
$\\{(X_{n+j},Y_{n+j})\\}_{j=1}^{m}$ are i.i.d. from
$Q:=\mathcal{P}_{X,O(1){\,|\,}T=0}$. The relationship between $P$ and $Q$ will
depend on the treatment assignment mechanism.

###  5.1 Warm-up: completely randomized experiments

As a warmup, consider a completely randomized experiment, where the treatment
assignments $T_{i}$ are i.i.d. draws from $\textrm{Bern}(\pi)$ for some
$\pi\in(0,1)$ independently from everything else. In this case,
$P=Q=\mathcal{P}_{X,O(1)}$ and it suffices to use conformal p-values without
weights. The testing procedure for randomized experiments has already been
studied in [JC22]. Take any _monotone_ nonconformity score $V(\cdot,\cdot)$
obtained from an independent training process, and compute
$V_{i}=V(X_{i},O_{i})=V(X_{i},O_{i}(1))$ for $i=1,\dots,n$, and
$\widehat{V}_{n+j}=V(X_{n+j},O_{n+j})=V(X_{n+j},O_{n+j}(0))$ for
$j=1,\dots,m$. Then construct conformal p-values (5) with $w(\cdot)\equiv 1$,
and run the BH procedure with these p-values at level $q$. FDR control is a
corollary of [JC22, Theorem 2.3], and we omit the proof.

######  Corollary 5.1.

The procedure above has FDR (12) at most $q$.

###  5.2 Stratified randomization and observational studies

We now consider a more general setting where the treatment assignment may
depend on the observed covariates, formalized as the following strong
ignorability condition [IR15].

######  Assumption 5.2 (Strong ignorability).

Under the joint distribution $\mathcal{P}$, it holds that
$(O(1),O(0))\perp\\!\\!\\!\perp T{\,|\,}X$. Equivalently, the treatment
assignments are independently generated from
$T_{i}\sim\textrm{Bern}(e(X_{i}))$, where $e\colon\mathcal{X}\to(0,1)$ is
known as the propensity score.

Assumption 5.2 is automatically satisfied in stratified randomization
experiments, where the treatment is randomized in a way that only depends on
the covariates, and the propensity score $e(\cdot)$ is known. In observational
studies where the treatment assignment mechanism is completely unknown,
Assumption 5.2 is standard for the identifiability of average treatment
effects [Ros02]. Inference on ITEs when this assumption is violated (i.e.,
when there is unmeasured confounding) is studied in [JRC23]; that said,
multiple testing under confounding needs additional techniques, and is beyond
the scope of the current work.

Under Assumption 5.2, the covariate shift condition (1) holds [LC21], as

$\displaystyle\frac{{\,\mathrm{d}}Q}{{\,\mathrm{d}}P}(x,y)=\frac{{\,\mathrm{d}}\mathcal{P}_{X,O(1){\,|\,}T=0}}{{\,\mathrm{d}}\mathcal{P}_{X,O(0){\,|\,}T=1}}(x,y)=\frac{{\,\mathrm{d}}\mathcal{P}_{X{\,|\,}T=0}}{{\,\mathrm{d}}\mathcal{P}_{X{\,|\,}T=1}}(x)=w(x):=\frac{\pi(1-e(x))}{(1-\pi)e(x)},$

where $e(x)=\mathcal{P}(T=1{\,|\,}X=x)$ is the propensity score, and
$\pi=\mathcal{P}(T=1)$ is the marginal probability of being treated. As such,
Algorithm 1 is readily applicable when $e(\cdot)$ is known.

######  Corollary 5.3.

Suppose $e(\cdot)$ is known, and Assumption 5.2 holds. Consider calibration
data $\\{(X_{i},O_{i}(1)\colon T_{i}=1\\}_{i=1}^{n}$, test data
$\\{X_{n+j}\colon T_{n+j}=1\\}_{j=1}^{m}$, thresholds $\\{O_{n+j}(0)\colon
T_{n+j}=0\\}_{j=1}^{m}$, any monotone score $V$, and weight function
$w(x)\propto(1-e(x))/{e(x)}$ as the input of Algorithm 1. Then any selection
procedure
$\mathcal{R}\in\\{\mathcal{R}_{\textnormal{hete}},\mathcal{R}_{\textnormal{homo}},\mathcal{R}_{\textnormal{dtm}}\\}$
has FDR at most $q$.

The propensity score function $e(\cdot)$ is unknown for observational data.
Under Assumption 5.2, we can estimate the propensity scores (hence the weight
function) using an independent training fold and plug this estimate into the
construction of p-values. As a corollary of Theorem 3.5, we can develop an FDR
bound for observational data.

######  Corollary 5.4.

In the setting of Corollary 5.3, take
$w(x)\propto(1-\widehat{e}(x))/{\widehat{e}(x)}$ as the input of Algorithm 1,
where $\widehat{e}(\cdot)$ is an estimate of $e(\cdot)$ which is independent
of the calibration and training data. Then any selection procedure
$\mathcal{R}\in\\{\mathcal{R}_{\textnormal{hete}},\mathcal{R}_{\textnormal{homo}},\mathcal{R}_{\textnormal{dtm}}\\}$
obeys

$\displaystyle\textnormal{FDR}\leq
q\cdot\mathbb{E}\bigg{[}\frac{\widehat{\gamma}^{2}}{1+q(\widehat{\gamma}^{2}-1)/m}\bigg{]},\qquad\widehat{\gamma}:=\sup_{x\in\mathcal{X}}\max\bigg{\\{}\frac{1-e(x)\widehat{e}(x)}{e(x)(1-\widehat{e}(x))},\frac{e(x)(1-\widehat{e}(x))}{1-e(x)\widehat{e}(x)}\bigg{\\}}.$

At a high level, the conformal p-values compare the observed outcome
$O_{n+j}(0)$ of the control units to the empirical distribution of $O_{i}(1)$
in the calibration data; the latter—after proper weighting—shows the typical
behavior of the counterfactual $O_{n+j}(1)$, and provides evidence for whether
it may be larger than $O_{n+j}(0)$, i.e., whether $\Delta_{n+j}$ is positive.

In predicting $O_{n+j}(1)$, we use the marginal distribution of $(X,O(1))$
from the calibration data, but ignore the information in the observed outcome
$O(0)$. Put differently, our inference for ITEs is valid regardless of how the
potential outcomes are coupled. Thus, the actual false discovery rate may be
lower than $q$. Indeed, it is observed in [JRC23] that the FDR for identifying
positive ITEs—even without adjusting for multiplicity—can sometimes be lower
than the nominal level. Next, we are to empirically investigate the FDR and
power of our method under various couplings of potential outcomes. A
theoretical understanding is left for future work.

###  5.3 Simulation studies

We design joint distributions of the variables $(X,O(1),O(0),T)$ with various
covariate distributions, regression functions, and couplings of potential
outcomes. The covariates $X_{i}\in\mathbb{R}^{10}$ ($p=10$) are obtained via
$X_{ij}=\Phi(X_{ij}^{0})$, $j=1,\dots,10$, where $X_{i}^{0}\in\mathbb{R}^{10}$
are i.i.d. $\mathcal{N}(0,\Sigma)$, and $\Phi(\cdot)$ is the CDF of a
$\mathcal{N}(0,1)$ random variable. We set $\Sigma=I_{p}$ for the independent
case and $\Sigma_{k,j}=0.9^{|k-j|}$ for the correlated case. 222Simulations
and real data analysis for ITE inference can be found at
https://github.com/ying531/conformal-selection. We consider three cases:

  * •

Setting 1: $O_{i}(0)=0.1\,\epsilon_{0,i}$,
$O_{i}(1)=\max\\{0,\mu_{1}(X_{i})+\sigma_{1}(X_{i})\,\epsilon_{1,i}\\}$, where
$\sigma_{1}(x)=0.2-\log x_{1}$ and
$\mu_{1}(x)={4}/{(1+e^{-12x_{1}-0.5})(1+e^{-12x_{2}-0.5})}$.

  * •

Setting 2 is a mixture of deterministic and stochastic ITEs. With probability
$0.1$, set $O_{i}(1)=0.1\,\epsilon_{0,i}-0.5$, and $O_{i}(0)=O_{i}(1)+0.05$;
otherwise, set them as in Setting 1.

  * •

Setting 3: $O_{i}(0)=\mu_{0}(X_{i})+0.1\,\epsilon_{0,i}$,
$O_{i}(1)=\max\\{0,\mu_{1}(X_{i})+\sigma_{1}(X_{i})\,\epsilon_{1,0}\\}$, where
$\sigma_{1}(\cdot)$ is as in Setting 1, and
$\mu_{0}(x)={2}/{(1+e^{-3x_{1}-0.5})(1+e^{-3x_{2}-0.5})}$,
$\mu_{1}(x)=0.1+1.5\,\mu_{0}(x)$.

Above, the variables ($\epsilon_{0,i},\epsilon_{1,i}$) are i.i.d. with
$\mathcal{N}(0,1)$ as marginals. We consider three coupling scenarios: (i)
independent, in which $\epsilon_{0,i}$ and $\epsilon_{1,i}$ are independent;
(ii) negative, in which $\epsilon_{1,i}=-\epsilon_{0,i}$; (iii) positive, in
which $\epsilon_{1,i}=\epsilon_{0,i}$. We generate treatment indicators via
$T_{i}\sim\textrm{Bern}(e(X_{i}))$ independently, where
$e(x)=(1+\textrm{Beta}_{2,4}(x_{1}))/{4}$, and $\textrm{Beta}_{2,4}(\cdot)$ is
the CDF of the Beta distribution with shape parameters $(2,4)$. Finally, we
set those $T_{i}=1$ as the calibration data ($n=250$), and those $T_{i}=0$ as
the test sample ($m=100$); as training data we select
$n_{\textnormal{train}}=750$ units, which are used to construct four
conformity scores:

  * •

reg: $V(x,y)=y-\widehat{\mu}(x)$, where $\widehat{\mu}(x)$ is an estimate of
$\mu_{1}(x)$ using regression forest from the grf R-package [ATW19].

  * •

cdf: $V(x,y)=\widehat{F}(x,y)$ [CWZ21], where $\widehat{F}(x,y)$ is an
estimate of $\mathcal{P}(O(1)\leq y{\,|\,}X=x)$. Here, we fit $\widehat{F}$
via inverting quantile regression forests from grf.

  * •

oracle: $V(x,y)=\mathcal{P}(O(1)\leq y{\,|\,}X=x)$, which is not computable
from data and is shown only for illustration.

  * •

cqr: $V(x,y)=y-\widehat{q}_{\beta}(x)$ [RPC19], where $\widehat{q}_{\beta}(x)$
is an estimate of the $\beta$-th quantile of $\mathcal{P}_{O(1){\,|\,}X=x}$,
$\beta\in\\{0.2,0.5,0.8\\}$, from the quantile regression forests in grf.

After generating the three folds of data, we fit $V(\cdot,\cdot)$ on the
training fold, and run four procedures on the calibration and training folds:
BH with weighted conformal p-values in (4) (named WBH in the plot), and
Algorithm 1 with three pruning options (named WCS.hete, WCS.homo, WCS.dtm,
respectively). The experiment is repeated for $N=1000$ independent runs.

Among the three, setting 2 is the most challenging: $10\%$ of the samples are
with slightly negative ITEs; they are likely to be falsely rejected since the
selection procedure ignores the coupling (it may reject a sample with a small
value of $O(0)$ even if $O(1)$ is smaller). We expect setting 1 to be the
least challenging as there is a strong signal in $O(1)$. Setting 3 is more
challenging than setting 1 since $O(1)$ and $O(0)$ are close to each other as
the regression functions are similar.

Figure 6: Empirical FDR for finding positive ITEs. The columns represent three
settings with either independent (ind) or correlated (corr) covariates. Each
row corresponds to a coupling of potential outcomes (positive, negative, and
no coupling). In each subplot, the $y$-axis is the nonconformity score, and
the $x$-axis is the empirical FDR. The red dashed lines indicate the nominal
level $q=0.10$.

##### FDR control.

The empirical FDR in all settings is plotted in Figure 6. We observe FDR
control for all methods, including the BH procedure. (We omit the results of
$\mathcal{R}_{\textnormal{dtm}}$ because it does not succeed in making any
selection.) Across settings with the same regression function, those with
correlated covariates often have higher FDR. The realized FDR also varies with
the scores: oracle incurs very low FDR, while its empirical counterpart cdf
has higher FDR. On the other hand, quantile regression based scores are robust
to the quantile level $\beta$ in the sense that cqr achieves similar FDR for
all choices of $\beta$.

The impact of coupling on FDR is less consistent across settings. The positive
coupling from setting 1 yields a low FDR, perhaps because it is difficult to
obtain sufficiently strong evidence which leads to fewer rejections (see the
power analysis below). The three couplings yield similar FDRs for setting 2.
For setting 3, positive coupling leads to the highest FDR.

##### Power.

The empirical power is defined as
$\textrm{Power}:=\mathbb{E}\Big{[}\frac{\sum_{j=1}^{m}\operatorname{\mathds{1}}\\{j\in\mathcal{R},O_{n+j}(1)>O_{n+j}(0)\\}}{\sum_{j=1}^{m}\operatorname{\mathds{1}}\\{O_{n+j}(1)>O_{n+j}(0)\\}}\Big{]}.$
As we see in Figure 7, the power is in general higher when the covariates are
correlated. The ITE depends on $(X_{1},X_{2})$, and in our design, higher
correlation between the first two entries leads to higher power. The power
also varies with the nonconformity scores. The regression-based reg leads to
lower power than other methods. This means that regression functions may be
underpowered in capturing the essential information for contrasting the
distributions of two potential outcomes. Quantile-regression based methods
(cqr) achieve similar power for different choices of $\beta$ in all settings.
Finally, fitting the conditional distribution function (cdf) is the most
powerful option in settings 1 and 3, but is far less powerful in setting 2.
Its oracle version (oracle) is also surprisingly less powerful; we observe
that the oracle cdf cannot create sufficiently small p-values, while the
estimated one is able to do so due to fluctuations caused by estimation
uncertainty.

Figure 7: Empirical power for finding positive ITEs. Everything else is as in
Figure 6.

The impact of coupling on the power is consistent across settings. In general,
negative coupling leads to the highest power, while positive coupling leads to
the lowest. This is because the conformal p-values compare the observed
outcome to the _marginal_ distribution of the counterfactuals, without
accounting for their joint distribution. Thus, when $O_{n+j}(0)$ is extremely
small, i.e., when we see a small p-value and $j\in\mathcal{R}$, under negative
coupling, it is more probable that $O_{n+j}(1)$ is relatively large, hence
leading to a true discovery. On the contrary, under positive coupling, a large
$O_{n+j}(1)$ often corresponds to a relatively large $O_{n+j}(0)$ and hence a
large conformal p-value, which may not be selected.

Finally, $\mathcal{R}_{\textnormal{dtm}}$ makes no selection. This may be due
to the threshold effect: since the selection requires $p_{j}\leq
q|\widehat{\mathcal{R}}_{j\to 0}|/m$ for test samples $j$ with the smallest
$|\widehat{\mathcal{R}}_{j\to 0}|$, the pruning step may exclude too many
candidates. We recommend using $\mathcal{R}_{\textnormal{hete}}$ and
$\mathcal{R}_{\textnormal{homo}}$ in practice.

##### Stability.

As aforementioned in Section 3, our method introduces extra randomness. It is
shown in Proposition 3.4 that this becomes asymptotically negligible, so that
our selection set is asymptotically the same as that the BH procedure yields.
We empirically evaluate the discrepancy in the selections, namely,
$|\mathcal{R}\Delta\mathcal{R}_{\textnormal{BH}}|/|\mathcal{R}_{\textnormal{BH}}|$,
for $\mathcal{R}$ being either $\mathcal{R}_{\textnormal{homo}}$ or
$\mathcal{R}_{\textnormal{hete}}$. The results are shown in Figure 8,
confirming our theory. Though we were not able to prove the asymptotic
equivalence of $\mathcal{R}_{\textnormal{hete}}$ and
$\mathcal{R}_{\textnormal{BH}}$ in the $m,n\to\infty$ regime, the discrepancy
is still small (despite being larger than that of
$\mathcal{R}_{\textnormal{homo}}$). All in all, the impact of additional
randomness seems acceptable in practice.

Figure 8: Discrepancy between selection rules
$|\mathcal{R}\Delta\mathcal{R}_{\textnormal{BH}}|/|\mathcal{R}_{\textnormal{BH}}|$.
Everything else is as in Figure 6.

The discrepancy for $\mathcal{R}_{\textnormal{dtm}}$ from
$\mathcal{R}_{\textnormal{BH}}$ is large in general is
$\mathcal{R}_{\textnormal{dtm}}$ is usually an empty set, and hence we did not
plot it. We conjecture that the stability of $\mathcal{R}_{\textnormal{dtm}}$
claimed in case (i) from Proposition 3.4 may not necessarily apply to the
moderately large sample sizes $m=100$ and $n=250$, and
$\mathcal{R}_{\textnormal{dtm}}$ may be far from
$\mathcal{R}_{\textnormal{BH}}$ for such sample sizes. In such cases, we
recommend using $\mathcal{R}_{\textnormal{homo}}$ or
$\mathcal{R}_{\textnormal{hete}}$ for better stability and power.

###  5.4 Real data analysis

We revisit the NSLM observational dataset from [CFM+19] and use our method to
detect multiple positive individual treatment effects. It is a semi-synthetic
observational dataset curated from a real randomized experiment, also analyzed
in [LC21, JRC23].

Figure 2 in [JRC23] plots the $\Gamma$-value, a measure of robustness of
positive ITEs against unmeasured confounding, versus individual covariates.
While units with larger $\Gamma$-values are of natural interest, properties of
inference hold on average, instead of over _selected_ units (e.g. those with
high $\Gamma$-values). We are to produce a similar plot (Figure 10) for
detecting positive ITEs, with guarantees over units that exhibit the strongest
evidence. (We do not consider the confounding issue, which may need nontrivial
extension of our techniques.)

We randomly subsample three disjoint folds of size $5000$, $1000$, and $391$
from the original dataset. The first is the training fold, and consists of
both treated and control units. The calibration data consists of all the
$n=997$ treated units in the second fold. The test data consists of all the
$m=256$ control units in the last fold. To deploy our procedure, we first
train a propensity score model $\widehat{e}(\cdot)$ using a regression forest
from the grf R-package, and set
$\widehat{w}(x)=(1-\widehat{e}(x))/{\widehat{e}(x)}$. We then use a quantile
regression forest from grf on the training data to fit $\widehat{q}_{0.5}(x)$,
the conditional median of $O(1)\mid X$, and set
$V(x,y)=y-\widehat{q}_{0.5}(x)$. Finally, we apply Algorithm 1 as in Theorem
3.1 to test for positive ITEs.

Figure 9 plots the empirical CDF of the weighted conformal p-values (5)
computed with $c_{n+j}:=O_{n+j}(0)$ on the test samples. These p-values are
stochastically smaller than Unif([0,1]), showing evidence for positive
treatment effects in general. Yet, to identify positive ITEs with rigorous
error control we must apply multiple testing ideas here.

Figure 9: Empirical CDF of the weighted conformal p-values for detecting
positive ITEs.

We observe
$\mathcal{R}_{\textnormal{BH}}=\mathcal{R}_{\textnormal{homo}}=\mathcal{R}_{\textnormal{hete}}$
while $\mathcal{R}_{\textnormal{dtm}}=\varnothing$ at FDR levels
$q\in\\{0.1,0.2,0.5\\}$. Figure 10 plots the weighted conformal p-values
versus school achievement levels of test units (each dot represents a test
unit) with red dots (from light to dark) identified as positive ITEs at
various FDR levels. Students in schools with moderate achievement levels
demonstrate the strongest evidence of benefiting from the treatment.

Figure 10: P-values and rejection sets for detecting positive ITEs among
control units.

##  6 Application to outlier detection

Finally, we study the extension of our framework for outlier detection, where
the calibration inliers may follow a distinct distribution compared with test
inliers. We discuss a new setup and apply a variant of Algorithm 1 to bank
marketing data.

###  6.1 Hypothesis-conditional FDR control

Assuming access to i.i.d. training data $\\{Z_{i}\\}_{i=1}^{n}$, we consider a
set of test data $\\{Z_{n+j}\\}_{j=1}^{m}$ for which the $Z_{n+j}$’s may only
be partially observed (e.g., if $Z=(X,Y)$ we would observe the features $X$
but not the response $Y$). We are interested in some null hypotheses
$\\{H_{j}\\}_{j=1}^{m}$ associated with the test samples. As before, whether
$H_{j}$ is true can be a random event depending on $Z_{n+j}$; this includes
our previous problem with $H_{j}\colon Y_{n+j}\leq c_{n+j}$.

####  6.1.1 Hypothesis-conditional covariate shift

######  Assumption 6.1 (Hypothesis-conditional covariate shift).

$\\{Z_{i}\\}_{i=1}^{n}\cup\\{Z_{n+j}\\}_{j=1}^{m}$ are mutually independent.
Also, conditional on the subset $\mathcal{H}_{0}\subset\\{1,\dots,m\\}$ of all
null hypotheses, it holds that
$\\{Z_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}P$,
and
$\\{Z_{n+j}\\}_{j\in\mathcal{H}_{0}}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}Q$,
where ${\,\mathrm{d}}Q/{\,\mathrm{d}}P(Z)=w(X)$ for some function
$w\colon\mathcal{X}\to\mathbb{R}^{+}$ and $X\subseteq Z$.

A special case of problem (2) obeys these assumptions.

######  Example 6.2 (Binary classification).

Set $Z=(X,Y)$. where $Y\in\\{0,1\\}$ is a binary response. Suppose the goal is
to find positive $Y_{n+j}$, e.g. an active drug or a qualified candidate. We
only observe the covariates $\\{X_{n+j}\\}_{j=1}^{m}$ for the test samples
$\\{(X_{n+j},Y_{n+j})\\}_{j=1}^{m}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}Q$.
Consider a reference dataset that only preserves positive samples among a set
of i.i.d. data from a covariate shifted distribution $P$ ; that is,
$\mathcal{D}_{\textnormal{calib}}=\\{Z_{i}\colon Y_{i}=0\\}$ where
$\\{(X_{i},Y_{i})\\}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}P$.
Although the (super-population) covariate shift (1) no longer applies,
conditional on $\mathcal{H}_{0}=\\{j\colon Y_{n+j}=0\\}$, it still holds that
$Z_{i}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}P_{Z{\,|\,}Y=0}$ for
$Z_{i}\in\mathcal{D}_{\textnormal{calib}}$, and
$Z_{n+j}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}Q_{Z{\,|\,}Y=0}$
for $j\in\mathcal{H}_{0}$.

The above example also applies to candidate screening with constant thresholds
$c_{n+j}\equiv c_{0}$. Setting
$\widetilde{Y}=\operatorname{\mathds{1}}\\{Y>c_{0}\\}\in\\{0,1\\}$ and
$\mathcal{D}_{\textnormal{calib}}=\\{Z_{i}\colon Y_{i}\leq c_{0}\\}$, all
arguments apply similarly to $\widetilde{Z}=(X,\widetilde{Y})$. However, this
does not necessarily apply when the thresholds $c_{n+j}$ are random variables
(especially when no ’threshold $c_{j}$’ is observed in the calibration data,
such as in the counterfactual inference problem we will study later).

Another application is outlier detection under covariate shifts.

######  Example 6.3 (Outlier detection).

We revisit outlier detection [BCL+21] while allowing for identifiable
covariate shifts between the calibration inliers and the test inliers. Given
$\\{Z_{i}\\}_{i=1}^{n}$ drawn i.i.d. from an unknown distribution $P$ and a
set of test data $\\{Z_{n+j}\\}_{j=1}^{m}$, we assume the inliers in the test
data are i.i.d. from a distribution $Q$ with
${\,\mathrm{d}}Q/{\,\mathrm{d}}P(Z)=w(Z)$ for a known function $w$, while
allowing outliers to be from arbitrary distributions. The covariate shift may
happen, for example, when inliers _were_ from $Q$ but the calibration set is
selected with preferences relying on $z$: for instance, one may include more
female users to balance the gender distribution when curating a reference
panel of normal transactions (inliers). In this case,
$\mathcal{H}_{0}=\\{j\colon Z_{n+j}\sim Q\\}$ is a deterministic set, and
Assumption 6.1 clearly holds.

The outlier detection example is closely related to identifying concept
drifts.

######  Example 6.4 (Concept shift detection).

Letting $Z=(X,Y)$ where $X\in\mathcal{X}$ is the family of covariates and $Y$
is the response, concept shift focuses on potential changes in the conditional
distribution of $Y$ given $X$. Given calibration data
$\\{Z_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}P$
and independent test data $\\{Z_{n+j}\\}_{j=1}^{m}$, [HL20] assume
$\\{Z_{n+j}\\}_{j=1}^{m}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}Q$,
and test for the global null
$H_{0}\colon{\,\mathrm{d}}Q/{\,\mathrm{d}}P(Z)=w(X)$ for some
$w\colon\mathcal{X}\to\mathbb{R}^{+}$. They achieve this by combining
independent p-values after sample splitting. Our framework can be used to test
individual concept drifts with dependent p-values. For instance, assume
$\\{X_{n+j}\\}_{j=1}^{m}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}Q_{X}$
for some unknown (but estimable) $Q_{X}$, we can test whether
$P_{Y_{n+j}{\,|\,}X_{n+j}}=P_{Y{\,|\,}X}$. The null hypotheses can be
formulated as $H_{j}\colon Z_{n+j}\sim Q$, where
${\,\mathrm{d}}Q/{\,\mathrm{d}}P(Z)=w(X)$ for some
$w\colon\mathcal{X}\to\mathbb{R}^{+}$ that is either known or can be estimated
well under proper conditions.

####  6.1.2 Multiple testing procedure

Our procedure for outlier detection under covariate shift is detailed in
Algorithm 2. This slightly modifies Algorithm 1 by removing the thresholds
(note differences in lines 1, 2, and 4). In the classification or constant
threshold problem (Example 6.2), it suffices to set $Z=X$ and leave out $Y$.

Algorithm 2 Hypothesis-conditional Weighted Conformalized Selection

0: Calibration data $\\{Z_{i}\\}_{i=1}^{n}$, test data
$\\{Z_{n+j}\\}_{j=1}^{m}$, weight function $w(\cdot)$, FDR target $q\in(0,1)$,
monotone nonconformity score
$V\colon\mathcal{X}\times\mathcal{Y}\to\mathbb{R}$, pruning method
$\in\\{\texttt{hete},\texttt{homo},\texttt{dtm}\\}$.

1: Compute $V_{i}=V(Z_{i})$ for $i=1,\dots,n$, and ${V}_{n+j}=V(Z_{n+j})$ for
$j=1,\dots,m$.

2: Construct weighted conformal p-values $\\{p_{j}\\}_{j=1}^{m}$ as in (5)
with $\widehat{V}_{n+j}$ replaced by $V_{n+j}$. \- First-step selection -

3: for $j=1,\dots,m$ do

4: Compute p-values $\\{{p}_{\ell}^{(j)}\\}$ as in (6) with
$\widehat{V}_{n+\ell}$ replaced by $V_{n+\ell}$ for all $\ell=1,\dots,m$.

5: (BH procedure) Compute $k^{*}_{j}=\max\big{\\{}k\colon 1+\sum_{\ell\neq
j}\operatorname{\mathds{1}}\\{{p}_{\ell}^{(j)}\leq qk/m\\}\geq k\big{\\}}$.

6: Compute $\widehat{\mathcal{R}}_{j\to 0}=\\{j\\}\cup\\{\ell\neq
j\colon{p}_{\ell}^{(j)}\leq qk^{*}_{j}/m\\}$.

7: end for

8: Compute the first-step selection set
$\mathcal{R}^{(1)}=\\{j\colon{p}_{j}\leq q|\widehat{\mathcal{R}}_{j\to
0}|/m\\}$. \- Second-step selection -

9: Compute $\mathcal{R}=\mathcal{R}_{\textrm{hete}}$ or
$\mathcal{R}=\mathcal{R}_{\textrm{homo}}$ or
$\mathcal{R}=\mathcal{R}_{\textrm{dtm}}$ as in Algorithm 1.

9: Selection set $\mathcal{R}$.

Algorithm 2 returns to the conventional perspective, where the null set is
deterministic, and the null p-values are dominated by Unif$([0,1])$. That is,
for $p_{j}$ constructed in Line 2 of Algorithm 2, it holds that

$\displaystyle\mathbb{P}(p_{j}\leq t{\,|\,}j\in\mathcal{H}_{0})\leq
t\quad\textrm{for all }t\in[0,1].$

After conditioning on $\mathcal{H}_{0}$, we no longer need to deal with the
randomness of the hypotheses and their interaction with the p-values. The only
issue is the mutual dependence among the p-values, which is addressed using a
similar idea as in our theoretical analysis of Algorithm 1.

Using calibration data obeying the covariate shift assumption, Algorithm 2
achieves a slightly stronger hypotheses-conditional FDR control. The proof of
Theorem 6.5 is in Appendix C.5.

######  Theorem 6.5.

Under Assumption 6.1, Algorithm 2 yields

$\displaystyle\mathbb{E}\bigg{[}\frac{|\mathcal{R}\cap\mathcal{H}_{0}|}{1\vee|\mathcal{R}|}\,\bigg{|}\,\mathcal{H}_{0}\bigg{]}\leq
q\cdot\frac{|\mathcal{H}_{0}|}{m}$

for any fixed $q\in(0,1)$, and each
$\mathcal{R}\in\\{\mathcal{R}_{\textnormal{homo}},\mathcal{R}_{\textnormal{hete}},\mathcal{R}_{\textnormal{dtm}}\\}$.

####  6.1.3 Comparison with Algorithm 1

In binary classification, or more generally, WCS with a constant threshold, we
have shown in Example 6.2 that Algorithm 2 yields FDR control. In this case,
Algorithms 1 and 2 differ in terms of (i) power, and (ii) distributional
assumptions. We elaborate on these distinctions.

First, Algorithm 2 only uses a subset of calibration data to construct
p-values, which leads to a power loss for specific choices of nonconformity
scores; see [JC22, Appendix A.1] for the i.i.d. case. Let us consider the
binary setting. Suppose we have access to a set of calibration data
$\\{(X_{i},Y_{i})\\}$ consisting of both $Y=1$ and $Y=0$ samples. As discussed
in Example 6.2, Assumption 2 holds if we only use data in the subset
$\mathcal{I}_{0}=\\{i\colon Y_{i}=0\\}$ as $\mathcal{D}_{\textnormal{calib}}$
in Algorithm 2. In contrast, Algorithm 1 uses all data points. Suppose we set
$V(x,y)=My-\widehat{\mu}(x)$ and $c_{n+j}\equiv 0$ in Algorithm 1, where
$\widehat{\mu}(\cdot)$ is a fitted point prediction, and
$M>2\sup_{x\in\mathcal{X}}|\widehat{\mu}(x)|$ is a sufficiently large
constant. Similarly, we set $V(x)=M-\widehat{\mu}(x)$ in Algorithm 2. This
construction ensures

$\displaystyle\inf_{x\in\mathcal{X}}V(x,1)=M-\sup_{x\in\mathcal{X}}\widehat{\mu}(x)>\sup_{x\in\mathcal{X}}|\widehat{\mu}(x)|\geq\sup_{x\in\mathcal{X}}V(x,0).$
(13)

Theorems 3.1 and 6.5 state that the FDR is controlled for both approaches.
However, letting $p_{j}$ denote the p-values constructed in Algorithm 1 and
$p_{j}^{\prime}$ denote those in Algorithm 2, we note that

$\displaystyle p_{j}$
$\displaystyle=\frac{\sum_{i\in\mathcal{I}_{0}}w(X_{i})\operatorname{\mathds{1}}{\\{V(X_{i},0)<V(X_{n+j},0)\\}}+w(X_{n+j})}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})}+\frac{\sum_{i\in\mathcal{I}_{1}}w(X_{i})\operatorname{\mathds{1}}{\\{V(X_{i},1)<V(X_{n+j},0)\\}}}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})}$
$\displaystyle=\frac{\sum_{i\in\mathcal{I}_{0}}w(X_{i})\operatorname{\mathds{1}}{\\{V(X_{i},0)<V(X_{n+j},0)\\}}+w(X_{n+j})}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})}$
$\displaystyle<\frac{\sum_{i\in\mathcal{I}_{0}}w(X_{i})\operatorname{\mathds{1}}{\\{V(X_{i},0)<V(X_{n+j},0)\\}}+w(X_{n+j})}{\sum_{i\in\mathcal{I}_{0}}w(X_{i})+w(X_{n+j})}=p_{j}^{\prime},$

where the second lines uses (13). That is, with this nonconformity score
(which is shown in [JC22] to be powerful), the p-values constructed in
Algorithm 1 is strictly smaller than those in Algorithm 2, leading to larger
rejection sets and higher power. Furthermore, in this case,

$\displaystyle\frac{p_{j}}{p_{j}^{\prime}}=\frac{\sum_{i\in\mathcal{I}_{0}}w(X_{i})+w(X_{n+j})}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})}$

is the weighted proportion of negative calibration samples. Thus, the power
gain of Algorithm 1 is more significant when there are more positive samples
in the test distribution.

Second, Algorithm 2 is suitable for dealing with imbalanced data, such as
those encountered in drug discovery. For instance, after selecting a subset of
molecules or compounds for virtual screening, the experimenter may discard a
few samples with $Y=0$ as she would deem them as uninteresting. This bias
would make Algorithm 1 inapplicable. However, if the decision to discard or
not does not depend on $X$, the shift between $P_{X{\,|\,}Y=0}$ and
$Q_{X{\,|\,}Y=0}$ remains the same as the covariate shift incurred by
selection into screening, and Algorithm 2 still provides reliable selection.

We finally illustrate the hypothesis-conditional variant from Algorithm 2 on
outlier detection tasks, where there is a covariate shift between calibration
inliers and test inliers. We focus on scenarios where the covariate shift is
known. For instance, if the demographics in normal financial transactions are
adjusted by rejection sampling to balance between male and female, urban and
rural users, and so on, then the covariate shift is given by the sampling
weights in the adjustment.333Results in this section can be reproduced at
https://github.com/ying531/conformal-selection.

###  6.2 Simulation studies

We adapt the simulation setting in [BCL+21] to a scenario with covariate shift
induced by rejection sampling. We fix a test sample size
$n_{\textnormal{test}}=1000$ and calibration sample size
$n_{\textnormal{calib}}=1000$.

At the beginning of the experiment, we sample a subset
$\mathcal{W}\subseteq\mathbb{R}^{50}$ with $|\mathcal{W}|=50$, where each
element in $\mathcal{W}$ is independently drawn from
$\textrm{Unif}([-3,3]^{50})$, and hold it as fixed afterwards. Fix a
proportion of outliers at $\rho_{\texttt{pop}}\in\\{0.1,0.2,\dots,0.5\\}$. The
number of outliers in the test data is
$n_{\textnormal{test}}\cdot\rho_{\texttt{pop}}$, where each of them is i.i.d.
generated as $X_{n+j}=\sqrt{a}V_{n+j}+W_{n+j}$ for signal strength $a$ varying
in the range $[1,4]$, $V_{n+j}\sim N(0,I_{50})$, and
$W_{n+j}\sim\textrm{Unif}(\mathcal{W})$. Following this, the test inliers are
i.i.d. generated as $X_{n+j}=V_{n+j}+W_{n+j}$, whose distribution is denoted
as $Q_{X}$. The calibration inliers are i.i.d. generated from $P_{X}$ with
${\,\mathrm{d}}Q_{X}/{\,\mathrm{d}}P_{X}(x)=w(x)\propto\sigma(x^{\top}\theta)$
($\sigma(\cdot)$ is the sigmoid), where $\theta\in\mathbb{R}^{50}$ and
$\theta_{j}=0.1\cdot\operatorname{\mathds{1}}\\{j\leq 5\\}$. We also generate
$n_{\textnormal{train}}=1000$ training sample from $P_{X}$. This setting
mimics a stylized scenario where the calibration inliers are collected in a
way such that the preference is prescribed by a logistic function of $X$, and
is known to the data analyst.

We train a one-class SVM $\widehat{\mu}(x)$ with rbf kernel using the scikit-
learn Python library, and apply Algorithm 2 with $Z_{i}=X_{i}$ at FDR level
$q=0.1$ using all the three pruning methods. All procedures are repeated
$N=1000$ times, and the FDPs and proportions of true discoveries are averaged
to estimate the FDR and power. We also evaluate the BH procedure applied to
p-values constructed in Line 4 of Algorithm 2. Figure 11 shows the empirical
FDR across runs. In line with Theorem 6.5, we see that the FDR is always below
$(1-\rho_{\texttt{pop}})q$. Also, BH applied to weighted conformal p-values
empirically controls the FDR and demonstrates comparable performance.

Figure 11: Empirical FDR averaged over $N=1000$ runs under increasing levels
$a$ of signal strength. Each subplot corresponds to one value of
$\rho_{\texttt{pop}}$. The red dashed lines are at the nominal level $q=0.1$.

Figure 12 plots the power of the four methods. Interestingly, although they
differ in FDR, especially when the signal strength is small, they achieve
nearly identical power, finding about the same number of true discoveries.

Figure 12: Empirical power averaged over $N=1000$ runs. Everything else is as
in Figure 11.

###  6.3 Real data application

We then apply Algorithm 2 to a bank marketing dataset [MCR14]; this dataset
has been used to benchmark outlier detection algorithms [PSvdH19, PSCH21]. In
short, each sample in the dataset represents a phone contact to a potential
client, whose demographic information such as age, gender, and marital status
is also recorded in the features. The outcome is a binary indicator $Y$, where
$Y=1$ represents a successful compaign. Positive samples are relatively rare,
accounting for about $10\%$ of all records.

In our experiments, we always use negative samples as the calibration data. We
find this dataset interesting because the classification and outlier detection
perspectives are somewhat blurred here. Viewing this as a classification task,
Example 6.2 becomes relevant, and our theory implies FDR control as long as
the distribution of negative samples (i.e., that of $X$ given $Y=0$) admits
the covariate shift. In the sense of finding positive responses (so as to
reach out to these promising clients), controlling the FDR ensures efficient
use of campaign resources. Alternatively, from an outlier detection
perspective, those $Y=1$ are outliers [PSvdH19] that deserve more attention.
Taking either of the two perspectives leads to the same calibration set and
the same guarantees following Section 6.1. Perhaps the only methodological
distinction is whether positive samples are leveraged in the training process.
Similar considerations also appeared in [LSS22], but the scenario therein is
more sophisticated than ours as they also use known outliers (positive
samples) in constructing p-values. As a classification problem, we may use
positive samples to train a classifier and produce nonconformity scores. As an
outlier detection problem, however, the widely-used one-class classification
relies exclusively on inliers (negative samples). We will evaluate the
performance of both (i.e., train the nonconformity score with or without
positive samples); see procedures cond_class and outlier below.

For comparison, we additionally evaluate Algorithm 1 which operates under a
covariate shift assumption on the joint distribution of $(X,Y)$ (see procedure
sup_class below); this is in contrast to Algorithm 2 which puts no assumption
on positive samples.

The total sample size is $N=41188$ with $4640$ positive samples. We use
rejection sampling to create the covariate shift between training/calibration
and test samples. In details, we use a subset of features
$X^{*}\in\mathbb{R}^{4}$ representing age and indicators of marriage, basic
four-year education, and housing. We first randomly select a subset of the
original dataset as the test data. Each sample enters the test data with
probability $e(x):=0.125\,\sigma({\theta^{\top}x^{*}})$, where
$\theta=(1,1/2,1/2,1/2)^{\top}$. This creates a covariate shift $w(x)\propto
e(x)/(1-e(x))$ between null samples in calibration and test folds, so that the
test set contains more senior, married, well-educated and housed people.

To test the two perspectives, we consider three procedures with FDR levels
$q\in\\{0.2,0.5,0.8\\}$:

  1. (i)

sup_class: We randomly split the data that is not in the test fold into two
equally-sized halves as the training and calibration folds,
$\mathcal{D}_{\textnormal{train}}$ and $\mathcal{D}_{\textnormal{calib}}$. We
use $\mathcal{D}_{\textnormal{train}}$ to train an SVM classifier using the
scikit-learn Python library with rbf kernel. Then, we apply Algorithm 1 with
$V(x,y)=100\cdot y-\widehat{\mu}(x)$, where $\widehat{\mu}(\cdot)$ is the
class probability output from the SVM classifier.

  2. (ii)

cond_class: With the same random split as in (i), the first half is used as
$\mathcal{D}_{\textnormal{train}}$, while we set
$\mathcal{D}_{\textnormal{calib}}$ as the set of negative samples in the
second half. We use $\mathcal{D}_{\textnormal{train}}$ to train an SVM
classifier in the same way as in (i) to obtain the predicted class probability
$\widehat{\mu}(\cdot)$. Lastly, we apply Algorithm 2 with
$V(x)=-\widehat{\mu}(x)$.

  3. (iii)

outlier: With the same split, we set $\mathcal{D}_{\textnormal{train}}$ and
$\mathcal{D}_{\textnormal{calib}}$ as the two folds after discarding all the
negative samples. We use $\mathcal{D}_{\textnormal{train}}$ to train a one-
class SVM classifier using the scikit-learn Python library with rbf kernel and
obtain the output $\widehat{\mu}(\cdot)$. We then apply Algorithm 2 with
$V(x)=-\widehat{\mu}(x)$.

To ensure consistency, in all procedures, the SVM-based classifier uses the
parameter gamma$=0.01$.444We find that this parameter leads to the highest
power for outlier detection, and is also a relatively powerful choice for
classification. FDPs and power (proportions of true discoveries) across
$N=1000$ independent runs are shown in Figures 13 and 14, respectively.

Figure 13: FDPs over $N=1000$ independent data splits. Each subplot
corresponds to one nominal FDR level shown by means of red dashed lines.

In Figure 13, we observe FDR control for all the methods in all settings.
Across the four testing methods, the FDPs are very similar, hence all of them
are reasonable choices. However, both the values and the variability of the
FDP vary across settings. The methods cond_class and outlier from Algorithm 2
do not use all of the error budget—note the $\pi_{0}$ factor in Theorem 6.5.
In contrast, Algorithm 1 leverages super-population structure (which is
present in this problem) and leads to tight FDPs and FDR around the target
level. In all cases, the outlier detection approach, where only negative
samples are used in the training process, has more variable FDPs. With this
dataset, the variability of FDP around the FDR is visible for $q=0.2$, but it
decreases as we increase the FDR level.

Figure 14: Proportion of true discoveries (empirical power). Everything else
is as in Figure 13.

In Figure 14, the power across the four multiple testing methods is similar.
However, we observe drastic differences among the three settings. While using
the same classifier, cond_class has lower power than sup_class, per our
discussion in Section 6.1.3. We also see that outlier has much lower power,
which may be due to the fact that the nonconformity score is not sufficiently
powerful in distinguishing inliers from outliers as the training process does
not use outliers. Our results suggest that in outlier detection problems, even
when we do not want to impose any distributional assumption on outliers (so
that Algorithm 1 becomes less reasonable), utilizing known outliers may be
helpful in obtaining better nonconformity scores.

Finally, we observe negligible difference between the rejection sets returned
by the BH procedure and Algorithm 1. There at most five distinct decisions
among $\approx 3700$ test samples.

## Acknowledgement

The authors thank John Cherian, Issac Gibbs, Kevin Guo, Kexin Huang, Jayoon
Jang, Lihua Lei, Shuangning Li, Zhimei Ren, Hui Xu, and Qian Zhao for helpful
discussions. E.C. and Y.J. were supported by the Office of Naval Research
grant N00014-20-1-2157, the National Science Foundation grant DMS-2032014, the
Simons Foundation under award 814641, and the ARO grant 2003514594.

## References

  * [AHBC17] Ernst Ahlberg, Oscar Hammar, Claus Bendtsen, and Lars Carlsson.  Current application of conformal prediction in drug discovery.  Annals of Mathematics and Artificial Intelligence, 81(1):145–154, 2017.
  * [ATW19] Susan Athey, Julie Tibshirani, and Stefan Wager.  Generalized random forests.  The Annals of Statistics, 47(2):1148–1178, 2019.
  * [BCL+21] Stephen Bates, Emmanuel Candès, Lihua Lei, Yaniv Romano, and Matteo Sesia.  Testing for outliers with conformal p-values.  arXiv preprint arXiv:2104.08279, 2021.
  * [BH95] Yoav Benjamini and Yosef Hochberg.  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  Journal of the Royal statistical society: series B (Methodological), 57(1):289–300, 1995.
  * [BHRZ23] Yajie Bao, Yuyang Huo, Haojie Ren, and Changliang Zou.  Selective conformal inference with fcr control.  arXiv preprint arXiv:2301.00584, 2023.
  * [BY01] Yoav Benjamini and Daniel Yekutieli.  The control of the false discovery rate in multiple testing under dependency.  Annals of statistics, pages 1165–1188, 2001.
  * [CCB19] Isidro Cortés-Ciriano and Andreas Bender.  Concepts and applications of conformal prediction in computational drug discovery.  arXiv preprint arXiv:1908.03569, 2019.
  * [CDLM21] Devin Caughey, Allan Dafoe, Xinran Li, and Luke Miratrix.  Randomization inference beyond the sharp null: Bounded null hypotheses and quantiles of individual treatment effects.  arXiv preprint arXiv:2101.09195, 2021.
  * [CFM+19] Carlos Carvalho, Avi Feller, Jared Murray, Spencer Woody, and David Yeager.  Assessing treatment effect variation in observational studies: Results from a data challenge.  Observational Studies, 5(2):21–35, 2019.
  * [CH18] Kristy A Carpenter and Xudong Huang.  Machine learning-based virtual screening and its applications to alzheimer’s drug discovery: a review.  Current pharmaceutical design, 24(28):3347–3358, 2018.
  * [CWZ21] Victor Chernozhukov, Kaspar Wüthrich, and Yinchu Zhu.  Distributional conformal prediction.  Proceedings of the National Academy of Sciences, 118(48), 2021.
  * [DDJ+21] Suresh Dara, Swetha Dhamercherla, Surender Singh Jadav, CH Babu, and Mohamed Jawed Ahsan.  Machine learning in drug discovery: A review.  Artificial Intelligence Review, pages 1–53, 2021.
  * [DHH+11] Mindy I Davis, Jeremy P Hunt, Sanna Herrgard, Pietro Ciceri, Lisa M Wodicka, Gabriel Pallares, Michael Hocker, Daniel K Treiber, and Patrick P Zarrinkar.  Comprehensive analysis of kinase inhibitor selectivity.  Nature biotechnology, 29(11):1046–1051, 2011.
  * [DWR21] Boyan Duan, Larry Wasserman, and Aaditya Ramdas.  Interactive identification of individuals with positive treatment effect while controlling false discoveries.  arXiv preprint arXiv:2102.10778, 2021.
  * [FL22] William Fithian and Lihua Lei.  Conditional calibration for false discovery rate control under dependence.  The Annals of Statistics, 50(6):3091–3118, 2022.
  * [FL23] Clara Fannjiang and Jennifer Listgarten.  Is novelty predictable?  arXiv preprint arXiv:2306.00872, 2023.
  * [FRTT12] Evanthia Faliagka, Kostas Ramantas, Athanasios Tsakalidis, and Giannis Tzimas.  Application of machine learning algorithms to an online recruitment system.  In Proc. International Conference on Internet and Web Applications and Services, pages 215–220. Citeseer, 2012.
  * [FSJB22] Adam Fisch, Tal Schuster, Tommi Jaakkola, and Regina Barzilay.  Conformal prediction sets with limited false positives.  In International Conference on Machine Learning, pages 6514–6532. PMLR, 2022.
  * [HFG+20] Kexin Huang, Tianfan Fu, Lucas M Glass, Marinka Zitnik, Cao Xiao, and Jimeng Sun.  Deeppurpose: A deep learning library for drug-target interaction prediction.  Bioinformatics, 2020.
  * [HL20] Xiaoyu Hu and Jing Lei.  A distribution-free test of covariate shift using conformal prediction.  arXiv preprint arXiv:2010.07147, 2020.
  * [Hua07] Ziwei Huang.  Drug discovery research: new frontiers in the post-genomic era.  John Wiley & Sons, 2007.
  * [IR15] Guido W. Imbens and Donald B. Rubin.  Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction.  Cambridge University Press, 2015.
  * [JC22] Ying Jin and Emmanuel J Candès.  Selection by prediction with conformal p-values.  arXiv preprint arXiv:2210.01408, 2022.
  * [JRC23] Ying Jin, Zhimei Ren, and Emmanuel J Candès.  Sensitivity analysis of individual treatment effects: A robust conformal inference approach.  Proceedings of the National Academy of Sciences, 120(6):e2214889120, 2023.
  * [KMLH17] Alexios Koutsoukas, Keith J Monaghan, Xiaoli Li, and Jun Huan.  Deep-learning: investigating deep neural networks hyper-parameters and comparison of performance to shallow methods for modeling bioactivity data.  Journal of cheminformatics, 9(1):1–13, 2017.
  * [Krs21] Damjan Krstajic.  Critical assessment of conformal prediction methods applied in binary classification settings.  Journal of Chemical Information and Modeling, 61:4823–4826, 10 2021\.
  * [LC21] Lihua Lei and Emmanuel J Candès.  Conformal inference of counterfactuals and individual treatment effects.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2021.
  * [LDG13] Antonio Lavecchia and Carmen Di Giovanni.  Virtual screening strategies in drug discovery: a critical review.  Current medicinal chemistry, 20(23):2839–2860, 2013.
  * [LSS22] Ziyi Liang, Matteo Sesia, and Wenguang Sun.  Integrative conformal p-values for powerful out-of-distribution testing with labeled outliers.  arXiv preprint arXiv:2208.11111, 2022.
  * [MCR14] Sérgio Moro, Paulo Cortez, and Paulo Rita.  A data-driven approach to predict the success of bank telemarketing.  Decision Support Systems, 62:22–31, 2014.
  * [MLMR22] Ariane Marandon, Lihua Lei, David Mary, and Etienne Roquain.  Machine learning meets false discovery rate.  arXiv preprint arXiv:2208.06685, 2022.
  * [NCBE14] Ulf Norinder, Lars Carlsson, Scott Boyer, and Martin Eklund.  Introducing conformal prediction in predictive modeling. a transparent and flexible alternative to applicability domain determination.  Journal of Chemical Information and Modeling, 54:1596–1603, 6 2014\.
  * [PPS+15] Sebastian Polak, Michael K Pugsley, Norman Stockbridge, Christine Garnett, and Barbara Wiśniowska.  Early drug discovery prediction of proarrhythmia potential and its covariates, 2015.
  * [PSCH21] Guansong Pang, Chunhua Shen, Longbing Cao, and Anton Van Den Hengel.  Deep learning for anomaly detection: A review.  ACM computing surveys (CSUR), 54(2):1–38, 2021.
  * [PSvdH19] Guansong Pang, Chunhua Shen, and Anton van den Hengel.  Deep anomaly detection with deviation networks.  In Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, pages 353–362, 2019.
  * [Ros02] P.R. Rosenbaum.  Observational Studies.  Springer Series in Statistics. Springer, 2002.
  * [RPC19] Yaniv Romano, Evan Patterson, and Emmanuel Candes.  Conformalized quantile regression.  Advances in neural information processing systems, 32, 2019.
  * [SAN+18] Fredrik Svensson, Natalia Aniceto, Ulf Norinder, Isidro Cortes-Ciriano, Ola Spjuth, Lars Carlsson, and Andreas Bender.  Conformal regression for quantitative structure–activity relationship modeling—quantifying prediction uncertainty.  Journal of Chemical Information and Modeling, 58(5):1132–1140, 2018\.
  * [SKM07] 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.
  * [SNB17] Fredrik Svensson, Ulf Norinder, and Andreas Bender.  Improving screening efficiency through iterative screening using docking and conformal prediction.  Journal of chemical information and modeling, 57(3):439–444, 2017\.
  * [SS16] Muhammad Ahmad Shehu and Faisal Saeed.  An adaptive personnel selection model for recruitment using domain-driven data mining.  Journal of Theoretical and Applied Information Technology, 91(1):117, 2016.
  * [SSKMM19] Sima Sajjadiani, Aaron J Sojourner, John D Kammeyer-Mueller, and Elton Mykerezi.  Using machine learning to translate applicant work history into predictors of performance and turnover.  Journal of Applied Psychology, 104(10):1207, 2019.
  * [STS04] John D Storey, Jonathan E Taylor, and David Siegmund.  Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(1):187–205, 2004.
  * [STW+09] Xiaojian Shao, Yingjie Tian, Lingyun Wu, Yong Wang, Ling Jing, and Naiyang Deng.  Predicting dna-and rna-binding proteins from sequences with kernel methods.  Journal of Theoretical Biology, 258(2):289–293, 2009.
  * [SZL+07] Juwen Shen, Jian Zhang, Xiaomin Luo, Weiliang Zhu, Kunqian Yu, Kaixian Chen, Yixue Li, and Hualiang Jiang.  Predicting protein–protein interactions based only on sequences information.  Proceedings of the National Academy of Sciences, 104(11):4337–4341, 2007.
  * [TFBCR19] Ryan J Tibshirani, Rina Foygel Barber, Emmanuel Candes, and Aaditya Ramdas.  Conformal prediction under covariate shift.  Advances in neural information processing systems, 32, 2019.
  * [VCC+19] Jessica Vamathevan, Dominic Clark, Paul Czodrowski, Ian Dunham, Edgardo Ferran, George Lee, Bin Li, Anant Madabhushi, Parantu Shah, Michaela Spitzer, et al.  Applications of machine learning in drug discovery and development.  Nature reviews Drug discovery, 18(6):463–477, 2019.
  * [VGS05] Vladimir Vovk, Alexander Gammerman, and Glenn Shafer.  Algorithmic learning in a random world.  Springer Science & Business Media, 2005.
  * [VW21] Vladimir Vovk and Ruodu Wang.  E-values: Calibration, combination and applications.  The Annals of Statistics, 49(3):1736–1754, 2021.
  * [WJR22] Lequn Wang, Thorsten Joachims, and Manuel Gomez Rodriguez.  Improving screening processes via calibrated subset selection.  arXiv preprint arXiv:2202.01147, 2022.
  * [WR20] Ruodu Wang and Aaditya Ramdas.  False discovery rate control with e-values.  arXiv preprint arXiv:2009.02824, 2020.

##  Appendix A Deferred results and details

###  A.1 Recap: Conformalied Selection

In this section, we briefly recall the conformalized selection framework of
[JC22]. As before, set $\mathcal{D}_{\textnormal{calib}}=\\{1,\dots,n\\}$ and
$\mathcal{D}_{\textnormal{test}}=\\{n+1,\dots,n+m\\}$. Screening for large
outcomes can be viewed as testing the null hypotheses

$\displaystyle H_{j}\colon Y_{n+j}\leq c_{n+j},\quad j=1,\dots,m,$

where selecting $j\in\mathcal{R}$ is equivalent to rejecting $H_{j}$. Perhaps
unconventionally, whether $H_{j}$ is true or not is here a random event. Yet,
we still construct _conformal_ p-values built upon the (split) conformal
inference framework [VGS05] to test these hypotheses.

Given a monotone score $V$ and thresholds $\\{c_{n+j}\\}_{j=1}^{m}$, we
compute $V_{i}=V(X_{i},Y_{i})$ for $i\in\mathcal{D}_{\textnormal{calib}}$,
$\widehat{V}_{n+j}=V(X_{n+j},c_{n+j})$, and then apply the Benjamini-Hochberg
(BH) procedure [BH95] to conformal p-values

$\displaystyle
p_{j}^{\textnormal{nw}}=\frac{\sum_{i=1}^{n}\operatorname{\mathds{1}}{\\{V_{i}<\widehat{V}_{n+j}\\}}+U_{j}(1+\sum_{i=1}^{n}\operatorname{\mathds{1}}{\\{V_{i}=\widehat{V}_{n+j}\\}})}{n+1},\quad
U_{j}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}\textrm{Unif}([0,1]).$
(14)

Under the assumption that $Q=P$, or $w(x)\equiv 1$ in (1),  each
$p_{j}^{\textnormal{nw}}$ is a valid p-value obeying (3). However, as shown
earlier, under non-trivial covariate shift, the random variables
$p_{j}^{\textnormal{nw}}$ are no longer valid p-values.555 A side result in
Appendix A.3 shows that running the BH procedure with
$\\{p_{j}^{\textnormal{nw}}\\}$ controls a weighted notion of FDR below the
target level $q$; this weighted notion downweights the error for $(X,Y)$
values that are less frequent in the calibration data.

###  A.2 Connection of our p-values and weighted conformal prediction

######  Remark A.1.

An intuitive interpretation for the p-values in (4) is that they are roughly
the critical confidence level when inverting weighted conformal prediction
intervals at the thresholds $\\{c_{n+j}\\}$. Indeed, the split weighted
conformal inference procedure [TFBCR19] (using $-V$ as the nonconformity score
for consistency) yields the one-sided prediction interval
$\widehat{C}(X_{n+j},1-\alpha)=[\eta(X_{n+j},1-\alpha),+\infty)$ for
$Y_{n+j}$, where

$\displaystyle-\eta(x,1-\alpha)=\textrm{Quantile}\big{(}1-\alpha;{\textstyle\sum_{i=1}^{n}}w_{i}(x)\delta_{-V_{i}}+w_{n+1}(x)\delta_{-\infty}\big{)},$

and $w_{i}(x)=\frac{w(X_{i})}{\sum_{i=1}^{n}w(X_{i})+w(x)}$, $i=1,\dots,n$,
and $w_{n+1}(x)=\frac{w(x)}{\sum_{i=1}^{n}w(X_{i})+w(x)}$ for any
$x\in\mathcal{X}$. Ignoring the tie-breaking random variables, we see that
$p_{j}$ is the smallest $\alpha$ such that $c_{n+j}<\eta(X_{n+j},1-\alpha)$.

###  A.3 Weighted FDR control for unweighted p-values

The following theorem shows the weighted FDR control when running BH($q$) with
unweighted conformal p-values while there is a covariate shift between
calibration and test samples.

######  Theorem A.2.

Suppose $V$ is monotone, $\\{(X_{i},Y_{i})\\}_{i=1}^{n}\sim P$ and
$\\{(X_{n+j},Y_{n+j})\\}_{j=1}^{m}\sim Q$ are mutually independent and obey
(1). Suppose
$\\{(X_{i},Y_{i})\\}_{i=1}^{n}\cup\\{(X_{n+\ell},c_{\ell})\\}_{\ell\neq
j}\cup\\{(X_{n+j},Y_{n+j})\\}$ are mutually independent for all $j$. For any
$q\in(0,1)$, let $\mathcal{R}_{\textnormal{nw}}$ be the rejection set of
BH($q$) procedure applied to the unweighted conformal $p$-values
$\\{p_{j}^{\textnormal{nw}}\\}$ in (14). Then

$\displaystyle\mathbb{E}\Bigg{[}\frac{\sum_{j=1}^{m}\operatorname{\mathds{1}}\\{j\in\mathcal{R}_{\textnormal{nw}},Y_{n+j}\leq
c_{n+j}\\}/w(X_{n+j})}{1\vee|\mathcal{R}_{\textnormal{nw}}|}\Bigg{]}\leq
q\cdot\mathbb{E}\bigg{[}\frac{n+1}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+1})}\bigg{]},$

where the expectation is over both the test and calibration data.

###  A.4 General theorem for asymptotic FDR control with BH procedure

We provide the general results for Theorem 2.5 below.

######  Theorem A.3.

Suppose $w(\cdot)$ is uniformly bounded by a fixed constant, the covariate
shift (1) holds, and $\\{c_{n+j}\\}_{j=1}^{m}$ are i.i.d. random variables.
Fix any $q\in(0,1)$, and let $\mathcal{R}$ be the output of the BH($q$)
procedure applied to $\\{p_{j}\\}_{j=1}^{m}$ in (4). Then the following two
results hold.

  1. (i)

For any fixed $m$, as $n\to\infty$, it holds that
$\limsup_{n\to\infty}\textnormal{FDR}\leq q$.

  2. (ii)

As $m,n\to\infty$, it holds that $\limsup_{m,n\to\infty}\textnormal{FDR}\leq
q$ if there exists some $t^{\prime}\in(0,1]$ such that
$G_{\infty}(t^{\prime})<q$, where
$G_{\infty}(t):=t/\mathbb{P}(F(V(X_{n+j},C_{n+j}),U_{j})\leq t)$ for any
$t\in[0,1]$, with the convention that $0/0=0$ and $t/0=\infty$ for $t>0$.
Here, for any $v\in\mathbb{R}$ and $u\in[0,1]$,

$\displaystyle F(v,u)=Q\big{(}V(X,Y)<v\big{)}+u\cdot Q\big{(}V(X,Y)=v\big{)},$
(15)

where the probability is taken with respect to the test distribution
$(X,Y)\sim Q$ obeying (1).

Furthermore, define $t^{*}=\sup\\{t\in[0,1]\colon G_{\infty}(t)\leq q\\}$.
Additionally suppose for any sufficiently small $\epsilon>0$, there exists
some $t\in(t^{*}-\epsilon,t)$ such that $G_{\infty}(t)<q$. Then the asymptotic
FDR of BH($q$) is $\frac{Q\\{F(V(X,C),U)\leq t^{*},Y\leq
C\\}}{Q\\{F(V(X,C),U)\leq t^{*}\\}}$, and the asymptotic power is
$Q\big{\\{}F(V(X,C),U)\leq t^{*}{\,|\,}Y>C\big{\\}}$, where all the
probabilities are induced by the test distribution $Q_{X,Y,C}$ and an
independent $U\sim$ Unif([0,1]).

###  A.5 Deferred simulation results for DTI prediction

Figure 15: Empirical FDR for DTI prediction with $q_{\textrm{pop}}=0.7$.
Details are the same as Figure 4. Figure 16: Empirical power for DTI with
$q_{\textrm{pop}}=0.7$. Details are the same as Figure 4.

##  Appendix B Proofs for weighted p-values

###  B.1 Proof of Lemma 2.2

###### Proof of Lemma 2.2.

For any $j\in\\{1,\dots,m\\}$, note that $V(X_{n+j},y)$ is increasing in $y$,
thus on the event $\\{Y_{n+j}\leq c_{n+j}\\}$, it holds that
$V_{n+j}\leq\widehat{V}_{n+j}$. We now define

$\displaystyle
p_{j}^{*}=\frac{\sum_{i=1}^{n}w(X_{i})\operatorname{\mathds{1}}{\\{V_{i}<{V}_{n+j}\\}}+(w(X_{n+j})+\sum_{i=1}^{n}w(X_{i})\operatorname{\mathds{1}}{\\{V_{i}={V}_{n+j}\\}})\cdot
U_{j}}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})},$

where $U_{j}$ is the same tie-breaking random variable as in $p_{j}$. That is,
$p_{j}^{*}$ is the oracle p-value we would get if we replace
$\widehat{V}_{n+j}$ by $V_{n+j}$ in the construction of $p_{j}$. Note that
$p_{j}^{*}$ is not computable; it is only for the purpose of proof. Thus, it
is clear that $p_{j}^{*}\leq p_{j}$ on the event $\\{Y_{n+j}\leq c_{n+j}\\}$.
This implies

$\displaystyle\mathbb{P}(p_{j}\leq t,~{}Y_{n+j}\leq
c_{n+j})\leq\mathbb{P}(p_{j}\leq t,~{}p_{j}^{*}\leq
p_{j})\leq\mathbb{P}(p_{j}^{*}\leq t).$

The right-handed side above is below $t$ using weighted exchangeability
introduced in [TFBCR19]. Formally, define the unordered set
$Z=[Z_{1},\dots,Z_{n},Z_{n+j}]$ for $Z_{i}=(X_{i},Y_{i})$, $i=1,\dots,n,n+j$,
and write $z$ as the realized value of $Z$, and $\mathcal{E}_{z}=\\{Z=z\\}$.
Under the covariate shift (1), the scores $\\{V_{1},\dots,V_{n},V_{n+j}\\}$
are weighted exchangeable, in the sense that

$\displaystyle
Z_{n+j}{\,|\,}\mathcal{E}_{z}~{}\sim~{}\sum_{i\in\\{1,\dots,n,n+j\\}}\delta_{z_{i}}\cdot\frac{w(x_{i})}{\sum_{i=1}^{n}w(x_{i})+w(x_{n+j})},$

i.e., given that the unordered set $Z$ equals $z$, the probability that
$Z_{n+j}$ taking on any value $(x,y)\in\\{z_{1},\dots,z_{n},z_{n+j}\\}$ is
proportional to the corresponding weight $w(x)$. This implies

$\displaystyle\mathbb{P}(p_{j}^{*}\leq t{\,|\,}\mathcal{E}_{z})=t,$

and thus $\mathbb{P}(p_{j}^{*}\leq t)=t$, which concludes the proof of Lemma
2.2. ∎

###  B.2 Proof of Theorem A.2

###### Proof of Theorem A.2.

For ease of illustration, we prove the result when there is no tie among all
nonconformity scores; similar results can be obtained with ties. We first use
a leave-one-out trick similar to [JC22]. Define the oracle conformal p-values

$\displaystyle\bar{p}_{j}^{\textnormal{nw}}=\frac{\sum_{i=1}^{n}\operatorname{\mathds{1}}{\\{V_{i}<{V}_{n+j}\\}}+U_{j}(1+\sum_{i=1}^{n}\operatorname{\mathds{1}}{\\{V_{i}={V}_{n+j}\\}})}{n+1},$

which is obtained by replacing the $\widehat{V}_{n+j}$ in (14) by
$V_{n+j}=V(X_{n+j},Y_{n+j})$. By the monotonicity of $V$, one has
$p_{j}^{\textnormal{nw}}\geq\bar{p}_{j}^{\textnormal{nw}}$ for any $j$ such
that $Y_{n+j}\leq c_{j}$. In the BH procedure, for any
$j\in\mathcal{R}_{\textnormal{nw}}$, replacing $p_{j}$ with a smaller value
does not change the rejection set. Therefore, letting $\bar{\mathcal{R}}_{j}$
be the rejection set of BH($q$) applied to
$\\{p_{1}^{\textnormal{nw}},\dots,p_{j-1}^{\textnormal{nw}},\bar{p}_{j}^{\textnormal{nw}},\dots,p_{m}^{\textnormal{nw}}\\}$,
we know that

$\displaystyle j\in\mathcal{R}_{\textnormal{nw}}~{}\textrm{and}~{}Y_{n+j}\leq
c_{n+j}\quad\Rightarrow\quad
j\in\bar{\mathcal{R}}_{j}~{}\textrm{and}~{}\bar{\mathcal{R}}_{j}=\mathcal{R}_{\textnormal{nw}}.$

We then have

$\displaystyle\mathbb{E}\Bigg{[}\frac{\sum_{j=1}^{m}\operatorname{\mathds{1}}\\{j\in\mathcal{R}_{\textnormal{nw}},Y_{n+j}\leq
c_{n+j}\\}/w(X_{n+j})}{1\vee|\mathcal{R}_{\textnormal{nw}}|}\Bigg{]}\leq\sum_{j=1}^{m}\mathbb{E}\Bigg{[}\frac{\operatorname{\mathds{1}}\\{j\in\bar{\mathcal{R}}_{j}\\}/w(X_{n+j})}{1\vee|\bar{\mathcal{R}}_{j}|}\Bigg{]}.$
(16)

Also note that $j\in\bar{\mathcal{R}}_{j}$ if and only if
$\bar{p}_{j}^{\textnormal{nw}}\leq q|\bar{\mathcal{R}}_{j}|/m$. Thus, for any
$j=1,\dots,m$,

$\displaystyle\mathbb{E}\Bigg{[}\frac{\operatorname{\mathds{1}}\\{j\in\bar{\mathcal{R}}_{j}\\}/w(X_{n+j})}{1\vee|\bar{\mathcal{R}}_{j}|}\Bigg{]}=\sum_{k=1}^{m}\frac{1}{k}\mathbb{E}\big{[}\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|=k\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}/w(X_{n+j})\big{]}$
$\displaystyle=\sum_{k=1}^{m}\frac{1}{k}\mathbb{E}\big{[}(\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}-\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k-1\\}})\cdot\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}/w(X_{n+j})\big{]}$
$\displaystyle\leq\frac{1}{m}\mathbb{E}\bigg{[}\frac{\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
q\\}}}{w(X_{n+j})}\bigg{]}+\sum_{k=1}^{m-1}\mathbb{E}\bigg{[}\frac{\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq\frac{qk}{m}\\}}}{k\cdot
w(X_{n+j})}-\frac{\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq\frac{q(k+1)}{m}\\}}}{(k+1)\cdot
w(X_{n+j})}\bigg{]}.$ (17)

We define the event $\mathcal{E}_{z}=\\{Z=z\\}$, where
$Z=[Z_{1},\dots,Z_{n},Z_{n+j}]$ is the unordered set for $Z_{i}=(X_{i},Y_{i})$
and $Z_{n+j}=(X_{n+j},Y_{n+j})$, and $z=[z_{1},\dots,z_{n},z_{n+j}]$ is the
unordered set of their realized values. By the covariate shift (1),
$\\{(X_{i},Y_{i})\\}_{i=1}^{n}$ and $(X_{n+j},Y_{n+j})$ are weighted
exchangeable [TFBCR19], hence

$\displaystyle
Z_{n+j}{\,|\,}\mathcal{E}_{z}~{}\sim~{}\sum_{i\in\\{1,\dots,n,n+j\\}}\delta_{z_{i}}\cdot\frac{w(x_{i})}{\sum_{i=1}^{n}w(x_{i})+w(x_{n+j})},$
(18)

where $\delta_{z_{i}}$ is the point mass on $z_{i}$. Therefore, by the tower
property,

$\displaystyle\mathbb{E}\Bigg{[}\frac{\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}}{k\cdot
w(X_{n+j})}-\frac{\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
q(k+1)/m\\}}}{(k+1)\cdot w(X_{n+j})}\,\bigg{|}\,\mathcal{E}_{z}\Bigg{]}$
$\displaystyle\leq\mathbb{E}\Bigg{[}\frac{\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}}{w(X_{n+j})}\cdot\frac{1}{k}\cdot\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,\bar{p}_{j}^{\textnormal{nw}},\mathcal{E}_{z}\big{)}\,\bigg{|}\,\mathcal{E}_{z}\Bigg{]}$
$\displaystyle\quad-\mathbb{E}\Bigg{[}\frac{\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
q(k+1)/m\\}}}{w(X_{n+j})}\cdot\frac{1}{k+1}\cdot\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,\bar{p}_{j}^{\textnormal{nw}},\mathcal{E}_{z}\big{)}\,\bigg{|}\,\mathcal{E}_{z}\Bigg{]}$

Note that the randomness in $U_{j}$ is independent of everything else, hence
conditioning on $(\bar{p}_{j}^{\textnormal{nw}},\mathcal{E}_{z})$ is the same
as conditioning on $(Z_{n+j},\mathcal{E}_{z})$. We also denote
$[1],[2],\dots,[n+1]$ as a permutation of $(1,2,\dots,n,n+j)$ such that
$V(X_{[1]},Y_{[1]})\leq V(X_{[2]},Y_{[2]})\leq\cdots\leq
V(X_{[n+1]},Y_{[n+1]})$. Conditional on $\mathcal{E}_{z}$, on the event
$\\{\bar{p}_{j}^{\textnormal{nw}}\leq qk/m\\}$, $Z_{n+j}$ can only take values
at $z_{i}$ for which $\frac{\sum_{\ell\neq
i}\operatorname{\mathds{1}}{\\{V(z_{\ell})<V(z_{k})\\}}}{n+1}\leq qk/m$; that
is, the possible values of $Z_{n+j}$ are in the ordered set
$\\{z_{[1]},\dots,z_{[i_{k}]}\\}$ for some $i_{k}$. Similarly, $Z_{n+j}$ can
only take values in $\\{z_{[1]},\dots,z_{[i_{k+1}]}\\}$ for $i_{k+1}\geq
i_{k}$ on the event $\\{\bar{p}_{j}^{\textnormal{nw}}\leq q(k+1)/m\\}$.
Therefore, we have

$\displaystyle\mathbb{E}\Bigg{[}\frac{\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}}{k\cdot
w(X_{n+j})}-\frac{\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
q(k+1)/m\\}}}{(k+1)\cdot w(X_{n+j})}\,\bigg{|}\,\mathcal{E}_{z}\Bigg{]}$
$\displaystyle\leq\mathbb{E}\Bigg{[}\sum_{\ell=1}^{i_{k}}\operatorname{\mathds{1}}{\\{Z_{n+j}=z_{[\ell]}\\}}\frac{\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}}{w(x_{[\ell]})}\cdot\frac{1}{k}\cdot\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[\ell]},\mathcal{E}_{z}\big{)}\,\bigg{|}\,\mathcal{E}_{z}\Bigg{]}$
$\displaystyle\quad-\mathbb{E}\Bigg{[}\sum_{\ell=1}^{i_{k+1}}\operatorname{\mathds{1}}{\\{Z_{n+j}=z_{[\ell]}\\}}\frac{\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
q(k+1)/m\\}}}{w(x_{[\ell]})}\cdot\frac{1}{k+1}\cdot\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[\ell]},\mathcal{E}_{z}\big{)}\,\bigg{|}\,\mathcal{E}_{z}\Bigg{]}$
$\displaystyle=\sum_{\ell=1}^{i_{k}}\frac{\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[\ell]},\mathcal{E}_{z}\big{)}}{k\cdot
w(x_{[\ell]})}\mathbb{E}\big{[}\operatorname{\mathds{1}}{\\{Z_{n+j}=z_{[\ell]}\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}\,\big{|}\,\mathcal{E}_{z}\big{]}$
$\displaystyle\quad-\sum_{\ell=1}^{i_{k+1}}\frac{\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[\ell]},\mathcal{E}_{z}\big{)}}{(k+1)\cdot
w(x_{[\ell]})}\mathbb{E}\big{[}\operatorname{\mathds{1}}{\\{Z_{n+j}=z_{[\ell]}\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
q(k+1)/m\\}}\,\big{|}\,\mathcal{E}_{z}\big{]}.$

By weighted exchangeability, for $i<i_{k}$ it holds deterministically that
$\bar{p}_{j}^{\textnormal{nw}}\leq qk/m$ and
$\mathbb{P}(Z_{n+j}=z_{[\ell]}{\,|\,}\mathcal{E}_{z})=\frac{w(x_{[\ell]})}{\sum_{i=1}^{n}w(x_{i})+w(x_{n+j})}$.
The randomness of $U_{j}\sim\textrm{Unif}(0,1)$ implies

$\displaystyle\mathbb{E}\big{[}\operatorname{\mathds{1}}{\\{Z_{n+j}=z_{[i_{k}]}\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}\,\big{|}\,\mathcal{E}_{z}\big{]}=(n+1)\cdot\bigg{(}\frac{qk}{m}-\frac{i_{k}-1}{n+1}\bigg{)}\cdot\frac{x_{[i_{k}]}}{\sum_{i=1}^{n}w(x_{i})+w(x_{n+j})}.$

As a result, we have

$\displaystyle\sum_{\ell=1}^{i_{k}}\frac{\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[\ell]},\mathcal{E}_{z}\big{)}}{k\cdot
w(x_{[\ell]})}\mathbb{E}\big{[}\operatorname{\mathds{1}}{\\{Z_{n+j}=z_{[\ell]}\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}\,\big{|}\,\mathcal{E}_{z}\big{]}$
$\displaystyle=\sum_{\ell=1}^{i_{k}-1}\frac{\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[\ell]},\mathcal{E}_{z}\big{)}}{k\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}+\frac{(n+1)\frac{qk}{m}-i_{k}+1}{k\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}\cdot\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[i_{k}]},\mathcal{E}_{z}\big{)}.$

Similar computation for $k+1$ leads to

$\displaystyle\mathbb{E}\Bigg{[}\frac{\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
qk/m\\}}}{k\cdot
w(X_{n+j})}-\frac{\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
q(k+1)/m\\}}}{(k+1)\cdot w(X_{n+j})}\,\bigg{|}\,\mathcal{E}_{z}\Bigg{]}$
$\displaystyle\leq\sum_{\ell=1}^{i_{k}-1}\frac{\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[\ell]},\mathcal{E}_{z}\big{)}}{k\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}-\sum_{\ell=1}^{i_{k+1}-1}\frac{\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[\ell]},\mathcal{E}_{z}\big{)}}{(k+1)\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}$
$\displaystyle\quad-\frac{(n+1)\frac{qk}{m}-i_{k}+1}{k\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}\cdot\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[i_{k}]},\mathcal{E}_{z}\big{)}$
$\displaystyle\quad-\frac{(n+1)\frac{q(k+1)}{m}-i_{k+1}+1}{(k+1)\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}\cdot\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[i_{k+1}]},\mathcal{E}_{z}\big{)}.$

Conditional on $\mathcal{E}_{z}$,
$\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[\ell]},\mathcal{E}_{z}\big{)}$ is increasing in
$\ell$; to see this, simply note that the p-values
$(p_{1}^{\textnormal{nw}},\dots,p_{j-1}^{\textnormal{nw}},\bar{p}_{j}^{\textnormal{nw}},p_{j+1}^{\textnormal{nw}},\dots,p_{m}^{\textnormal{nw}})$
are coordinate-wisely increasing in $\ell$, hence the size of the rejection
set $\bar{\mathcal{R}}_{j}$ is decreasing in $\ell$. Thus, replacing the index
$\ell$ in the two summations by $i_{k}$, we have

$\displaystyle\mathbb{E}\Bigg{[}\frac{\operatorname{\mathds{1}}{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}_{\\{p_{j}^{*}\leq qk/m\\}}}{k\cdot
w(X_{n+j})}-\frac{\operatorname{\mathds{1}}_{\\{|\bar{\mathcal{R}}_{j}|\leq
k\\}}\operatorname{\mathds{1}}_{\\{p_{j}^{*}\leq q(k+1)/m\\}}}{(k+1)\cdot
w(X_{n+j})}\,\bigg{|}\,\mathcal{E}_{z}\Bigg{]}$
$\displaystyle\leq\sum_{\ell=1}^{i_{k}-1}\frac{\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[i_{k}]},\mathcal{E}_{z}\big{)}}{k\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}-\sum_{\ell=1}^{i_{k+1}-1}\frac{\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[i_{k}]},\mathcal{E}_{z}\big{)}}{(k+1)\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}$
$\displaystyle\quad+\frac{(n+1)\frac{qk}{m}-i_{k}+1}{k\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[i_{k}]},\mathcal{E}_{z}\big{)}$
$\displaystyle\quad+\frac{(n+1)\frac{q(k+1)}{m}-i_{k+1}+1}{(k+1)\cdot(\sum_{i=1}^{n}w(x_{i})+w(x_{n+j}))}\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[i_{k}]},\mathcal{E}_{z}\big{)}$
$\displaystyle=\frac{\mathbb{P}\big{(}|\bar{\mathcal{R}}_{j}|\leq
k\,\big{|}\,Z_{n+j}=z_{[i_{k}]},\mathcal{E}_{z}\big{)}}{\sum_{i=1}^{n}w(x_{i})+w(x_{n+j})}\Big{(}\frac{(n+1)qk/m}{k}-\frac{(n+1)q(k+1)/m}{k+1}\Big{)}=0.$

For any $j$, the first term in (B.2) can be calculated as

$\displaystyle\mathbb{E}\big{[}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
q\\}}/w(X_{n+j})\,\big{|}\,\mathcal{E}_{z}\big{]}$
$\displaystyle=\sum_{\ell=1}^{i_{m}}\frac{\mathbb{P}(\bar{p}_{j}^{\textnormal{nw}}\leq
q{\,|\,}Z_{n+j}=z_{[\ell]},\mathcal{E}_{z})\cdot\mathbb{P}(Z_{n+j}=z_{[\ell]}{\,|\,}\mathcal{E}_{z})}{w(x_{[\ell]})}$
$\displaystyle=\frac{i_{m}-1+(n+1)q-i_{m}+1}{\sum_{i=1}^{n}w(x_{i})+w(x_{n+j})}=\frac{(n+1)q}{\sum_{i=1}^{n}w(x_{i})+w(x_{n+j})},$

where the second equality uses (18), and the definition of
$\bar{p}_{j}^{\textnormal{nw}}$. Returning to (16) and (B.2), we obtain

$\displaystyle\mathbb{E}\Bigg{[}\frac{\sum_{j=1}^{m}\operatorname{\mathds{1}}\\{j\in\mathcal{R}_{\textnormal{nw}},Y_{n+j}\leq
c_{n+j}\\}/w(X_{n+j})}{1\vee|\mathcal{R}^{\textnormal{nw}}|}\Bigg{]}$
$\displaystyle\leq\frac{1}{m}\sum_{j=1}^{m}\mathbb{E}\big{[}\operatorname{\mathds{1}}{\\{\bar{p}_{j}^{\textnormal{nw}}\leq
q\\}}/w(X_{n+j})\big{]}=\mathbb{E}\bigg{[}\frac{(n+1)q}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+1})}\bigg{]},$

where the last equality uses the fact that the covariates of test data are
i.i.d. This concludes the proof of Theorem A.2. ∎

###  B.3 Proof of Non-PRDS

###### Proof of Proposition 2.4.

Without loss of generality, we consider the oracle p-values

$\displaystyle
p_{j}^{*}=\frac{\sum_{i=1}^{n}w(X_{i})\operatorname{\mathds{1}}{\\{V_{i}<{V}_{n+j}\\}}+(w(X_{n+j})+\sum_{i=1}^{n}w(X_{i})\operatorname{\mathds{1}}{\\{V_{i}={V}_{n+j}\\}})\cdot
U_{j}}{\sum_{i=1}^{n}w(X_{i})+w(X_{n+j})},$ (19)

which coincide with (4) when $Y_{n+j}=c_{n+j}$.

We prove by constructing a counterexample and computing the distribution
functions of p-values to show the violation of the PRDS property. We set
$w(x)=x$ for $P=\textrm{Unif}([1/2,3/2])$, let $Q$ have density function
$q(x)=x\operatorname{\mathds{1}}\\{1/2\leq x\leq 3/2\\}$, and $V(x,y)=-w(x)$.
Take $n=2$ and any arbitrary $m\geq 2$. If the weighted conformal p-values are
PRDS, we would have

$\displaystyle F(s,t):=\mathbb{P}(p_{1}\leq s{\,|\,}p_{2}\geq t)$

is decreasing in $t\in(0,1)$ for any fixed $s\in(0,1)$. However, we will show
$F(\frac{1}{10},\frac{3}{5})<F(\frac{1}{10},\frac{9}{10})$. We will compute
the two values exactly; they can also be approximated to arbitrary accuracy
via Monte Carlo.

For notational simplicity, we write $W_{i}=w(X_{i})$ for $i=1,2$, and
$\widetilde{W}_{j}=w(X_{2+j})$ for $j=1,2$. By construction,
$W_{1},W_{2},\widetilde{W}_{1},\widetilde{W}_{2}$ are mutually independent
with $W_{1},W_{2}\sim P$, $\widetilde{W}_{1},\widetilde{W}_{2}\sim Q$. As $P$
and $Q$ do not have point mass, there are no ties among the four scores almost
surely. By the symmetry of $W_{1},W_{2}$, for any $s,t\in(0,1)$, we have

$\displaystyle F(s,t)=(1-t)^{-1}\cdot\mathbb{P}(p_{1}\leq s,p_{2}\geq
t)=2(1-t)^{-1}\cdot\mathbb{P}(p_{1}\leq s,p_{2}\geq t,W_{1}<W_{2}).$

Since $v(x,y)=-w(x)$, the two weighted conformal p-values are equivalently

$\displaystyle
p_{i}=\frac{W_{1}\operatorname{\mathds{1}}\\{W_{1}>\widetilde{W}_{i}\\}+W_{2}\operatorname{\mathds{1}}\\{W_{2}>\widetilde{W}_{i}\\}+U_{i}\widetilde{W}_{i}}{W_{1}+W_{2}+\widetilde{W}_{i}},\quad
i=1,2.$

When $W_{1}<W_{2}$, the form of $p_{i}$ can be divided into three cases:

  1. (i)

When $\widetilde{W}_{i}<W_{1}<W_{2}$, one has
$p_{i}=\frac{W_{1}+W_{2}+U_{i}\widetilde{W}_{i}}{W_{1}+W_{2}+\widetilde{W}_{i}}$,
and $p_{i}\in[2/3,1]$ holds deterministically.

  2. (ii)

When $W_{1}<\widetilde{W}_{i}<W_{2}$, one has
$p_{i}=\frac{W_{2}+U_{i}\widetilde{W}_{i}}{W_{1}+W_{2}+\widetilde{W}_{i}}$. As
$W_{i},\widetilde{W}_{i}\in[1/2,3/2]$, it holds that $p_{i}\in[1/3,6/7]$.

  3. (iii)

When $W_{1}<W_{2}<\widetilde{W}_{i}$, one has
$p_{i}=\frac{U_{i}\widetilde{W}_{i}}{W_{1}+W_{2}+\widetilde{W}_{i}}$, and
$p_{i}\in[0,3/5]$.

Intuitively, the reason for violating PRDS in this case is that, when
$p_{1}\leq 1/10$ is small, it must always take the form
$p_{1}=\frac{U_{1}\cdot\widetilde{W}_{1}}{W_{1}+W_{2}+\widetilde{W}_{1}}$ as
in case (iii). Meanwhile, when $p_{2}\geq 3/5$ is relatively large, it may be
as cases (i) or (ii). In that two cases, roughly, a larger $p_{2}$ means
relatively larger $W_{1},W_{2}$, hence smaller $p_{1}$.