# ZAP: $Z$-value Adaptive Procedures for False Discovery Rate Control with
Side Information

Dennis Leung  School of Mathematics and Statistics, University of Melbourne
<EMAIL_ADDRESS>and  Wenguang Sun  Center for Data Science,
Zhejiang University<EMAIL_ADDRESS>
###### Abstract.

Adaptive multiple testing with covariates is an important research direction
that has gained major attention in recent years. It has been widely recognized
that leveraging side information provided by auxiliary covariates can improve
the power of false discovery rate (FDR) procedures. Currently, most such
procedures are devised with $p$-values as their main statistics. However, for
two-sided hypotheses, the usual data processing step that transforms the
primary statistics, known as $z$-values, into $p$-values not only leads to a
loss of information carried by the main statistics, but can also undermine the
ability of the covariates to assist with the FDR inference. We develop a
$z$-value based covariate-adaptive (ZAP) methodology that operates on the
intact structural information encoded jointly by the $z$-values and
covariates. It seeks to emulate the oracle $z$-value procedure via a working
model, and its rejection regions significantly depart from those of the
$p$-value adaptive testing approaches. The key strength of ZAP is that the FDR
control is guaranteed with minimal assumptions, even when the working model is
misspecified. We demonstrate the state-of-the-art performance of ZAP using
both simulated and real data, which shows that the efficiency gain can be
substantial in comparison with p-value based methods. Our methodology is
implemented in the R package zap.

###### Key words and phrases:

multiple testing, false discovery rate, $z$-value, beta mixture, side
information

###### 2000 Mathematics Subject Classification:

62H05

##  1\. Introduction

In modern scientific studies, a ubiquitous task is to test a multitude of two-
sided hypotheses regarding the presence of nonzero effects. The problem of
multiple testing with covariates has received much recent attention, as
leveraging contextual information beyond what is offered by the main
statistics can enhance both the power and interpretability of existing false
discovery rate (FDR; Benjamini and Hochberg, 1995) methods. This has marked a
gradual paradigm shift from the Benjamini-Hochberg (BH) procedure and its
immediate variants (e.g. Benjamini and Hochberg, 2000, Storey, 2002) that are
based solely on the $p$-values. For instance, in the differential analysis of
RNA-sequencing data, the average read depths across samples can provide useful
side information alongside individual $p$-values, and incorporating such
information promises to improve the efficiency of existing methods. The
importance of this direction has been reflected by its intense research
activities; see Lei and Fithian (2018), Zhang and Chen (2020), Yurko et al.
(2020), Chen et al. (2017), Ignatiadis et al. (2016), Li and Barber (2019),
Boca and Leek (2018) for an incomplete list of related works. In contrast with
BH and its variants that apply a universal threshold to all $p$-values, these
methods boil down to setting varied $p$-value thresholds that are adaptive to
the covariate information.

This seemingly natural modus operandi, which involves using $p$-values as the
basic building blocks, however, is suboptimal for two-sided testing; the
$p$-values involved are typically formed via a data reduction step, which
applies a non-bijective transformation to “primary” test statistics such as
the $z$-values, $t$-statistics (Ritchie et al., 2015) or Wald statistics (Love
et al., 2014). Sun and Cai (2007) and Storey et al. (2007) argued that
reducing $z$-values to two-sided $p$-values may lead to substantial loss of
information, particularly when the $z$-values exhibit distributional
asymmetry. The main thrust of this article is to reveal a new source of
information loss in the context of covariate-adaptive multiple testing, and to
develop a $\bf z$-value covariate-adaptive methodology, which we call “ZAP”
for short, that bypasses the data reduction step. As illustrated in Section
2.2, the interactive relationship between the $z$-values and the covariates
can capture structural information that can be exploited for more testing
power. However, this interactive information may be undercut, and in some
scenarios, completely forgone when converting the $z$-values to $p$-values.
Hence, the data reduction step not only leads to a loss of information carried
by the main statistics, but also undermines the ability of the covariates to
assist with the FDR inference.

Few works on covariate-adaptive testing have pursued the $z$-value direction
for two-sided testing since combining the $z$-values and covariates poses an
additional layer of challenges. Existing $z$-value based procedures either
make strong assumptions on the data generating model (Scott et al., 2015), or
are not robust for handling multi-dimensional covariate data (Cai et al.,
2019). By contrast, ZAP retains the merits of $z$-value based methods and
avoids the information loss from “collapsing into $p$-values”, without relying
on strong assumptions nor forgoing robustness. It faithfully preserves the
interactive structure between the primary statistics and covariates as a
starting point for inference, and is deployed with a working model, whose
potential misspecification will not invalidate the FDR control.

Our contribution is twofold. First, ZAP represents a $z$-value based,
covariate-adaptive testing framework that attains state-of-the-art power
performance under minimal assumptions, filling an important gap in the
literature. Second, in light of a plethora of $p$-value based covariate-
adaptive methods that have emerged in recent years, our study explicates new
sources of information loss in data processing, which provides new insights
and gives caveats for conducting covariate-adaptive inference in practical
settings.

The rest of the paper is structured as follows. Section 2 states the problem
formulation and describes the high-level ideas of ZAP. Section 3 formally
introduces our two data-driven methods of ZAP and their implementation
details. Numerical results based on both simulated and real data are presented
in Section 4. Section 5 concludes the article with a discussion of open
issues.

##  2\. Problem Formulation and Basic Framework

###  2.1. The problem statement

Suppose we are interested in making inference of $m$ real-valued effects
$\mu_{i}$, $i=1,\dots,m$, and for each $i$, we observe a primary statistic
$Z_{i}\in\mathbb{R}$ (“$z$-value”) and an auxiliary covariate
$X_{i}\in\mathbb{R}^{p}$ that can be multivariate. We consider a multiple
testing problem where the goal is to identify nonzero effects or,
equivalently, determine the values of the indicators

(2.1) $H_{i}\equiv I(\mu_{i}\neq 0)=\begin{dcases*}1&if $\mu_{i}\neq 0$\\\
0&otherwise\end{dcases*}.$

Assume that the triples $\\{H_{i},Z_{i},X_{i}\\}_{i=1}^{m}$ are independent
and identically distributed, and the data are described by the following
mixture model:

(2.2) $Z_{i}|X_{i}=x\sim f_{x}(z)\equiv
f(z|x)=(1-w_{x})f_{0}(z)+w_{x}f_{1,x}(z),$

where $w_{x}\equiv P(H_{i}=1|X_{i}=x)$ is the conditional probability of
having a non-zero effect given $X_{i}=x$ and $f_{1,x}(z)\equiv
f(z|H_{i}=1,X_{i}=x)$ is the conditional density under the alternative.
$f_{0}$ denotes the null density, which is invariant to the covariate value.
In this article we assume $f_{0}(z)\equiv\phi(z)$, the density of a
$\mathcal{N}(0,1)$ variable111This can be easily achieved via the composite
transformation $\Phi^{-1}\circ G_{0}(\cdot)$ if the primary test statistic has
a known null distribution function $G_{0}(\cdot)$, e.g. a t-distribution,
where $\Phi(\cdot)$ is the standard normal distribution function.. In contrast
with Scott et al. (2015) which assumes a fixed alternative density, i.e.
$f_{1,x}\equiv f_{1}$, the data generating model in (2.2) provides a more
general framework for multiple testing with covariates by allowing both
$w_{x}$ and $f_{1,x}(z)$ to vary in $x$.

Let $\mathcal{R}\subset\\{1,\dots,m\\}$ be the set of hypotheses rejected by a
multiple testing procedure. In large-scale testing problems, the widely used
FDR is defined as

$\mbox{FDR}=\mathbb{E}\left[\frac{V}{R\vee 1}\right],$

where $V=\sum_{i=1}^{m}(1-H_{i})I(i\in\mathcal{R})\text{ and
}R=\sum_{i=1}^{m}I(i\in\mathcal{R})$ are respectively the number of false
positives and the number of rejections. Throughout, $\mathbb{E}[\cdot]$
denotes an expectation operator with respect to the joint distribution of
$\\{H_{i},Z_{i},X_{i}\\}_{i=1}^{m}$, and $\mathbb{E}[\cdot|\cdot]$ denotes a
conditional expectation that should be self-explanatory from its context. The
ratio ${V}/(R\vee 1)$ is known as the false discovery proportion (FDP). The
power of a testing procedure can be evaluated using the expected number of
true discoveries $\mbox{ETD}=\mathbb{E}[R-V]$ or the true positive rate

$\mbox{TPR}=\mathbb{E}\left[\frac{R-V}{(\sum_{i=1}^{m}H_{i})\vee 1}\right].$

Our goal is to devise a powerful procedure that can control the FDR under a
pre-specified level $\alpha\in(0,1)$.

###  2.2. Information loss in covariate-adaptive testing

A two-sided $p$-value is formed by the non-bijective transformation
$P_{i}=2\Phi(-|Z_{i}|)$, where $\Phi$ is the cumulative distribution function
of a $\mathcal{N}(0,1)$ variable. We call a testing procedure _z-value based_
if it makes rejection decisions based on the full dataset
$\\{Z_{i},X_{i}\\}_{i=1}^{m}$, and _p-value based_ if it does so only based on
the reduced dataset $\\{P_{i},X_{i}\\}_{i=1}^{m}$. This section presents
examples to illustrate that the interactive structure between $Z_{i}$ and
$X_{i}$ is generally not preserved by transforming into $p$-values; the
associated information loss can lead to decreased power in the FDR inference.

Consider Model (2.2), and suppose $X_{i}\sim\mbox{Unif}(-1,1)$. Our study
examines three situations:

  1. Example 2.1

Asymmetric alternatives:
$f(z|x)=\frac{8-x}{10}f_{0}(z)+\frac{x+2}{10}\phi(z-1.5)$.

  2. Example 2.2

Unbalanced covariate effects on the non-null proportions:

$f(z|x)=0.8f_{0}(z)+\frac{1-x}{10}\phi(z+1.5)+\frac{1+x}{10}\phi(z-1.5).$

  3. Example 2.3

Unbalanced covariate effects on the alternative means:

$f(z|x)=0.9f_{0}(z)+0.1\phi(z-1.5\;\text{sgn}(x)),\;\text{where}\;\text{sgn}(x)=I(x\geq
0)-I(x<0).$

We investigate two approaches to FDR analysis for these examples that
respectively reject hypotheses with suitably small posterior probabilities
$\\{P(H_{i}=0|P_{i},X_{i})\\}_{i=1}^{m}$ and
$\\{P(H_{i}=0|Z_{i},X_{i})\\}_{i=1}^{m}$. The latter probabilities are assumed
to be known by an oracle222In practice these posterior probabilities are
unknown.. In the literature the $z$-value based quantity
$P(H_{i}=0|Z_{i},X_{i})$ is also called the conditional local false discovery
rate (CLfdr, Efron, 2008, Cai and Sun, 2009). It is known that the _optimal_
$p$-value and $z$-value based procedures, which maximize true discoveries
subject to false discovery constraints, have the respective forms

$\boldsymbol{\delta}^{\mathcal{P}}=\bigl{\\{}I[P(H_{i}=0|P_{i},X_{i})\leq
t_{\mathcal{P}}]\bigr{\\}}_{i=1}^{m}\text{ and
}\boldsymbol{\delta}^{\mathcal{Z}}=\bigl{\\{}I[P(H_{i}=0|Z_{i},X_{i})\leq
t_{\mathcal{Z}}]\bigr{\\}}_{i=1}^{m},$

where the rejection decisions are expressed by indicators, and the thresholds
$t_{\mathcal{P}}$ and $t_{\mathcal{Z}}$ are calibrated such that the nominal
FDR level is exactly $\alpha$; see Appendix A.1 for a review. In our
comparisons we choose suitable thresholds such that the FDR of both methods is
exactly $0.1$, and their powers are reported as the TPR empirically computed
by 150 repeated experiments for $m=1000$:

  * Example 2.1: $\mbox{TPR}_{\boldsymbol{\delta}^{\mathcal{P}}}=4.4\%$; $\mbox{TPR}_{\boldsymbol{\delta}^{\mathcal{Z}}}=11.7\%$.

  * Example 2.2: $\mbox{TPR}_{\boldsymbol{\delta}^{\mathcal{P}}}=3.4\%$; $\mbox{TPR}_{\boldsymbol{\delta}^{\mathcal{Z}}}=5.5\%$.

  * Example 2.3: $\mbox{TPR}_{\boldsymbol{\delta}^{\mathcal{P}}}=0.6\%$; $\mbox{TPR}_{\boldsymbol{\delta}^{\mathcal{Z}}}=2.6\%$.

Apparently, ${\boldsymbol{\delta}^{\mathcal{Z}}}$ is more powerful than
${\boldsymbol{\delta}^{\mathcal{P}}}$.

To understand the differences in power, we first remark that either oracle
procedure essentially amounts to one by which $i$ is rejected if and only if

(2.3) $Z_{i}\in\mathcal{S}(X_{i})\subset\mathbb{R}$

for some rejection region $\mathcal{S}(\cdot)$ on the $z$-value scale that is
a function of the covariate value; the theoretical derivation is sketched in
Appendix A.2. Let $\mathcal{S}^{\mathcal{P}}(x)$ and
$\mathcal{S}^{\mathcal{Z}}(x)$ denote the respective rejection regions of
$\boldsymbol{\delta}^{\mathcal{P}}$ and $\boldsymbol{\delta}^{\mathcal{Z}}$ on
the $z$-value scale for a given covariate value $x$, which are plotted for the
three examples in Figure 2.1. On the left panel, both
$\mathcal{S}^{\mathcal{P}}(x)$ and $\mathcal{S}^{\mathcal{Z}}(x)$ enlarge as
the covariate value increases, suggesting that the covariates are informative
for both methods. The information loss leading to the lesser power of
$\boldsymbol{\delta}^{\mathcal{P}}$ in Example 2.1 is intrinsically within the
main statistics when converting $z$-values to $p$-values (Sun and Cai, 2007,
Storey et al., 2007). By contrast, the middle and right panels show that
$\mathcal{S}^{\mathcal{Z}}(x)$ changes with $x$, while
$\mathcal{S}^{\mathcal{P}}(x)$ is completely insensitive to the changes in
$x$; see Appendix A.2 for the relevant calculations. Hence, the covariates are
only informative for $\boldsymbol{\delta}^{\mathcal{Z}}$. This fundamental
phenomenon reveals that upon reduction to $p$-values, the information loss not
only can occur internally within the main statistics, but also externally due
to the failure of $\boldsymbol{\delta}^{\mathcal{P}}$ in fully capturing the
original interactive information between $Z_{i}$ and $X_{i}$. When the latter
interactive structure represents the bulk of the information provided by the
covariates for testing, reduction to $p$-values can substantially undermine
the covariates’ ability to assist with inference.

Figure 2.1.  Comparisons of rejection regions. For each covariate value
$x\in(-1,1)$, the rejection regions $\mathcal{S}^{\mathcal{P}}(x)$ and
$\mathcal{S}^{\mathcal{Z}}(x)$ are respectively marked in green and red on the
$z$-value scale. The overlapped region is depicted in yellow.

###  2.3. The ZAP framework and a preview of contributions

The previous examples motivate us to focus on $z$-value adaptive procedures to
avoid information loss. This naturally boils down to pursuing the oracle
procedure $\boldsymbol{\delta}^{\mathcal{Z}}$ in some shape or form, which
presents unique challenges. Existing $z$-value based works such as Scott et
al. (2015) and Cai et al. (2019) are built directly upon
$\\{P(H_{i}=0|Z_{i},X_{i})\\}_{i=1}^{m}$, the CLfdr statistics, which
unfortunately involves unknown quantities that can be difficult to estimate in
the presence of covariates. Commonly used algorithms may not produce desired
estimates, and even lead to invalid FDR procedures if the modelling
assumptions are violated. That the theory on FDR control critically depends on
the quality of these estimates has limited the scope and applicability of
these works.

We aim to develop a new class of $z$-value adaptive (ZAP) procedures that are
assumption-lean, robust and capable of effectively exploiting the interactive
information between $Z_{i}$ and $X_{i}$. The key idea is to emulate the oracle
procedure $\boldsymbol{\delta}^{\mathcal{Z}}$ while circumventing the direct
estimation of $P(H_{i}=0|Z_{i},X_{i})$. Next we first outline the key steps
(ranking and thresholding) of our framework and then provide a preview of its
contributions.

In the first _ranking_ step, we introduce the new concept of _assessor
functions_ , which can be estimated from the data based on a working model, to
construct a new sequence of significance indices $\\{T_{i}\\}_{i=1}^{m}$ as
proxies for $\\{P(H_{i}=0|Z_{i},X_{i})\\}_{i=1}^{m}$. While many potential
working models can be used, in this work we focus on a class of beta-mixture
models that are carefully defined on a bijective transformation of the
$z$-values and particularly suitable for two-sided testing (Section 3.2). In
the second _thresholding_ step, ZAP calibrates a threshold along the ranking
produced by $\\{T_{i}\\}_{i=1}^{m}$. The essential idea is to count the number
of false rejections by any candidate threshold value with a “mirroring”
sequence of the rejected significance indices, which can be created via either
simulation (Algorithm 1) or partial data masking (Algorithm 2). The key
strength of ZAP over the methods in Scott et al. (2015) and Cai et al. (2019)
is that it seeks to emulate the oracle $z$-value procedure while avoiding a
direct substitution of $P(H_{i}|Z_{i},X_{i})$ with its estimate. ZAP is
assumption-lean and robust in the sense that it is provably valid for FDR
control under misspecifications of the working model. We stress that the
resulting rejection regions of ZAP significantly depart from those of
$p$-value adaptive methods, including the closely related CAMT (Zhang and
Chen, 2020) and AdaPT (Lei and Fithian, 2018). Our simulation and real data
studies show that the resulting gain in testing power can be substantial.

##  3\. Data-Driven ZAP Procedures

This section develops the framework of ZAP and its data-driven algorithms for
covariate-adaptive FDR inference. Section 3.1 introduces the concept of an
assessor function and a prototype procedure inspired by the oracle $z$-value
procedure. The assessor function can be constructed based on a working beta-
mixture model, which is proposed in Section 3.2. Sections 3.3 and 3.4 lay out
two variants of data-driven ZAP procedures and establish their theoretical
properties. Further implementation details are discussed in Section 3.5.

###  3.1. Preliminaries: oracle $z$-value procedure, assessor function and a
prototype ZAP algorithm

To facilitate the development of a working model, we consider the following
lossless transformation: $U_{i}=\Phi(Z_{i})$. The transformed statistic
$U_{i}$ is referred to as a _$u$ -value_333Despite the similarity in their
constructions, the $u$-values should _not_ be treated as $p$-values for one-
sided tests, which are not the subject matter of this work., which, according
to (2.2), obeys the induced mixture model

(3.1) $U_{i}|X_{i}=x\sim h_{x}(u)\equiv
h(u|x)=(1-w_{x})h_{0}(u)+w_{x}h_{1,x}(u),$

with $h_{0}(u)$ and $h_{1,x}(u)\equiv h(u|H_{i}=1,X_{i}=x)$ respectively being
the null $\mbox{Unif}(0,1)$ and conditional alternative densities. An optimal
FDR procedure (Sun and Cai, 2007, Cai and Sun, 2009, Heller and Rosset, 2021)
is a thresholding rule based on the conditional local false discovery rates
(CLfdr)

(3.2) $\mbox{CLfdr}_{i}\equiv
P(H_{i}=0|Z_{i},X_{i})=P(H_{i}=0|U_{i},X_{i})=\frac{(1-w_{X_{i}})h_{0}(U_{i})}{h_{X_{i}}(U_{i})},i=1,\dots,m.$

Since each $\mbox{CLfdr}_{i}$ is a function of $U_{i}$ conditional on $X_{i}$,
we let $\text{CLfdr}_{x}(u):(0,1)\rightarrow(0,1)$ be the corresponding
function defined on the $u$-value scale for a given realized covariate value
$x$. Related data-driven CLfdr procedures involve first estimating the CLfdr
statistics, and second determining a threshold for them using, for example,
step-wise algorithms (Sun and Cai, 2007), randomized rules (Basu et al., 2018)
or linear programming (Heller and Rosset, 2021). However, the first estimation
step poses significant challenges as it boils down to a hard _density
regression_ problem (Dunson et al., 2007). For example, to estimate
$f_{x}(\cdot)$ (or equivalently $h_{x}(\cdot)$), a line of works (Scott et
al., 2015, Tansey et al., 2018, Deb et al., 2021) proceeds by assuming a fixed
alternative density, i.e.

(3.3) $f_{1,x}(z)\equiv f_{1}(z),$

to make way for the application of an EM algorithm. If the assumption fails to
hold, the CLfdr statistics can be poorly estimated and lead to both invalid
FDR control and adversely affected power. The non-parametric CARS procedure
developed in Cai et al. (2019) does not require the assumption in (3.3).
However, it still employs the CLfdr statistics as its basic building blocks,
which are estimated with kernel density methods. Due to the curse of
dimensionality, the methodology becomes unstable in the presence of
multivariate covariates, which has limited its applicability.

By contrast, ZAP strives to sensibly emulate the oracle procedure without
heavy reliance on the quality of the CLfdr estimates, which is its key
strength. To motivate our data-driven procedures in the next sections, we
shall first discuss a prototype ZAP procedure to illustrate two key steps of
our testing framework: (a) how to combine $Z_{i}$ (or equivalently $U_{i}$)
and $X_{i}$ for assessing the significance of hypotheses; and (b) how to
threshold the new significance indices.

Step (a) involves the construction of an _assessor function_ 444Or simply
known as an assessor. $a_{x}(u):(0,1)\rightarrow(0,1)$, which seeks to
approximate the $\text{Clfdr}_{x}(u)$ function to integrate the information in
both the $u$-value and covariate. For the present assume that $a_{x}(\cdot)$
is pre-determined. Let $T_{i}\equiv a_{X_{i}}(U_{i})$ be the new significance
index for $i$ and $c_{i}(t)\equiv P(T_{i}\leq t|H_{i}=0,X_{i})$ be its null
distribution conditional on $X_{i}$. Assume that $c_{i}(\cdot)$ is continuous
and strictly increasing555Both are true as the consequences of the way we will
construct $a_{x}(\cdot)$; see the discussion after Lemma C.2 in Appendix C.,
and denote its inverse by $c_{i}^{-1}(\cdot)$. All hypotheses will then be
ordered according to the $T_{i}$’s, with a smaller $T_{i}$ indicating a more
significant hypothesis.

In Step (b), we aim to determine a threshold for the $T_{i}$’s to control the
FDR. This involves the construction of a conservative FDP estimator for any
candidate threshold $t$ by the Barber-Candes (BC) method (Arias-Castro et al.,
2017, Barber and Candès, 2015):

(3.4) $\widehat{\mbox{FDP}}(t)\equiv\frac{1+\\#\\{i:S_{i}\geq
1-c_{i}(t)\\}}{\\#\\{i:T_{i}\leq t\\}\vee
1}=\frac{1+\\#\\{i:T^{\mathfrak{m}}_{i}\leq t\\}}{\\#\\{i:T_{i}\leq t\\}\vee
1},$

where, given that $i$ is a true null, $S_{i}\equiv c_{i}(T_{i})$ is the
probability of realizing a smaller significance index and
$T^{\mathfrak{m}}_{i}\equiv c_{i}^{-1}(1-S_{i})$ is the _mirror statistic_
that “reflects” $T_{i}$’s position in the distribution $c_{i}$. Define

(3.5)
$\hat{t}(\alpha)\equiv\max\\{t\in(0,t_{\max}]:\widehat{\mbox{FDP}}(t)\leq\alpha\\},$

where $t_{\max}\equiv\max\\{t:c_{i}(t)\leq 0.5\text{ for all }i\\}$. It
follows from Barber and Candès (2019, Lemma 1) that a procedure which rejects
$i$ whenever $T_{i}\leq\hat{t}(\alpha)$ controls the FDR at level $\alpha$;
see Appendix B. Importantly, the FDR is controlled under the desired level
$\alpha$ whether $a_{x}(\cdot)$ is a good approximation of
$\text{CLfdr}_{x}(\cdot)$ or not.

However, the assessor $a_{x}(\cdot)$, which is taken as pre-determined thus
far, is to be estimated from the observed data in practice. This leads to
additional difficulties in both methodological and theoretical developments;
for one thing, the theory in Barber and Candès (2019) cannot be directly
applied to prove the FDR controlling property. Section 3.2 discusses a working
beta-mixture model, whose parameters can be estimated from the observed data
and subsequently used to construct a data-driven assessor
$\hat{a}_{x}(\cdot)$. From there we can test the hypotheses in a data-driven
manner by either implementing the prototype procedure directly using
$\hat{a}_{x}(\cdot)$ as if it is pre-determined (Section 3.3), or mimic the
prototype procedure in a more nuanced manner by leveraging the partial data
masking technique in Lei and Fithian (2018) (Section 3.4). These two variants
of ZAP have their own relative strengths and weaknesses: the direct approach
offers asymptotic FDR control under suitable regularity conditions, and is
both computationally and power efficient, while the data masking approach
offers finite-sample FDR control but is computationally intensive and
moderately less powerful in practice; see Appendix H.6 for a discussion on
aspects of their computational costs.

###  3.2. A beta-mixture model

We now develop a working model to approximate (3.1), which will be
subsequently used to construct the assessor. We propose to capture the overall
shape of $h_{x}(u)$ using a three-component mixture:

(3.6)
$h_{x}(u)=(1-\pi_{l,x}-\pi_{r,x})h_{0}(u)+\pi_{l,x}h_{l,x}(u)+\pi_{r,x}h_{r,x}(u),$

where, given $X_{i}=x$, $\pi_{l,x}$ and $\pi_{r,x}$ respectively denote the
mixing probabilities that $\mu_{i}<0$ and $\mu_{i}>0$ 666The different symbols
$\pi_{l,x}$ and $\pi_{r,x}$ are used in the working model. In the true data
generating model (3.1), the mixing probability is denoted $w_{x}$., and
$h_{l,x}$ and $h_{r,x}$ respectively represent the densities of the negative
and positive effects (on the left and right sides of the null). Our working
model assumes that $\pi_{l,x}$ and $\pi_{r,x}$ are multinomial probabilities
with regression parameter vectors $\theta_{l}$ and $\theta_{r}$:

$\pi_{l,x}=\frac{\exp(\tilde{x}^{T}\theta_{l})}{1+\exp(\tilde{x}^{T}\theta_{r})+\exp(\tilde{x}^{T}\theta_{l})},\;\pi_{r,x}=\frac{\exp(\tilde{x}^{T}\theta_{r})}{1+\exp(\tilde{x}^{T}\theta_{r})+\exp(\tilde{x}^{T}\theta_{l})},$

where $\tilde{x}=(1,x^{T})^{T}$ is the intercept-augmented covariate vector.
Further, $h_{l,x}$ and $h_{r,x}$ are chosen to be beta densities with
regression parameters $\beta_{l}$ and $\beta_{r}$:

$h_{l,x}(u)=\frac{1}{B(k_{l,x},\gamma_{l})}u^{k_{l,x}-1}(1-u)^{\gamma_{l}-1},h_{r,x}(u)=\frac{1}{B(\gamma_{r},k_{r,x})}u^{\gamma_{r}-1}(1-u)^{k_{r,x}-1},$

where $k_{l,x}=\\{1+\exp(-\tilde{x}^{T}\beta_{l})\\}^{-1}$ and
$k_{r,x}=\\{1+\exp(-\tilde{x}^{T}\beta_{r})\\}^{-1}$, for two _fixed_ shape
parameters $\gamma_{l}$ and $\gamma_{r}$. $h_{l,x}$ and $h_{r,x}$ are
respectively left-leaning (right-skewed) and right-leaning (left-skewed)
functions. We require that $\gamma_{l}>2$ and $\gamma_{r}>2$ to ensure that
both are strictly monotone and _convex_ , and thus provide a reasonable
approximation to the underlying true density in practice; see Lemma C.1 in
Appendix C for a precise result. The exact choices for
$\\{\gamma_{l},\gamma_{r}\\}$ will be further discussed in Section 3.5. The
working model may be generalized to capture non-linearity in $x$ using, say,
spline functions.

Beta mixtures have long been identified as a flexible modeling tool for
variables taking values in the unit interval; see Pounds and Morris (2003), Ji
et al. (2005), Parker and Rothenberg (1988), Markitsis and Lai (2010),
Migliorati et al. (2018), Ferrari and Cribari-Neto (2004) for related works.
In the context of covariate-adaptive multiple testing, Lei and Fithian (2018)
and Zhang and Chen (2020) employ a two-component beta-mixture model for the
$p$-values that consists of a uniform and another left-leaning beta component.
Our working model defined on the $u$-value scale can be viewed as a natural
extension of these works to capture important patterns in the $u$-value
distribution associated with two-sided covariate-adaptive testing.

The assessor can be constructed as the $\text{CLfdr}_{x}(\cdot)$ function with
respect to our working model (3.6). Since $h_{0}\equiv 1$, it follows that

(3.7)
$a_{x}(u)\equiv\frac{1-\pi_{l,x}-\pi_{r,x}}{1-{\pi}_{l,x}-{\pi}_{r,x}+{\pi}_{l,x}{h}_{l,x}(u)+{\pi}_{r,x}{h}_{r,x}(u)},\quad
0<u<1.$

The corresponding data-driven assessor is denoted by $\hat{a}_{x}(u)$ if the
parameters $\\{\theta_{l},\theta_{r},\beta_{l},\beta_{r}\\}$ are estimated
from the data for its construction.

###  3.3. Asymptotic ZAP

We now develop a direct data-driven version of the prototype algorithm in
Section 3.1. To construct $\hat{a}_{x}(u)$, we first obtain the maximum
likelihood estimates (MLE) of the unknown regression parameters
$\\{\theta_{l},\theta_{r},\beta_{l},\beta_{r}\\}$ with the data
$\\{U_{i},X_{i}\\}_{i=1}^{m}$; the EM algorithm for their computations are
provided in Appendix F.1. Denote $\hat{T}_{i}\equiv\hat{a}_{X_{i}}(U_{i})$,
and let $\hat{c}_{i}(\cdot)$ be its null distribution by treating
$\hat{a}_{X_{i}}(\cdot)$ as if it is pre-determined. With
$\hat{S}_{i}\equiv\hat{c}_{i}(\hat{T}_{i})$, the estimated mirror statistics
are correspondingly defined as
$\hat{T}^{\mathfrak{m}}_{i}\equiv\hat{c}_{i}^{-1}(1-\hat{S}_{i})$, which can
be computed numerically by performing quantile estimation. The FDP for a
candidate threshold $t$ can be estimated as

(3.8)
$\widehat{\mbox{FDP}}_{asymp}(t)\equiv\frac{1+\\#\\{i:\hat{T}^{\mathfrak{m}}_{i}\leq
t\\}}{\\#\\{i:\hat{T}_{i}\leq t\\}\vee 1}.$

Define $\hat{t}_{asymp}(\alpha)\equiv\sup\\{0\leq t\leq
1:\widehat{\mbox{FDP}}_{asymp}(t)\leq\alpha\\}$, and reject $i$ whenever
$\hat{T}_{i}\leq\hat{t}_{asymp}(\alpha)$. In practice, it suffices to consider
only the values of $\hat{T}_{1},\dots,\hat{T}_{m}$ as candidate thresholds.
This procedure is summarized in Algorithm 1.

1

2Construct $\hat{a}_{X_{i}}(\cdot)$’s using the MLEs obtained via the EM
algorithm in Appendix F.1 and compute $\hat{T}_{i}=\hat{a}_{X_{i}}(U_{i})$ for
each $i$.

3Compute the mirror statistics $\\{\hat{T}^{\mathfrak{m}}_{i}\\}_{i=1}^{m}$:
(i) Generate i.i.d. realizations $u_{1},\dots,u_{N}$ from $\mbox{Unif}(0,1)$
for a large $N$.  (ii) For each $i$, evaluate
$\hat{a}_{X_{i}}(u_{1}),\dots,\hat{a}_{X_{i}}(u_{N})$ to simulate the null
distribution $\hat{c}_{i}(\cdot)$. Compute $\hat{T}^{\mathfrak{m}}_{i}$ via
e.g. quantile() in R.

4Order $\\{\hat{T}_{i}\\}_{i=1}^{m}$ as
$\hat{T}_{(1)}\leq\dots\leq\hat{T}_{(m)}$. Reject $i$ if
$T_{i}\leq\hat{T}_{(k)}$, where
$k=\max\left\\{l\in\\{1,\dots,m\\}:\frac{1+\\#\\{i:\hat{T}^{\mathfrak{m}}_{i}\leq\hat{T}_{(l)}\\}}{l\vee
1}\leq\alpha\right\\}.$

Algorithm 1 Asymptotic ZAP

The main theory requires the following classical assumption from the
literature on misspecified models (White, 1981, 1982):

######  Assumption 1 (Existence of a unique maximizer).

The expected log-likelihood

$\mathbb{E}\log[(1-\pi_{l,X_{i}}-\pi_{r,X_{i}})+\pi_{l,X_{i}}h_{l,X_{i}}(U_{i})+\pi_{r,X_{i}}h_{r,X_{i}}(U_{i})]$

of the beta-mixture model (3.6) has a unique maximum at
$\\{\theta^{*},\beta^{*}\\}$ over $\theta\in\Theta$ and
$\beta\in\boldsymbol{B}$ for compact spaces $\Theta$ and $\boldsymbol{B}$,
where $\theta\equiv(\theta_{l}^{T},\theta_{r}^{T})^{T}$ and
$\beta\equiv(\beta_{l}^{T},\beta_{r}^{T})^{T}$. The expectation is taken with
respect to the true joint distribution of $\\{H_{i},Z_{i},X_{i}\\}$.

Together with Assumptions 2 \- 3 in Appendix D.1, which are standard
regularity and strong-law conditions, we can prove the following asymptotic
FDR controlling property.

######  Theorem 3.1.

Let $\hat{a}_{x}(\cdot)$ be constructed with the MLE
$\\{\hat{\theta},\hat{\beta}\\}=\underset{\theta,\beta}{\mathrm{argmax}}\sum_{i=1}^{m}\log
h_{X_{i}}(U_{i})$ of the beta-mixture model (3.6). Under Assumptions 1-3, the
procedure that rejects $i$ whenever $\hat{T}_{i}\leq\hat{t}_{asymp}(\alpha)$
controls the FDR asymptotically in the sense that
$\limsup_{m\rightarrow\infty}\mbox{FDR}\leq\alpha.$

We highlight two aspects of this result. First, it doesn’t require the
estimated assessor function to be a good proxy for $\text{CLfdr}_{x}(\cdot)$.
Hence, its theory is more attractive than that of Cai et al. (2019), which
requires consistent CLfdr estimates to ensure asymptotic FDR control. Second,
to establish Glivenko-Cantelli results (Lemma D.5) for the following three
empirical processes

$m^{-1}\sum_{i=1}^{m}I\left(\hat{T}_{i}\leq t\right),\quad
m^{-1}\sum_{i=1}^{m}(1-H_{i})I\left(\hat{T}_{i}\leq t\right),\quad
m^{-1}\sum_{i=1}^{m}I\left\\{\hat{S}_{i}\geq 1-\hat{c}_{i}(t)\right\\},$

typical of similar asymptotic analyses (Zhang and Chen, 2020, Storey et al.,
2004), we heavily utilize the convexity properties (Lemma C.2) of the
functional form in (3.7) to uniformly control the deviations of the estimated
assessors $\hat{a}_{X_{i}}(\cdot)$ from the assessors $a_{X_{i}}^{*}(\cdot)$
constructed with the population parameters $\\{\theta^{*},\beta^{*}\\}$; the
delicate techniques involved may be of independent interest.

###  3.4. Finite-sample ZAP

This section introduces an alternative ZAP procedure that offers finite-sample
control of the FDR. The operation again involves approximating the CLfdr
statistics via an assessor function. However, the thresholding step is based
on a more nuanced approach to FDP estimation inspired by the p-value method
AdaPT (Lei and Fithian, 2018). In this approach, multiple testing is conducted
in an iterative manner, where data are initially partially masked and then
gradually revealed at steps $t=0,1,\ldots$, with the thresholds sequentially
updated based on the revealed data at each step. In what follows, if
$g_{1}(\cdot)$ and $g_{2}(\cdot)$ are two functions defined on the same space,
$g_{1}\preceq g_{2}$ means $g_{1}(x)\leq g_{2}(x)$ for all $x$ in that space.
If $C$ is a constant, $g_{1}\preceq C$ means $g_{1}(x)\leq C$ for all $x$.
Similarly we can define $g_{1}\succeq g_{2}$ and $g_{1}\succeq C$.

Since the inherent convexity777Refer to Lemma C.2 of the assessor functional
form in (3.7) suggests that rejecting hypotheses with small values of
$T_{i}=a_{X_{i}}(U_{i})$ amounts to rejecting extreme $u$-values near 0 or 1,
our iterative algorithm emulates this essential operational characteristic of
the prototype procedure. We first divide the covariate values into a left and
right group based on the observed $u$-values:

$\mathcal{X}_{l}=\\{X_{i}:U_{i}\leq 0.5\\}\text{ and
}\mathcal{X}_{r}=\\{X_{i}:U_{i}>0.5\\}.$

At each step $t=0,1,\dots$, let $s_{l,t}:\mathcal{X}_{l}\rightarrow[0,0.25]$
and $s_{r,t}:\mathcal{X}_{r}\rightarrow[0.75,1]$ denote two corresponding
thresholding functions, and define the candidate rejection set
$\mathcal{R}_{t}\equiv\mathcal{R}_{l,t}\cup\mathcal{R}_{r,t}$, where

(3.9) $\mathcal{R}_{l,t}\equiv\\{i:U_{i}\leq s_{l,t}(X_{i})\wedge 0.5\\}\text{
and }\mathcal{R}_{r,t}\equiv\\{i:U_{i}\geq s_{r,t}(X_{i})\vee 0.5\\}.$

Let $\mathcal{A}_{t}\equiv\mathcal{A}_{l,t}\cup\mathcal{A}_{r,t}$ be the
corresponding set of “accepted” hypotheses, where

$\mathcal{A}_{l,t}\equiv\\{i:0.5-s_{l,t}(X_{i})\leq U_{i}\leq 0.5\\}\text{ and
}\mathcal{A}_{r,t}\equiv\\{i:0.5<U_{i}\leq 1.5-s_{r,t}(X_{i})\\}.$

Intuitively, $|\mathcal{A}_{l,t}|$ estimates the number of false rejections in
the left candidate rejection set $\mathcal{R}_{l,t}$: Given $H_{i}=0$ and
$U_{i}\leq 0.5$, the events $\\{U_{i}<s_{l,t}(X_{i})\\}$ and
$\\{U_{i}>0.5-s_{l,t}(X_{i})\\}$ are equally likely. The FDP of possibly
rejecting $\mathcal{R}_{t}$ at step $t$ can then be estimated as

(3.10)
$\widehat{\mbox{FDP}}_{finite}(t)=\frac{1+|\mathcal{A}_{t}|}{|\mathcal{R}_{t}|\vee
1}.$

If $\widehat{\text{FDP}}_{finite}(t)\leq\alpha$, the algorithm terminates and
the hypotheses in $\mathcal{R}_{t}$ are rejected. Otherwise, the algorithm
proceeds to the next step $t+1$ and updates the two thresholding functions
under two restrictions. First, it must be that $s_{l,t+1}\preceq s_{l,t}$ and
$s_{r,t+1}\succeq s_{r,t}$; this ensures that $\mathcal{R}_{t}$ shrinks in
size as $t$ increases. Second, $s_{l,t+1}$ and $s_{r,t+1}$ must be updated
based on the knowledge of $|\mathcal{R}_{t}|$, $|\mathcal{A}_{t}|$ and the
partially masked data $\\{\tilde{U}_{t,i},X_{i}\\}_{i=1}^{m}$ only, where

(3.11) $\widetilde{U}_{t,i}\equiv\begin{dcases*}U_{i}&if
$U_{i}\not\in\mathcal{A}_{t}\cup\mathcal{R}_{t}$\\\
\\{U_{i},\widecheck{U}_{i}\\}&if
$U_{i}\in\mathcal{A}_{t}\cup\mathcal{R}_{t}$\end{dcases*}$

is a singleton or a two-element set depending on whether $i$ is in the
“masked” set $\mathcal{A}_{t}\cup\mathcal{R}_{t}$, and $\widecheck{U}_{i}$ is
the “reflection” of $U_{i}$ about the “middle” axis at $u=0.25$ or $u=0.75$,
depending on which group (left or right) $U_{i}$ belongs to:

$\widecheck{U}_{i}\equiv(1.5-U_{i})I(U_{i}>0.5)+(0.5-U_{i})I(U_{i}\leq 0.5).$

For example, if the underlying $U_{i}$ is $0.1$ and $i$ is masked at step $t$,
the algorithm can only update for $s_{l,t+1}$and $s_{r,t+1}$ with the partial
knowledge that $U_{i}$ is either $0.1$ or its reflection value $0.4$.
Algorithm 4 in Appendix F.2 describes one such updating scheme which first
applies an EM algorithm (Appendix F.3) acting only on the partially masked
data to estimate the beta-mixture model (3.6). With the latter subsequently
used to produce estimated $\text{CLfdr}_{i}$’s for comparing the significance
of the hypotheses in the masked set $\mathcal{A}_{t}\cup\mathcal{R}_{t}$, the
thresholding functions are updated in such a way that the masked $i$ deemed to
be the least significant will have its $u$-value revealed at step $t+1$.
Figure 3.1 illustrates how the data $\\{U_{i},X_{i}\\}_{i=1}^{m}$ are
partitioned into $\mathcal{A}_{t}$, $\mathcal{R}_{t}$ and the unmasked set
$\\{1,\dots,m\\}\backslash\\{\mathcal{A}_{t}\cup\mathcal{R}_{t}\\}$ at a given
step $t$, based on Example 2.2 in Section 2.2. In particular, we remark that
the algorithm cannot tell the true data point from a given red-pink (blue-
cyan) pair in the plot $(b)$ where the reflection points
$\\{\widecheck{U}_{i}\\}_{i\in\mathcal{A}_{t}\cup\mathcal{R}_{t}}$ are also
shown.

Figure 3.1.  Illustration of Algorithm 2 at a step $t$ based on Example 2.2,
$m=2000$. (a): The red, blue and grey are points in the respective sets
$\mathcal{R}_{t}$, $\mathcal{A}_{t}$ and
$\\{1,\dots,m\\}\backslash\\{\mathcal{A}_{t}\cup\mathcal{R}_{t}\\}$, where
$\widehat{\text{FDP}}_{finite}(t)=193/389\approx 0.5$. (b): The reflections of
the points in $\mathcal{R}_{t}$ and $\mathcal{A}_{t}$ are respectively shown
in pink and cyan. A pink $\widecheck{U}_{i}$ below (above) $0.5$ is the
reflection of a red $U_{i}$ with the same covariate value about the middle
axis at $u=0.25$ ($u=0.75$) of the left (right) group; the cyan are the
reflections of the blue.

The steps described above are summarized in Algorithm 2, whose finite-sample
FDR controlling property is stated in Theorem 3.2.

Data: $\\{U_{i},X_{i}\\}_{i=1}^{m}$

Input: FDR level $\alpha$, initial thresholding functions $s_{l,0}\preceq 0.5$
and $s_{r,0}\succeq 0.5$;

1

2for  _t= 0,1 …,_ do

3 Compute $\widehat{\mbox{FDP}}_{finite}(t)$ in (3.10);

4 if  _$\widehat{\mbox{FDP}}_{finite}(t) >\alpha$_ then

5 Update $s_{l,t+1}$ and $s_{r,t+1}$ while respecting the two conditions in
Theorem 3.2. E.g. Apply Algorithm 4 in Appendix F.2;

6

7 else

8 Record $\mathcal{R}_{t}$; break;

9

10 end if

11

12 end for

Output: Reject all hypotheses in $\mathcal{R}_{t}$.

Algorithm 2 Finite-sample ZAP

######  Theorem 3.2 (Finite-sample FDR control).

Under the conditions that  (i) $s_{r,t+1}\succeq s_{r,t}$and $s_{l,t+1}\preceq
s_{l,t}$ and  (ii) $s_{r,t+1}$, and $s_{l,t+1}$ are updated based on
$|\mathcal{R}_{t}|$, $|\mathcal{A}_{t}|$ and
$\\{\tilde{U}_{t,i},X_{i}\\}_{i=1}^{m}$ only,  Algorithm 2 controls the FDR
under $\alpha$ for finite samples. Specifically, we have

$\mathbb{E}\left[FDP\Big{|}\\{H_{i},X_{i}\\}_{i=1}^{m}\right]\leq\alpha.$

Lastly, we highlight a crucial difference between Algorithm 2 and AdaPT in the
present context. Operating on the two-sided $p$-values, AdaPT proceeds
iteratively with a single thresholding function $s_{t}$ defined on
$\\{X_{i}\\}_{i=1}^{m}$ such that $s_{t}\preceq 0.5$, and the ratio
$\frac{1+|\\{i:P_{i}\geq 1-s_{t}(X_{i})\\}|}{1\vee|\\{i:P_{i}\leq
s_{t}(X_{i})\\}|}$ is used as an FDP estimator for the candidate rejection set
$\\{i:P_{i}\leq s_{t}(X_{i})\\}$. It is easy to see that

(3.12) $P_{i}\leq s_{t}(X_{i})\iff U_{i}\leq s_{t}(X_{i})/2\text{ or
}U_{i}\geq 1-s_{t}(X_{i})/2.$

Hence, on the $u$-value scale, AdaPT always adopts symmetric rejection regions
about $u=0.5$. By contrast, Algorithm 2 employs two different thresholding
functions $s_{l,t}$ and $s_{r,t}$, which allow for asymmetric rejection
regions, and therefore provides additional flexibility to fully capitialize on
covariate information for two-sided tests. As seen in Figure 3.1, the pattern
of the candidate rejection points in red agrees with the middle panel of
Figure 2.1; as the covariate increases from $-1$ to $1$, the algorithm’s
rejection priorities change from the $u$-values near $0$ to those near $1$.

###  3.5. Implementation details

An R package zap for our two data-driven methods is available on
https://github.com/dmhleung/zap, and we shall discuss further details of their
implementation.

For both data-driven procedures, the shape parameters
$\\{\gamma_{l},\gamma_{r}\\}$ of the working model need to be pre-specified
before running the EM algorithms. While requiring $\gamma_{l},\gamma_{r}>2$
ensures a convex shape for the three-component beta-mixture density (Lemma
C.1), we recommend choosing $(\gamma_{l},\gamma_{r})=(4,4)$ as a default,
which has yielded consistently good performance in our numerical studies.

To illustrate the effectiveness of our recommendation, we simulate 8000 i.i.d.
$z$-values $Z_{1},\dots,Z_{8000}$ from the normal mixture model

(3.13) $0.78f_{0}(z)+0.15\phi(z+1.5)+0.07\phi(z-2)$

without any covariates. The histogram of the corresponding $u$-values is
plotted in Figure 3.2(a), overlaid with the true underlying density function,
as well as estimated densities of the beta mixture (3.6) fitted with
regression intercepts only, where $(\gamma_{l},\gamma_{r})$ is respectively
fixed at $(1,1)$ and $(4,4)$. When modeling p-values with a two-component beta
mixture, Lei and Fithian (2018) and Zhang and Chen (2020) set an analogous
shape parameter to be $1$, so $(\gamma_{l},\gamma_{r})=(1,1)$ would be a
seemingly natural choice to extend their model for two-sided tests. Both
fitted densities visually coincide with the true density, attesting to the
flexibility of beta mixtures for modeling data on the unit interval. However,
the estimated component probabilities differ significantly for
$(\gamma_{l},\gamma_{r})=(1,1)$ vs $(\gamma_{l},\gamma_{r})=(4,4)$. In Figure
3.2(b), we present the estimated quantities pertaining to the non-null
components. We can see that setting $(\gamma_{l},\gamma_{r})=(1,1)$ has
drastically overestimated the left and right non-null probabilities, whereas
setting $(\gamma_{l},\gamma_{r})=(4,4)$ provides good approximations to the
truths.

Figure 3.2.  (a): Histogram of $U_{i}=\Phi(Z_{i})$ generated by (3.13). The
red curve is the true density; while the blue and green curves respectively
correspond to the estimated densities of our beta-mixture model with
$(\gamma_{l},\gamma_{r})$ set at $(1,1)$ and $(4,4)$. (b): The long dashed
blue curve is the entire non-null component of the solid blue estimated
density in $(a)$. The two dotted blue curves are the left-leaning and right-
leaning non-null components that add up to the long dashed blue. The green
curves are constructed analogously with respect to the solid green density in
$(a)$. The legend shows the estimated probabilities for the left and right-
leaning components, where the subscript “$x$” is omitted from $\hat{\pi}_{l}$
and $\hat{\pi}_{r}$ as the fit uses intercepts only. (c): Plot of the left-
leaning beta density $B(k_{l},\gamma_{l})^{-1}u^{k_{l}-1}(1-u)^{\gamma_{l}-1}$
for $k_{l}=0.1$ and different values of $\gamma_{l}$.

To gain insight into why larger shape parameters are preferred, in Figure
3.2(c) we plot the density of a left-leaning beta density

$B(k_{l},\gamma_{l})^{-1}u^{k_{l}-1}(1-u)^{\gamma_{l}-1}$

for different values of $\gamma_{l}$ and a fixed $k_{l}=0.1$, which supposedly
captures the negative effects in two-sided tests. We can see that small values
of $\gamma_{l}$ tend to yield a density component that slants in the middle of
the unit interval. As a result, when added to another right-leaning beta
density for the positive effects with a similar but mirroring shape, it gives
rise to an overall non-null density component with a large plateau in the
middle of the interval $(0,1)$ akin to the U-shaped blue curve in Figure
3.2$(b)$. This inflates the non-null component probability estimates. In
contrast, larger values of $\gamma_{l}$ and $\gamma_{r}$ effectively mitigate
the issue by rendering sharply convex non-null component densities like the
purple and pink curves in Figure 3.2(c), avoiding overestimation of the non-
null probabilities and leading to a better estimated overall non-null density
component like the U-shaped green curve in Figure 3.2$(b)$. More setups are
experimented in Appendix G; the careful choice of
$(\gamma_{l},\gamma_{r})=(4,4)$ produces reasonable probability estimates
throughout.

Other aspects of implementation are as follows. For the asymptotic method
(Algorithm 1), since a large $N$ allows us to compute the mirror statistics up
to arbitrary precision, we evaluate at $N=50000$ uniform realizations by
default. For the iterative finite-sample method (Algorithm 2), we set the
initial thresholding functions as $s_{l0}\equiv 0.2$ and $s_{r0}\equiv 0.8$,
but other values close to $0.25$ and $0.75$ tend to be equally effective. We
also update the thresholding functions every $\lceil m/100\rceil$ steps.
Ideally one would want to update at every step along the way to reveal the
masked $u$-values sooner. However, it is more practical to carry out
intermittent updates since the EM component involved in Algorithm 4 is
computationally costly. Lastly, one can also perform feature selection at any
step if $X_{i}$ is multivariate, as long as it is done properly based on the
masked data, akin to what was suggested by Lei and Fithian (2018, Section
4.2). We have not performed this step for simplicity.

##  4\. Numerical studies

We conduct numerical studies to gauge the performance of ZAP alongside other
methods on both simulated and real data. For expositional considerations, here
we only limit the comparisons to a selection of representative FDR methods.
This makes the ensuing graphs (Figures 4.1-4.2) less cluttered with lines and
easier to read. Comparisons with more methods in the literature are included
in Appendices H and I, but the basic conclusions do not change. Here, we
consider:

  1. (a)

ZAP (asymp): Algorithm 1 with specifications described in Section 3.5.

  2. (b)

ZAP (finite): Algorithm 2 with specifications described in Section 3.5.

  3. (c)

CAMT: the covariate-adaptive multiple testing method (Zhang and Chen, 2020).

  4. (d)

AdaPT: the adaptive $p$-value thresholding method (Lei and Fithian, 2018).
Their working model is updated based on the EM algorithm for every $\lceil
m/100\rceil$ steps; other default specifications are chosen based on the R
package adaptMT.

  5. (e)

IHW: Independent hypothesis weighting (Ignatiadis et al., 2016, Ignatiadis and
Huber, 2021), based on the default provided by the R package IHW, version
1.16.0, which provides asymptotic FDR control (Ignatiadis and Huber, 2021,
Proposition 1). This method only handles univariate covariates.

  6. (f)

FDRreg: false discovery rate regression method (Scott et al., 2015). The
theoretical null $\mathcal{N}(0,1)$ has been used.

  7. (g)

AdaPT-$\text{GMM}_{g}$: A $z$-value based variant of AdaPT by Chao and Fithian
(2021) which is also a data masking procedure similar to ZAP (finite); unlike
ZAP (finite), it instead employs a covariate-dependent Gaussian mixture
working model for the $z$-values, which is updated with an EM algorithm for
every $\lceil m/100\rceil$ steps. Other default specifications are chosen as
in the R package AdaPTGMM.

Among them, ZAP, AdaPT-$\text{GMM}_{g}$, and FDRreg are $z$-valued based,
while all other methods are $p$-value based.

###  4.1. Simulated data

We simulate data to test $m=5000$ hypotheses. Two-dimensional covariates
$X_{i}=(X_{1i},X_{2i})^{T}$, $i=1,\ldots,m$, are independently generated from
the bivariate normal distribution $\mathcal{N}\left\\{(\begin{smallmatrix}0\\\
0\end{smallmatrix}),(\begin{smallmatrix}1/2&0\\\
0&1/2\end{smallmatrix})\right\\}.$ Conditional on
$X_{i}=x\equiv(x_{1},x_{2})^{T}$, $Z_{i}$ is generated with a normal mixture
density

(4.1)
$(1-w_{l,x}-w_{r,x})f_{0}(z)+w_{l,x}\phi\left(\frac{z-\mu_{l,x}}{\sigma}\right)+w_{r,x}\phi\left(\frac{z-\mu_{r,x}}{\sigma}\right),$

where $w_{l,x}$ and $w_{r,x}$ 888These data generating probabilities
$\\{w_{l,x},w_{r,x}\\}$ should again be distinguished from
$\\{\pi_{l,x},\pi_{r,x}\\}$ in the working model (3.6). are probabilities that
control the sparsity levels of negative and positive effects, $\mu_{l,x}<0$
and $\mu_{r,x}>0$ are negative and positive non-null normal means, and
$\sigma^{2}$ is the variance of the alternative components. The covariate-
adjusted overall non-null density is then given by

(4.2)
$f_{1,x}(z)=\frac{w_{l,x}\phi(\frac{z-\mu_{l,x}}{\sigma})+w_{r,x}\phi(\frac{z-\mu_{r,x}}{\sigma})}{w_{l,x}+w_{r,x}}.$

We shall allow $\\{w_{l,x},w_{r,x},\mu_{l,x},\mu_{r,x}\\}$ to depend on $x$ in
different ways to induce the simulation setups below, which can be considered
as more realistic versions of the stylized examples in Section 2.2. We fix
$\sigma^{2}=1$ in this section, while other values of $\sigma^{2}$ will be
explored in Appendix H, which also contains more extensive simulation studies
as well as a comparison of the computational efficiency of the different
methods. Note that the sum $X_{\bullet i}\equiv X_{1i}+X_{2i}$ of the
covariate components is $\mathcal{N}(0,1)$-distributed, and $x_{\bullet}\equiv
x_{1}+x_{2}$ will denote a realized value of it below.

  1. Setup 1 (Asymmetric alternatives). The quantities in (4.1) are

$w_{r,x}=\frac{1}{1+\exp(-\eta-\zeta x_{\bullet})},\ \
\mu_{r,x}=\frac{2\varepsilon}{1+\exp(-\zeta x_{\bullet})},\ \ w_{l,x}=0,\ \
\mu_{l,x}=0,$

with the simulation parameters ranging as

$\zeta\in\\{0,0.5,1\\},\quad\varepsilon\in\\{1.3,1.5,1.7,1.9,2.1\\}\text{ and
}\eta=-2.$

Since $w_{l,x}=0$, all the non-null statistics come from the right centered
alternative density $\phi(z-\mu_{r,x})$. We briefly explain the simulation
parameters. $\varepsilon$ is an effect size parameter. Generally, $\zeta$
controls the _informativeness_ of the covariates in relation to both the non-
null probabilities and alternative means: when $\zeta>0$, a greater value of
$x_{\bullet}$ makes the signals denser and stronger (i.e. $w_{r,x}$ and
$\mu_{r,x}$ become larger). The value of $\eta$ controls the sparsity levels.
For example, when the covariates are non-informative at $\zeta=0$, setting
$\eta=-2$ yields a baseline signal proportion of roughly $12\%$, i.e.
$w_{r,x}=w_{l,x}+w_{r,x}=11.9\%$. Note that $f_{1,x}$ in (4.2) varies in $x$
given the dependence of $\mu_{r,x}$ on $x$, so (3.3) is an invalid assumption.

  2. Setup 2 (Unbalanced covariate effects on the non-null proportions). Let

$w_{r,x}=\frac{\exp(\zeta x_{\bullet})}{\exp(-\eta)+\exp(-\zeta
x_{\bullet})+\exp(\zeta x_{\bullet})},\;w_{l,x}=\frac{\exp(-\zeta
x_{\bullet})}{\exp(-\eta)+\exp(-\zeta x_{\bullet})+\exp(\zeta x_{\bullet})},$

$\mu_{rx}=\varepsilon$ and $\mu_{lx}=-\varepsilon$. We fix $\eta=-2.5$ and
vary other parameters in the range

$\zeta\in\\{0,0.7,1\\},\quad\varepsilon\in\\{1.3,1.5,1.7,1.9,2.1\\}.$

Only $w_{l,x}$ and $w_{r,x}$ depend on the covariate value: for $\zeta>0$,
$w_{r,x}$ increases and $w_{l,x}$ decreases as $x_{\bullet}$ increases, and
vice versa as $x_{\bullet}$ decreases. In consideration of (4.2), the
conditional non-null density $f_{1,x}(z)$ will change sharply in shape from
concentrating on negative $z$-values to concentrating on positive $z$-values
as $x_{\bullet}$ increases from being negative to positive. This relationship
provides important structural information which can be leveraged for enhancing
the power. However, if one collapses the $z$-values into two-sided $p$-values,
then the analogous conditional $p$-value density is less likely to capture
drastic changes in $x_{\bullet}$, since both very negative and positive
$x_{\bullet}$ can correspond to very small $p$-values, making the interactive
relationship between the $p$-values and the covariates less pronounced.
Intuitively, this would lead to power loss of $p$-value based methods. The
choice of $\eta$ corresponds to a baseline signal proportion of roughly $14\%$
when $\zeta=0$.

  3. Setup 3 (Unbalanced covariate effects on the alternative means). Let

$w_{r,x}=\frac{1/2}{1+\exp(-\eta)},w_{l,x}=\frac{1/2}{1+\exp(-\eta)},\mu_{r,x}=\frac{2\varepsilon}{1+\exp(-\zeta
x_{\bullet})},\mu_{l,x}=\frac{-2\varepsilon}{1+\exp(\zeta x_{\bullet})}.$

The simulation parameters range as

$\zeta\in\\{0,1.5,3\\},\quad\varepsilon\in\\{1.3,1.5,1.7,1.9,2.1\\}\text{ and
}\eta=-2.$

Our choice of $\eta$ corresponds to the a baseline signal proportion of
roughly $12\%$ when $\zeta=0$. When the covariates are informative
($\zeta>0$), $\mu_{r,x}$ and $\mu_{l,x}$ respectively become more positive and
less negative as $x_{\bullet}$ increases. Such a directional relationship can
be potentially exploited by ZAP for power improvement. However, if one
collapses the $z$-values into $p$-values, then under $H_{i}=1$ both very
positive and negative values of $X_{\bullet i}$ can imply a small $P_{i}$, and
the interactive relationship between the main statistic $P_{i}$ and auxiliary
statistic $X_{i}$ will be much weakened.

Figure 4.1.  FDR and TPR performances of different methods under Setup 1 - 3 .
All methods are applied at a targeted FDR level of $0.05$ (the horizontal
black lines). The x-axes show the values of $\varepsilon$. non-informative,
more informative and most informative correspond to different values of
$\zeta$ from the smallest to the largest.

We apply the seven methods at the nominal FDR level 0.05. Since IHW can only
handle univariate covariates, it is applied with $X_{\bullet i}$, which is an
effective summary covariate in all three setups. The simulation results are
reported in Figure 4.1, where the empirical FDR and TPR levels of different
methods are computed based on 150 repetitions. The following observations can
be made:

  1. (a)

Asymptotic ZAP, depicted in blue, is in general more powerful than finite-
sample ZAP, depicted in red. This is likely attributable to the latter’s
information loss from the “$u$-value masking” step, which can be perceived as
a necessary trade-off for the strong finite-sample FDR controlling property.
Even so, the asymptotic ZAP has its FDR controlled under or around the nominal
level $0.05$.

  2. (b)

Both asymptotic and finite-sample ZAP methods demonstrate superior
performances over the $p$-value based methods (CAMT, AdaPT and IHW). The gains
in power become more substantial when the covariates become more informative.
The $z$-value based method AdaPT-$\text{GMM}_{g}$ has quite comparable power
to the finite-sample ZAP in these setups too, and we will provide more
discussion on the comparison between the them in Section 5.

  3. (c)

The covariate-adjusted non-null density (4.2) depends on $x$ for all three
setups, so FDRreg, which makes the conflicting assumption in (3.3), is
possibly invalid for FDR control. This is indeed observed for a number of
settings. Moreover, it can’t match the power of ZAP in Setup 2, likely because
the assumption $f_{1,x}(z)\equiv f_{1}(z)$ itself obstructs the interactive
information between the $z$-values and the covariates to be utilized.

  4. (d)

In Setup 3, ZAP only has slightly discernible power advantage over the other
methods when the covariates are informative. In fact, no current FDR methods
can demonstrate near-optimal power under this setup, as shown in Appendix H.
In Section 5, we discuss possible future work that might address this.

###  4.2. Real data

This section investigates the performance of ZAP using several publicly
available real datasets summarized in Table 1. Three datasets (bottomly,
airway, hippo) are generated by RNA sequencing (RNA-Seq) experiments for
detecting differential expressions in transcriptomes, where the primary
statistic $Z_{i}$ measures the observed difference in the expression level of
a gene under two experimental conditions. Meanwhile, an auxiliary covariate,
the average normalized read count for each gene, is collected alongside the
primary data. The datasets bottomly and airway have been analyzed by the works
of Ignatiadis et al. (2016), Lei and Fithian (2018), Zhang and Chen (2020)
with the methods IHW, AdaPT and CAMT respectively. The more recent data set
hippo (Harris et al., 2019) is generated by the cutting-edge single-cell RNA
(sc-RNA) sequencing technology to study differential expressions in mouse
hippocampus. For all datasets above, we have adopted the standard data pre-
processing step, which filters out genes with excessively low read counts
across samples before further downstream analyses (Chen et al., 2016) such as
model fitting and multiple testing. This is a common practice among
bioinformaticians for a number of reasons; see Appendix I.1 for more
discussion. The fourth dataset is based on the experiments in Smith and Kohn
(2008) and Kelly et al. (2010), where each $Z_{i}$ is a normalized test
statistic that, for a given pair of neurons in the primary visual cortex,
measures how synchronous their spike trains are, and Scott et al. (2015) has
applied FDRreg to it for detecting neural interactions. Correspondingly, each
such hypothesis has two covariates: the distance and the correlation of the
“tuning curves” between the two activated neurons. We have named this dataset
scott for short.

Name | # tests |  Brief description
---|---|---
bottomly | 11484 |  DE in striatum for the two mouse strains C57BL/6J(B6) and DBA/2J(D2); bulk RNA-seq (Bottomly et al., 2011).
airway | 20941 |  DE in human airway smooth muscle cell lines in response to dexamethasone; bulk RNA-seq (Himes et al., 2014).
hippo | 15000 |  DE in mouse hippocampus in response to enzymatic dissociation in comparison to standard tissue homogenization; scRNA-seq (Harris et al., 2019).
scott | 7004 |  Synchronous firing of pairs of neurons, based on neuron recordings in the primary visual cortex of an anesthetized monkey in response to visual stimuli (Scott et al., 2015).
Table 1. Description of four real datasets. “# tests” shows the number of
tests for each dataset after any necessary data pre-processing. DE =
Differential Expression.  Figure 4.2. $\bf(a)-(d)$ plot the numbers of
rejections for different methods across datasets, against targeted FDR level
at $0.01,0.05,0.1,0.15,0.2$. $\bf(e)$ and $\bf(f)$ are respectively the
histograms of the “$u$-values” for the hippo and scott datasets.

For the RNA-seq datasets, all the methods that accommodate multivariate
covariates (CAMT, AdaPT, FDRreg, AdaPT-$\text{GMM}_{g}$ and the two methods of
ZAP) are applied with the log mean normalized read count expanded by a natural
cubic spline basis with 4 interior knots, using the ns function in the R
package splines (with its df argument set to $6$), and IHW is applied with the
original log mean normalized read count as it can only handle a univariate
covariate. As there are two covariates for the neural dataset scott, IHW is
not applied, and following what was done in Scott et al. (2015), the
multivariate methods are applied with each of the two covariates expanded by a
B-spline basis using the bs function in R with its argument df set to 3, which
results in six expanded covariates in total.

The number of rejections for the various methods are shown in Figure
4.2$(a)$-$(d)$, and ZAP has attained top power performances in general. For
the dataset bottomly, FDRreg shows moderately more rejections than other
methods, but the power gain may be due to overflow in FDR resulting from
making the assumption in (3.3). For the dataset airway, it is surprising to
see the power of AdaPT-$\text{GMM}_{g}$ drops quite a bit beyond the target
FDR level $0.1$. The histogram plots in Figure I.2 of Appendix I show that the
$u$-values are almost symmetrically distributed for these two datasets, which
suggests that $p$-value and $z$-value based methods tend to have comparable
power, unless the reduction to $p$-values fails to capture the interactive
information between the $z$-values and covariates. For the datasets hippo and
scott, the histograms, which are shown in Figure 4.2$(e)$ and $(f)$, show that
the $u$-value distribution is asymmetric. This can explain why the $z$-value
based methods exhibits conspicuous power improvement over the $p$-value
methods in Figure 4.2 $(c)$ and $(d)$, which are in agreement with the
intuition that the former are more capable of exploiting the distributional
asymmetry. Similar to what we observed in the simulation studies, the
asymptotic ZAP tends to reject more hypotheses than the finite-sample ZAP,
whose FDR controlling validity is based on fewer assumptions (Theorem 3.2).
For the scott dataset, the increase in power of the finite-sample ZAP trails
off a bit beyond the target FDR level of $0.15$. Overall, the real data
analyses affirm that $z$-value based approaches to covariate-adaptive testing
promise to better exploit the full data $\\{Z_{i},X_{i}\\}_{i=1}^{m}$ to boost
testing power. In Section 5, we will offer a general comment on the power of
data masking procedures.

##  5\. Discussion

We have introduced ZAP, which is a $z$-value based covariate-adaptive testing
framework that offers control of the FDR under minimal assumptions. The main
thrust of our proposal is to avoid the common data reduction step of forming
two-sided $p$-values used by most other covariate-adaptive methods in the
recent literature, so as to preserve the intact structural information in the
data as a starting point to devise more powerful procedures. While there is no
“one-size-fits-all” solution to all FDR analysis problems, as is seen in the
extensive simulation studies in the recent paper of Korthauer et al. (2019),
we believe the current form of ZAP is a competitive choice in many covariate-
adaptive testing situations.

Although the finite-sample ZAP offers exact FDR control (Theorem 3.2), the
requisite data masking machinery incurs some information loss as a trade-off.
This often injects extra instability into its power performance; for one
thing, a working model estimated based on the partially masked data does not
necessarily fit the original full data as well, particularly when it is
flexible enough to overfit the partially masked data. On the other hand, a
reasonably flexible working model is needed to better capture the true
underlying data generating mechanism, which will be translated into power.
Prior works on data masking have also mentioned similar power-instability
issues; see Lei and Fithian (2018, Appendix A.3) and Lei et al. (2017, Section
3.2) for instance. Hence, when the number of hypotheses $m$ is large enough
(e.g. $m\geq 5000$) to justify the large-sample theory, we generally advocate
the use of the asymptotic version of ZAP.

We have chosen to operate on the $u$-value scale, leveraging a covariate-
dependent beta mixture as our working model. A natural alternative choice is
to employ some form of Gaussian mixture models defined on the original
$z$-value scale, which is also popular among researchers for performing FDR
testing in different applications (McLachlan et al., 2006, Nguyen et al.,
2018). Regardless of the model choice, since the asymptotic ZAP aims to
directly mimic the optimal oracle procedure, for it to be powerful the
assessor function $a_{x}(\cdot)$ should be a reasonably accurate proxy for the
true conditional local false discovery rate function
$\text{CLfdr}_{x}(\cdot)$. We have had success in this regard with our
carefully calibrated beta-mixture model; as explained in Section 3.5, setting
$(\gamma_{l},\gamma_{r})=(4,4)$ generally renders reasonably close estimates
for the (non-)null probabilities, which are critical for constructing a more
accurate assessor function with the formula in (3.7) with respect to the
working model. In unreported simulations, we have also experimented with the
normal mixture model, but with discouraging results; namely, it is very hard
to consistently obtain good estimates of the mixing probabilities, even with
the full data. This is not a surprise since Gaussian deconvolution is known to
be a very difficult problem; see Fan (1991) for a theoretical discussion, or
Scott et al. (2015, Section 2.1) for an illustration with a simple simulation
setup.

Chao and Fithian (2021) appeared at around the same time as our present paper.
(We thank an anonymous referee for pointing out this work.) Being most
comparable to our finite-sample ZAP method, their AdaPT-$\text{GMM}_{g}$ is a
$z$-value based data masking algorithm that does leverage a working Gaussian
mixture density with $K$ components of the form

(5.1)
$Z_{i}|X_{i}=x\sim\sum_{k=1}^{K}\pi_{k,x}\phi\left(\frac{z-\mu_{k}}{\sigma_{k}}\right)$

to iteratively reveal the masked $z$-values, where each $\pi_{k,x}$ is a
component probability depending on $x$ based on a model of choice (e.g. a
neural network), and $\mu_{k}$ and $\sigma_{k}$ are component mean and
standard deviation to be estimated. Unlike our current version of the finite-
sample ZAP which uses estimated $\text{CLfdr}_{i}$’s as the measuring scores
to decide which one of the masked $i$’s is the least promising and to be
revealed in the next step (Algorithm 4 in the Appendix F.2),
AdaPT-$\text{GMM}_{g}$, in our style of notation, uses conditional
probabilities of the type

(5.2) $P\Bigl{(}|\widecheck{Z}_{i}|>|Z_{i}|\;\ \Big{|}\;\
X_{i},\\{Z_{i},\widecheck{Z}_{i}\\}\Bigr{)}$

as scores to compare the masked $i$’s, where $\\{Z_{i},\widecheck{Z}_{i}\\}$
is a masked two-element set containing the value of $Z_{i}$ and its reflection
that is defined similarly to $\\{U_{i},\widecheck{U}_{i}\\}$ in (3.11) but on
the $z$-value scale. Intuitively, a masked $i$ with the largest such
conditional probability is considered the least promising to contain a signal
and should be revealed in the next step. Here, a Gaussian mixture model is
still a sensible working model to use because the quantity in (5.2) can be
reasonably estimated without relying on accurate deconvolution; see Chao and
Fithian (2021) for more details. In this vein, it will be interesting to
explore a mixture-of-experts working model (Chamroukhi and Huynh, 2019, Nguyen
et al., 2016), which is essentially (5.1) but with the means $\mu_{k}$ also
depending on covariates, particularly for a situation like Setup 3 in Section
4.1 where the covariate is informative via influencing the alternative means.
We leave this to future research that may be opportune in other instances.

## References

  * Arias-Castro et al. (2017) Arias-Castro, E., Chen, S., et al. (2017).  “Distribution-free multiple testing.”  Electronic Journal of Statistics, 11(1): 1983–2001.
  * Barber and Candès (2015) Barber, R. F. and Candès, E. J. (2015).  “Controlling the false discovery rate via knockoffs.”  The Annals of Statistics, 43(5): 2055–2085.
  * Barber and Candès (2019) — (2019).  “A knockoff filter for high-dimensional selective inference.”  The Annals of Statistics, 47(5): 2504–2537.
  * Basu et al. (2018) Basu, P., Cai, T. T., Das, K., and Sun, W. (2018).  “Weighted false discovery rate control in large-scale multiple testing.”  Journal of the American Statistical Association, 113(523): 1172–1183.
  * Benjamini and Hochberg (1995) Benjamini, Y. and Hochberg, Y. (1995).  “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.
  * Benjamini and Hochberg (2000) — (2000).  “On the adaptive control of the false discovery rate in multiple testing with independent statistics.”  Journal of educational and Behavioral Statistics, 25(1): 60–83.
  * Boca and Leek (2018) Boca, S. M. and Leek, J. T. (2018).  “A direct approach to estimating false discovery rates conditional on covariates.”  PeerJ, 6: e6035.
  * Bottomly et al. (2011) Bottomly, D., Walter, N. A., Hunter, J. E., Darakjian, P., Kawane, S., Buck, K. J., Searles, R. P., Mooney, M., McWeeney, S. K., and Hitzemann, R. (2011).  “Evaluating gene expression in C57BL/6J and DBA/2J mouse striatum using RNA-Seq and microarrays.”  PloS one, 6(3): e17820.
  * Cai and Sun (2009) Cai, T. T. and Sun, W. (2009).  “Simultaneous Testing of Grouped Hypotheses: Finding Needles in Multiple Haystacks.”  J. Amer. Statist. Assoc., 104: 1467–1481.
  * Cai et al. (2019) Cai, T. T., Sun, W., and Wang, W. (2019).  “Covariate-assisted ranking and screening for large-scale two-sample inference.”  In Royal Statistical Society, volume 81.
  * Chamroukhi and Huynh (2019) Chamroukhi, F. and Huynh, B.-T. (2019).  “Regularized maximum likelihood estimation and feature selection in mixtures-of-experts models.”  Journal de la société française de statistique, 160(1): 57–85.
  * Chao and Fithian (2021) Chao, P. and Fithian, W. (2021).  “AdaPT-GMM: Powerful and robust covariate-assisted multiple testing.”  arXiv preprint arXiv:2106.15812.
  * Chen et al. (2017) Chen, X., Robinson, D. G., and Storey, J. D. (2017).  “The functional false discovery rate with applications to genomics.”  Biostatistics.
  * Chen et al. (2016) Chen, Y., Lun, A. T., and Smyth, G. K. (2016).  “From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline.”  F1000Research, 5.
  * Deb et al. (2021) Deb, N., Saha, S., Guntuboyina, A., and Sen, B. (2021).  “Two-component mixture model in the presence of covariates.”  Journal of the American Statistical Association, 1–35.
  * Dunson et al. (2007) Dunson, D. B., Pillai, N., and Park, J.-H. (2007).  “Bayesian density regression.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(2): 163–183.
  * Efron (2008) Efron, B. (2008).  “Simultaneous inference: When should hypothesis testing problems be combined?”  Ann. Appl. Stat., 2: 197–223.
  * Fan (1991) Fan, J. (1991).  “On the optimal rates of convergence for nonparametric deconvolution problems.”  The Annals of Statistics, 1257–1272.
  * Ferrari and Cribari-Neto (2004) Ferrari, S. and Cribari-Neto, F. (2004).  “Beta regression for modelling rates and proportions.”  Journal of applied statistics, 31(7): 799–815.
  * Harris et al. (2019) Harris, R. M., Kao, H.-Y., Alarcon, J. M., Hofmann, H. A., and Fenton, A. A. (2019).  “Hippocampal transcriptomic responses to enzyme-mediated cellular dissociation.”  Hippocampus, 29(9): 876–882.
  * Heller and Rosset (2021) Heller, R. and Rosset, S. (2021).  “Optimal control of false discovery criteria in the two-group model.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 83(1): 133–155.
  * Himes et al. (2014) Himes, B. E., Jiang, X., Wagner, P., Hu, R., Wang, Q., Klanderman, B., Whitaker, R. M., Duan, Q., Lasky-Su, J., Nikolos, C., et al. (2014).  “RNA-Seq transcriptome profiling identifies CRISPLD2 as a glucocorticoid responsive gene that modulates cytokine function in airway smooth muscle cells.”  PloS one, 9(6): e99625.
  * Ignatiadis and Huber (2021) Ignatiadis, N. and Huber, W. (2021).  “Covariate powered cross-weighted multiple testing.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 83(4): 720–751.
  * Ignatiadis et al. (2016) Ignatiadis, N., Klaus, B., Zaugg, J. B., and Huber, W. (2016).  “Data-driven hypothesis weighting increases detection power in genome-scale multiple testing.”  Nature methods, 13(7): 577–580.
  * Ji et al. (2005) Ji, Y., Wu, C., Liu, P., Wang, J., and Coombes, K. R. (2005).  “Applications of beta-mixture models in bioinformatics.”  Bioinformatics, 21(9): 2118–2122.
  * Kelly et al. (2010) Kelly, R. C., Smith, M. A., Kass, R. E., and Lee, T. S. (2010).  “Local field potentials indicate network state and account for neuronal response variability.”  Journal of computational neuroscience, 29(3): 567–579.
  * Korthauer et al. (2019) Korthauer, K., Kimes, P. K., Duvallet, C., Reyes, A., Subramanian, A., Teng, M., Shukla, C., Alm, E. J., and Hicks, S. C. (2019).  “A practical guide to methods controlling false discoveries in computational biology.”  Genome biology, 20(1): 1–21.
  * Lei and Fithian (2018) Lei, L. and Fithian, W. (2018).  “AdaPT: an interactive procedure for multiple testing with side information.”  J. R. Stat. Soc. Ser. B. Stat. Methodol., 80(4): 649–679.
  * Lei et al. (2017) Lei, L., Ramdas, A., and Fithian, W. (2017).  “STAR: A general interactive framework for FDR control under structural constraints.”  arXiv preprint arXiv:1710.02776.
  * Leung (2022) Leung, D. (2022).  “$Z$-value Directional False Discovery Rate Control with Data Masking.”  arXiv preprint arXiv:2201.05828.
  * Li and Barber (2019) Li, A. and Barber, R. F. (2019).  “Multiple testing with the structure-adaptive Benjamini–Hochberg algorithm.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 81(1): 45–74.
  * Love et al. (2014) Love, M. I., Huber, W., and Anders, S. (2014).  “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.”  Genome biology, 15(12): 1–21.
  * Markitsis and Lai (2010) Markitsis, A. and Lai, Y. (2010).  “A censored beta mixture model for the estimation of the proportion of non-differentially expressed genes.”  Bioinformatics, 26(5): 640–646.
  * McLachlan et al. (2006) McLachlan, G. J., Bean, R., and Jones, L. B.-T. (2006).  “A simple implementation of a normal mixture approach to differential gene expression in multiclass microarrays.”  Bioinformatics, 22(13): 1608–1615.
  * Migliorati et al. (2018) Migliorati, S., Di Brisco, A. M., Ongaro, A., et al. (2018).  “A new regression model for bounded responses.”  Bayesian Analysis, 13(3): 845–872.
  * Nguyen et al. (2016) Nguyen, H. D., Lloyd-Jones, L. R., and McLachlan, G. J. (2016).  “A universal approximation theorem for mixture-of-experts models.”  Neural computation, 28(12): 2585–2593.
  * Nguyen et al. (2018) Nguyen, H. D., Yee, Y., McLachlan, G. J., and Lerch, J. P. (2018).  “False discovery rate control under reduced precision computation for analysis of neuroimaging data.”  arXiv preprint arXiv:1805.04394.
  * Parker and Rothenberg (1988) Parker, R. and Rothenberg, R. (1988).  “Identifying important results from multiple statistical tests.”  Statistics in medicine, 7(10): 1031–1043.
  * Pounds and Morris (2003) Pounds, S. and Morris, S. W. (2003).  “Estimating the occurrence of false positives and false negatives in microarray studies by approximating and partitioning the empirical distribution of p-values.”  Bioinformatics, 19(10): 1236–1242.
  * Resnick (2019) Resnick, S. (2019).  A probability path.  Springer.
  * Ritchie et al. (2015) Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., and Smyth, G. K. (2015).  “limma powers differential expression analyses for RNA-sequencing and microarray studies.”  Nucleic acids research, 43(7): e47–e47.
  * Scott et al. (2015) Scott, J. G., Kelly, R. C., Smith, M. A., Zhou, P., and Kass, R. E. (2015).  “False discovery rate regression: an application to neural synchrony detection in primary visual cortex.”  Journal of the American Statistical Association, 110(510): 459–471.
  * Smith and Kohn (2008) Smith, M. A. and Kohn, A. (2008).  “Spatial and temporal scales of neuronal correlation in primary visual cortex.”  Journal of Neuroscience, 28(48): 12591–12603.
  * Storey (2002) Storey, J. D. (2002).  “A direct approach to false discovery rates.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(3): 479–498.
  * Storey et al. (2007) Storey, J. D., Dai, J. Y., and Leek, J. T. (2007).  “The optimal discovery procedure for large-scale significance testing, with applications to comparative microarray experiments.”  Biostatistics, 8(2): 414–432.
  * Storey et al. (2004) Storey, J. D., Taylor, J. E., and Siegmund, D. (2004).  “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.
  * Sun and Cai (2007) Sun, W. and Cai, T. T. (2007).  “Oracle and adaptive compound decision rules for false discovery rate control.”  Journal of the American Statistical Association, 102(479): 901–912.
  * Tansey et al. (2018) Tansey, W., Koyejo, O., Poldrack, R. A., and Scott, J. G. (2018).  “False discovery rate smoothing.”  Journal of the American Statistical Association, 113(523): 1156–1171.
  * Tian et al. (2021) Tian, Z., Liang, K., and Li, P. (2021).  “A powerful procedure that controls the false discovery rate with directional information.”  Biometrics, 77(1): 212–222.
  * Varadhan and Roland (2008) Varadhan, R. and Roland, C. (2008).  “Simple and globally convergent methods for accelerating the convergence of any EM algorithm.”  Scandinavian Journal of Statistics, 35(2): 335–353.
  * White (1981) White, H. (1981).  “Consequences and detection of misspecified nonlinear regression models.”  Journal of the American Statistical Association, 76(374): 419–433.
  * White (1982) — (1982).  “Maximum likelihood estimation of misspecified models.”  Econometrica: Journal of the Econometric Society, 1–25.
  * Yurko et al. (2020) Yurko, R., G?Sell, M., Roeder, K., and Devlin, B. (2020).  “A selective inference approach for false discovery rate control using multiomics covariates yields insights into disease risk.”  Proceedings of the National Academy of Sciences, 117(26): 15028–15035.
  * Zhang and Chen (2020) Zhang, X. and Chen, J. (2020).  “Covariate Adaptive False Discovery Rate Control with Applications to Omics-Wide Multiple Testing.”  Journal of the American Statistical Association, (just-accepted): 1–31.

##  Appendix A Oracle procedures

###  A.1. Optimality of $\boldsymbol{\delta}^{\mathcal{P}}$ and
$\boldsymbol{\delta}^{\mathcal{Z}}$

In this section we briefly review the optimality properties of the oracle
procedures $\boldsymbol{\delta}^{\mathcal{P}}$ and
$\boldsymbol{\delta}^{\mathcal{Z}}$ in Section 2.2 and how their thresholds
$t_{\mathcal{P}}$ and $t_{\mathcal{Z}}$ are determined. To streamline the
discussion, we will let $\\{M_{i}\\}_{i=1}^{m}$ denote the set of _main
statistics_ , where it can either be that $M_{i}=P_{i}$ or $M_{i}=Z_{i}$ for
all $i$, depending on whether $p$-value or $z$-value based methods are
considered. Given the data $\\{M_{i},X_{i}\\}_{i=1}^{m}$, it is well-known
that _optimal_ procedures, which aim to maximize true discoveries subject to
false discovery constraints, should operate by rejecting $i$ if its
corresponding posterior probability $P(H_{i}=0|M_{i},X_{i})$ falls below a
data-dependent threshold $t_{\mathcal{M}}$. We will use
$\boldsymbol{\delta}^{\mathcal{M}}$ (in a similar way as
$\boldsymbol{\delta}^{\mathcal{P}}$ or $\boldsymbol{\delta}^{\mathcal{Z}}$) to
denote the procedure that thresholds the quantities
$\\{P(H_{i}=0|M_{i},X_{i})\\}_{i=1}^{m}$ with $t_{\mathcal{M}}$.

There are subtly different ways to define “optimality”, depending on the
particular false discovery (e.g. FDR, mFDR, pFDR) and power (e.g. TPR, ETD,
mFNR) measures used, which may lead to different ways of setting
$t_{\mathcal{M}}$. A most recent result of Heller and Rosset (2021, Theorem
3.1) suggests that among all the testing procedures that are functions of
$\\{M_{i},X_{i}\\}_{i=1}^{m}$, the $t_{\mathcal{M}}$ that renders an ETD-
maximing $\boldsymbol{\delta}^{\mathcal{M}}$ with $FDR\leq\alpha$ can be found
by solving an integer optimization problem (Heller and Rosset, 2021, Theorem
3.1). For our purpose, by letting $L_{(1)}<\dots<L_{(m)}$ be the order
statistics of the posterior probabilities
$\\{P(H_{i}=0|M_{i},X_{i})\\}_{i=1}^{m}$, we have considered the
computationally simpler optimal procedure first proposed in Sun and Cai (2007)
which takes $t_{\mathcal{M}}=L_{(j)}$, where

(A.1)
$j\equiv\max\left\\{i^{\prime}\in\\{1,\dots,m\\}:\frac{\sum_{i=1}^{i^{\prime}}L_{(i)}}{i^{\prime}}\leq\alpha\right\\}.$

This procedure has $FDR\leq\alpha$ because for any procedure that produces a
rejection set $\mathcal{R}$ based on $\\{M_{i},X_{i}\\}_{i=1}^{m}$, its FDR
can be written as

$FDR=\mathbb{E}\left[\underbrace{\mathbb{E}\left[\frac{\sum_{i\in\mathcal{R}}(1-H_{i})}{|\mathcal{R}|\vee
1}\Big{|}\\{M_{i},X_{i}\\}_{i=1}^{m}\right]}_{\mathbb{E}[FDP|\\{M_{i},X_{i}\\}_{i=1}^{m}]}\right]=\mathbb{E}\left[\frac{\sum_{i\in\mathcal{R}}P(H_{i}=0|M_{i},X_{i})}{|\mathcal{R}|\vee
1}\right].$

Conditional on any instance of the data $\\{M_{i},X_{i}\\}_{i=1}^{m}$,
$\boldsymbol{\delta}^{\mathcal{M}}$ prioritizes rejections of the hypotheses
that are least likely to be true nulls, all the while controlling the
conditional FDR $\mathbb{E}[FDP|\\{M_{i},X_{i}\\}_{i=1}^{m}]$ below $\alpha$
by setting $t_{\mathcal{M}}=L_{(j)}$, as the ratio
$\sum_{i=1}^{i^{\prime}}L_{(i)}/i^{\prime}$ in (A.1) is precisely the
conditional FDR of rejecting the most promising $i^{\prime}$ hypotheses. As a
result, its controls the FDR under $\alpha$ since the conditional version
$\mathbb{E}[FDP|\\{M_{i},X_{i}\\}_{i=1}^{m}]$ is always not larger than
$\alpha$.

We now give a more precise account of the optimal property of the prior
procedure. Another popular measure of type 1 errors is the marginal FDR
(mFDR), which for any rejection set $\mathcal{R}$ is the ratio

$mFDR\equiv\frac{\mathbb{E}[V]}{\mathbb{E}[R]},$

where $V$ and $R$ are defined as in the main text. For $\alpha\in(0,1)$, it is
known that among all the procedures based on $\\{M_{i},X_{i}\\}_{i=1}^{m}$
with $mFDR\leq\alpha$, the ETD-maximizing procedure is the one given by
$\boldsymbol{\delta}^{\mathcal{M}}$ that sets
$t_{\mathcal{M}}=\lambda^{*}_{\mathcal{M}}$, where

(A.2)
$\lambda^{*}_{\mathcal{M}}\equiv\sup\left\\{\lambda\in(0,1]:\frac{\sum_{i}\mathbb{E}[(1-H_{i})I(P(H_{i}=0|M_{i},X_{i})\leq\lambda)]}{\sum_{i}\mathbb{E}[I(P(H_{i}=0|M_{i},X_{i})\leq\lambda)]}\leq\alpha\right\\};$

see Sun and Cai (2007) and Heller and Rosset (2021). In practice, since
$\lambda^{*}_{\mathcal{M}}$ could be tricky to obtain even with oracle
knowledge, this optimal procedure for mFDR control is often approximated by
our computationally handy version with $t_{\mathcal{M}}=L_{(j)}$ above when
$m$ is large, as is the case with most FDR analyses. Their asymptotic
equivalence can be shown by standard arguments, such as those in Sun and Cai
(2007, Section 4). The aforementioned references provide a more detailed
exposition.

###  A.2. Rejection regions for Examples 2.1-2.3

Consider the p-value conditional mixture density

(A.3) $P_{i}|X_{i}=x\sim g_{x}(p)\equiv
g(p|x)=(1-w_{x})g_{0}(p)+w_{x}g_{1,x}(p),$

induced by (2.2), where $g_{0}\equiv 1$ is the uniform null density of
$P_{i}$, and $g_{1,x}(p)\equiv g(p|H_{i}=1,X_{i}=x)$ is the conditional
alternative density of $P_{i}$. The rejection regions
$\mathcal{S}^{\mathcal{P}}(x)$ and $\mathcal{S}^{\mathcal{Z}}(x)$ in Figure
2.1 are derived based on the threshold $\lambda^{*}_{\mathcal{M}}$ in (A.2),
where

$\mathcal{S}^{\mathcal{Z}}(x)\equiv\left\\{z:P(H_{i}=0|Z_{i}=z,X_{i}=x)\leq\lambda^{*}_{\mathcal{Z}}\right\\}=\left\\{z:\frac{w_{x}f_{1,x}(z)}{(1-w_{x})f_{0}(z)}\leq\frac{1-\lambda^{*}_{\mathcal{Z}}}{\lambda^{*}_{\mathcal{Z}}}\right\\}$

and

$\mathcal{S}^{\mathcal{P}}(x)\equiv\left\\{z:P(H_{i}=0|P_{i}=p,X_{i}=x)\leq\lambda^{*}_{\mathcal{P}}\right\\}=\left\\{z:\frac{w_{x}g_{1,x}(2\Phi(-|z|))}{(1-w_{x})g_{0}(2\Phi(-|z|))}\leq\frac{1-\lambda^{*}_{\mathcal{P}}}{\lambda^{*}_{\mathcal{P}}}\right\\}.$

$\lambda^{*}_{\mathcal{P}}$ and $\lambda^{*}_{\mathcal{Z}}$ are
$\lambda^{*}_{\mathcal{M}}$ defined with $M_{i}=P_{i}$ and $M_{i}=Z_{i}$ for
all $i$ respectively. Since the examples are relatively simple, these regions
can be found by numerical means, and we will derive
$\mathcal{S}^{\mathcal{Z}}(x)$ in Example Example 2.2 as an illustration: It
is the set

$\left\\{z:\frac{0.2\left[\frac{1-x}{2}\exp(-\frac{(z+\mu)^{2}}{2})+\frac{1+x}{2}\exp(-\frac{(z-\mu)^{2}}{2})\right]}{0.8\exp(-\frac{z^{2}}{2})}\leq\frac{1-\lambda^{*}_{\mathcal{Z}}}{\lambda^{*}_{\mathcal{Z}}}\right\\},$

where $\mu\equiv 1.5$. By setting
$\frac{w_{x}f_{1,x}(z)}{(1-w_{x})f_{0}(z)}=\frac{1-\lambda^{*}_{\mathcal{Z}}}{\lambda^{*}_{\mathcal{Z}}}$,
one arrives at the equation

$\frac{1+x}{2}\exp(2\mu
z)-4\frac{1-\lambda^{*}_{\mathcal{Z}}}{\lambda^{*}_{\mathcal{Z}}}\exp\left(\frac{\mu^{2}}{2}\right)\exp(\mu
z)+\frac{1-x}{2}=0.$

in $z$. To solve for a solution $z^{*}$, we can apply the formula for the
solutions of a quadratic equation to get

$\exp(\mu
z^{*})=\frac{4\frac{1-\lambda_{\mathcal{Z}}^{*}}{\lambda_{\mathcal{Z}}^{*}}\exp(\frac{\mu^{2}}{2})\pm\sqrt{16\left(\frac{1-\lambda_{\mathcal{Z}}^{*}}{\lambda_{\mathcal{Z}}^{*}}\right)^{2}\exp(\mu^{2})-(1-x)(1+x)}}{1+x},$

which in turn implies the two boundary points

$z^{*}=\frac{1}{\mu}\log\left\\{\frac{4\frac{1-\lambda_{\mathcal{Z}}^{*}}{\lambda_{\mathcal{Z}}^{*}}\exp(\frac{\mu^{2}}{2})\pm\sqrt{16\left(\frac{1-\lambda_{\mathcal{Z}}^{*}}{\lambda_{\mathcal{Z}}^{*}}\right)^{2}\exp(\mu^{2})-(1-x)(1+x)}}{1+x}\right\\}$

for the red regions in the middle panel of Figure 2.1 as a function of $x$.

We now explain why $\mathcal{S}^{\mathcal{P}}(x)$ doesn’t change with $x$ in
Examples 2.2 and 2.3. From the conditional $p$-value density (A.3), one can
see that

$P(H_{i}=0|P_{i},X_{i})=\frac{(1-w_{X_{i}})g_{0}(P_{i})}{g_{X_{i}}(P_{i})}=\frac{(1-w_{X_{i}})g_{0}(P_{i})}{(1-w_{X_{i}})g_{0}(P_{i})+w_{X_{i}}g_{1,X_{i}}(P_{i})}.$

Hence, if $w_{x}$ and $g_{1,x}$ do not depend on $x$, it is apparent that
$\mathcal{S}_{x}^{\mathcal{P}}$ will not vary in $x$. Simple calculations can
show that $w_{x}=0.2$ and $w_{x}=0.1$ for Examples 2.2 and 2.3 respectively,
and

$g_{1,x}(p)=\frac{\phi(\Phi^{-1}(p/2)-1.5)+\phi(-\Phi^{-1}(p/2)-1.5)}{2\phi(-\Phi^{-1}(p/2))}=\frac{\phi(-|z|-1.5)+\phi(|z|-1.5)}{2\phi(|z|)}$

for both examples; all of these quantities do not depend on $x$.

##  Appendix B Proof for the prototype method

We will prove the FDR validity of the prototype procedure in Section 3.1. The
false discovery proportion of the prototype testing procedure, which
thresholds the test statistics $T_{i}$’s with the threshold
$\hat{t}_{\alpha}\equiv\hat{t}(\alpha)$, can be written as

$\displaystyle FDP$ $\displaystyle=\frac{\\#\\{i\text{ null:
}T_{i}\leq\hat{t}_{\alpha}\\}}{1\vee\\#\\{T_{i}\leq\hat{t}_{\alpha}\\}}$
$\displaystyle=\frac{\\#\\{i\text{ null:
}T_{i}\leq\hat{t}_{\alpha}\\}}{1+\\#\\{i\text{ null: }S_{i}\geq
1-c_{i}(\hat{t}_{\alpha})\\}}\underbrace{\frac{1+\\#\\{i\text{ null:
}S_{i}\geq
1-c_{i}(\hat{t}_{\alpha})\\}}{1\vee\\#\\{T_{i}\leq\hat{t}_{\alpha}\\}}}_{\leq\alpha\text{
by the defintion of our procedure}}$
$\displaystyle\leq\alpha\frac{\\#\\{i\text{ null:
}T_{i}\leq\hat{t}_{\alpha}\\}}{1+\\#\\{i\text{ null: }S_{i}\geq
1-c_{i}(\hat{t}_{\alpha})\\}}.$

We only have to show that $\mathbb{E}\left[\frac{\\#\\{i\text{ null:
}T_{i}\leq\hat{t}_{\alpha}\\}}{1+\\#\\{i\text{ null: }S_{i}\geq
1-c_{i}(\hat{t}_{\alpha})\\}}\right]$ is bounded by 1 using the stopping time
argument from (Barber and Candès, 2019).

Without loss of generality we will assume the true nulls are the first $m_{0}$
hypotheses. For each $i\in\\{1,\dots,m_{0}\\}$, define

$\widecheck{S}_{i}=\begin{dcases*}S_{i}&when $S_{i}\leq 0.5$\\\ 1-S_{i}&when
$S_{i}>0.5$\end{dcases*}$

and $\widecheck{T}_{i}\equiv c_{i}^{-1}(\widecheck{S}_{i})$. Using the order
statistics $\widecheck{T}_{(1)}\leq\dots\leq\widecheck{T}_{(m_{0})}$ of
$\widecheck{T}_{1},\dots,\widecheck{T}_{m_{0}}$, we moreover let $B_{i}\equiv
I(S_{(i)}>0.5)$ for $i=1,\dots,m_{0}$, where the order of $S_{(i)}$’s here is
inherited from the order of the $\widecheck{T}_{(i)}$’s, rather than the
magnitudes of the $S_{i}$’s themselves. Let $1\leq J\leq m_{0}$ be the index
such that

$\widecheck{T}_{(1)}\leq\dots\leq\widecheck{T}_{(J)}\leq\hat{t}_{\alpha}<\widecheck{T}_{(J+1)}\leq\dots\leq\widecheck{T}_{(m_{0})},$

we then have

$\frac{\\#\\{i\text{ null: }T_{i}\leq\hat{t}_{\alpha}\\}}{1+\\#\\{i\text{
null: }S_{i}\geq 1-c_{i}(\hat{t}_{\alpha})\\}}=\frac{\\#\\{i\text{ null:
}T_{i}\leq\hat{t}_{\alpha}\\}}{1+\\#\\{i\text{ null:
}c_{i}^{-1}(1-S_{i})\leq\hat{t}_{\alpha}\\}}\\\
=\frac{(1-B_{1})+\dots+(1-B_{J})}{1+B_{1}+\dots+B_{J}}=\frac{1+J}{1+B_{1}+\dots+B_{J}}-1,$

considering that $\hat{t}_{\alpha}$ must be less than
$t_{\max}:=\max\\{t:c_{i}(t)\leq 0.5\text{ for all }i\\}$. Hence it amounts to
showing $E[\frac{1+J}{1+B_{1}+\dots+B_{J}}]\leq 2$. This final step can be
shown by applying Barber and Candès (2019, Lemma 1), since _conditional on_
(i) $\widecheck{T}_{(1)}\dots\widecheck{T}_{(m_{0})}$and  (ii)
$\\{T_{i}:i\text{ is non null}\\}$ , $B_{1},\dots,B_{m_{0}}$ are independent
Bernoulli$(0.5)$ random variables, and $J$ can be seen as a stopping time in
reverse time with respect to the filtrations
$\\{\mathcal{F}_{j}\\}_{j=1}^{m_{0}}$, where
$\mathcal{F}_{j}\equiv\\{B_{1}+\dots+B_{j},B_{j+1},\dots,B_{m_{0}}\\}$.

##  Appendix C Properties of the working model

In this section we will develop some properties of the beta-mixture model in
Section 3.2 and the assessor functions it induces. To simplify notation, we
will use $\pi_{li},\pi_{ri},k_{li},k_{ri},h_{li},h_{ri}$ to respectively
denote the quantities and functions $\pi_{l,X_{i}}$, $\pi_{r,X_{i}}$,
$k_{l,X_{i}}$, $k_{r,X_{i}}$, $h_{l,X_{i}}$, $h_{r,X_{i}}$ from Model (3.6)
when the observed covariate $X_{i}$ is used, where the underlying parameters
$\\{\theta_{l},\theta_{r},\beta_{l},\beta_{r}\\}$ are unspecified but common
for all $i=1,\dots,m$. Likewise, we also use

(C.1) $a_{i}(u)\equiv
a_{X_{i}}(u)=\frac{1-\pi_{li}-\pi_{ri}}{(1-\pi_{li}-\pi_{ri})+\pi_{li}h_{li}(u)+\pi_{ri}h_{ri}(u)}$

to denote the assessor function constructed with them, and $T_{i}$, $S_{i}$
and $c_{i}(\cdot)$ will denote the test statistics and null distribution
function based on $a_{i}(\cdot)$ as in Section 3.1.

First, the following lemma states properties concerning the left and right
alternative functions $h_{li}$ and $h_{ri}$.

######  Lemma C.1 (Properties of the non-null component densities).

For $\gamma_{l}>2$, $h_{li}(\cdot)$ is a strictly convex and strictly
decreasing function with the properties

$\lim_{u\rightarrow 0}h_{li}(u)=\infty\text{ and }\lim_{u\rightarrow
1}h_{li}(u)=0.$

Similarly, for $\gamma_{r}>2$, $h_{r_{i}}(\cdot)$ is a strictly convex and
strictly increasing function with the properties

$\lim_{u\rightarrow 0}h_{ri}(u)=0\text{ and }\lim_{u\rightarrow
1}h_{ri}(u)=\infty.$

###### Proof of Lemma C.1.

It suffices to show the facts for $h_{li}$ since those for $h_{ri}$ can be
proven exactly analogously. Recall that

$h_{li}(u):=B(k_{li},\gamma_{l})^{-1}u^{k_{li}-1}(1-u)^{\gamma_{l}-1},$

where for brevity we have suppressed the dependence on $X_{i}$ in notations.
Differentiating with respect to $u$ we get

(C.2)
$h_{li}^{\prime}(u)=B(k_{li},\gamma_{l})^{-1}[(k_{li}-1)u^{k_{li}-2}(1-u)^{\gamma_{l}-1}-({\gamma_{l}-1})(1-u)^{\gamma_{l}-2}u^{k_{li}-1}],$

which, given $\gamma_{l}>2$ (actually $\gamma_{l}>1$ is sufficient), can be
seen to be always negative for any $u\in(0,1)$ and hence proves that $h_{li}$
is strictly decreasing. For convexity, we differentiate one more time to get

$\displaystyle h_{li}^{\prime\prime}(u)$
$\displaystyle=B(k_{l},\gamma_{l})(k_{li}-1)\left[(k_{li}-2)u^{k_{li}-3}(1-u)^{\gamma_{l}-1}-u^{k_{li}-2}(\gamma_{l}-1)(1-u)^{\gamma_{l}-2}\right]$
$\displaystyle\quad-B(k_{l},\gamma_{l})(\gamma_{l}-1)\left[(2-\gamma_{l})(1-u)^{\gamma_{l}-3}u^{k_{li}-1}+(k_{li}-1)u^{k_{li}-2}(1-u)^{\gamma_{l}-2}\right]$
$\displaystyle=B(k_{l},\gamma_{l})\underbrace{u^{k_{li}-3}(1-u)^{\gamma_{l}-3}}_{>0}\times$
$\displaystyle\qquad\left[\underbrace{(k_{li}-1)(k_{li}-2)(1-u)^{2}}_{>0}-\underbrace{2(k_{li}-1)(\gamma_{l}-1)u(1-u)}_{>0}+(\gamma_{l}-1)(\gamma_{l}-2)u^{2}\right],$

where in the last equality, the positive terms are positive since
$0<k_{li}<1$. As such, $h_{li}^{\prime\prime}(u)$ is strictly positive for all
$u\in(0,1)$ as long as $\gamma_{l}>2$, which proves the strict convexity of
$h_{li}$. ∎

A closer inspection of the proof above will reveal that for
$\gamma_{l}\in(1,2)$, $h_{li}$ may not even be convex, and the same is true
for $h_{ri}$. Hence we have required that $\gamma_{l},\gamma_{r}>2$ in our
model. To facilitate the proof in later sections we will also define the
_reciprocal_ assessor function

(C.3) $b_{i}(u)\equiv 1/a_{i}(u).$

By the properties of $h_{li}$ and $h_{ri}$ in Lemma C.1, one can readily
conclude the following lemma, which will help us develop some useful facts
later:

######  Lemma C.2 (Properties of the reciprocal assessor).

The reciprocal assessor function defined in (C.3) (for
$\gamma_{l},\gamma_{r}>2$) is strictly convex and smooth, with the property
that

(C.4) $\lim_{u\rightarrow 0}b_{i}(u)=\lim_{u\rightarrow 1}b_{i}(u)=\infty.$

Hence, there exists a unique minimal $\underline{u}_{i}$ such that

$b_{i}(\underline{u}_{i})<b_{i}(u)\text{ for all }u\in(0,1).$

$T_{i}\equiv a_{i}(U_{i})$ as a random variable has the range
$(0,a_{i}(\underline{u}_{i})]$ in light of Lemma C.2. By construction,
$a_{i}(\cdot)$’s level sets can only be of Lebesgue measure 0, so
$c_{i}(\cdot)$ is continuous under the uniform null distribution of $U_{i}$.
With the strict convexity of $b_{i}(\cdot)$, one can also conclude that
$c_{i}(\cdot)$, its null distribution function, is invertible (or
equivalently, strictly increasing), since no interval in the range
$(0,a_{i}(\underline{u}_{i})]$ will have zero measure under the law of
$a_{i}(U_{i})$ induced by the uniform null distribution of $U_{i}$ by the
intermediate value theorem. The smooth “bowl” shape of $b_{i}(\cdot)$ also
implies that, for any $t\in(0,a_{i}(\underline{u}_{i})]$, the event
$\\{T_{i}>t\\}$ is equivalent to $U_{i}$ taking values in a certain sub-
interval of $(0,1)$. One can define two smooth functions to describe this
fact:

######  Definition C.1 (Expression for the event $\\{T_{i}>t\\}$).

For each $i$,
$\omega_{iL}:(0,a_{i}(\underline{u}_{i})]\rightarrow(0,\underline{u}_{i}]$ and
$\omega_{iR}:(0,a_{i}(\underline{u}_{i})]\rightarrow[\underline{u}_{i},1)$ are
respectively two smooth functions such that for any
$t\in(0,a_{i}(\underline{u}_{i})]$,

$\\{T_{i}>t\\}=\\{\omega_{iL}(t)<U_{i}<\omega_{iR}(t)\\}=\\{S_{i}>c_{i}(t)\\},$

with $\omega_{iL}(\cdot)$ and $\omega_{iR}(\cdot)$ being strictly increasing
and strictly decreasing, respectively. Since $S_{i}$ and $U_{i}$ are uniformly
distributed when $H_{i}=0$,

$\omega_{iR}(t)-\omega_{iL}(t)=1-c_{i}(t).$

Moreover, $b_{i}(\omega_{iL}(t))=b_{i}(\omega_{iR}(t))=1/t$.

Of course the variable $S_{i}=c_{i}(T_{i})$ has the range $(0,1]$, and we can
define functions to describe events of the form $\\{S_{i}>s\\}$ similar to
Definition C.1:

######  Definition C.2 (Expression for the event $\\{S_{i}>s\\}$).

For each $i$, $\psi_{iL}:(0,1]\rightarrow(0,\underline{u}_{i}]$ and
$\psi_{iR}:(0,1]\rightarrow[\underline{u}_{i},1)$ are respectively two smooth
functions such that for any $s\in(0,1]$,

$\\{S_{i}>s\\}=\\{\psi_{iL}(s)<U_{i}<\psi_{iR}(s)\\}=\\{T_{i}>c_{i}^{-1}(s)\\},$

with $\psi_{iL}(\cdot)$ and $\psi_{iR}(\cdot)$ being strictly increasing and
strictly decreasing, respectively. Since $S_{i}$ and $U_{i}$ are uniformly
distributed when $H_{i}=0$,

(C.5) $\psi_{iR}(s)-\psi_{iL}(s)=1-s.$

Moreover, $b_{i}(\psi_{iL}(s))=b_{i}(\psi_{iR}(s))=1/c_{i}^{-1}(s)$.

##  Appendix D Proof for the asymptotic method

Before proving Theorem 3.1, we remark that the theorem is established by
assuming that the mirror statistic $\hat{T}^{\mathfrak{m}}_{i}$ is the _exact_
reflection of $\hat{T}_{i}$ under the null distribution $\hat{c}_{i}$. In
practice, $\hat{T}^{\mathfrak{m}}_{i}$ can be determined up to arbitrary
precision in Algorithm 1 as long as the number of uniform realizations $N$ is
set to be very large, as recommended in Section 3.5.

We will make heavy use of the notation and results in Appendix C. We will also
use $a^{*}_{i}(\cdot)$, $b^{*}_{i}(\cdot)$, $c^{*}_{i}(\cdot)$, $T_{i}^{*}$
and $S_{i}^{*}$ to denote the respective functions and statistics when
$\\{\theta,\beta\\}$ is taken to be the pair $\\{\theta^{*},\beta^{*}\\}$ to
construct $a_{i}(\cdot)$, $b_{i}(\cdot)$, $c_{i}(\cdot)$, $T_{i}$ and $S_{i}$.
Similarly, the quantities and functions appearing in Lemma C.2 and Definitions
C.1 and C.2 all have their “star” versions: $\underline{u}^{*}_{i}$,
$\omega^{*}_{iL}(\cdot)$, $\omega^{*}_{iR}(\cdot)$, $\psi^{*}_{iL}(\cdot)$ and
$\psi^{*}_{iR}(\cdot)$. Generally speaking, $C,c>0$ will denote unspecified
universal constants required for the asymptotic arguments in this section.

###  D.1. Additional assumptions for Theorem 3.1

######  Assumption 2 (Regularity conditions).

  1. (i)

$\max_{i}\|X_{i}\|_{\infty}\leq C$ almost surely for some universal constant
$C>0$, where $\|\cdot\|_{\infty}$ indicates the sup norm.

  2. (ii)

Let $[\delta_{1},\delta_{2}]$ be any fixed compact interval in $(0,1)$. For
each $i$, let $V_{1i}$ and $V_{2i}$ be two measurable subsets in
$(\delta_{1},\delta_{2})$. Then for large enough $m$,

$\left|\frac{1}{m}\sum_{i=1}^{m}P(U_{i}\in V_{1i}|X_{i})-P(U_{i}\in
V_{2i}|X_{i})\right|\leq C(\delta_{1},\delta_{2})\max_{1\leq i\leq
m}\lambda\left(V_{i1}\Delta V_{i2}\right),$

where $\lambda(\cdot)$ is the Lebesgue measure, $C(\delta_{1},\delta_{2})$ is
a constant that may depends on $\delta_{1}$ and $\delta_{2}$, and
$\mathcal{S}_{1}\Delta\mathcal{S}_{2}$ is their symmetric difference for any
two sets $\mathcal{S}_{1},\mathcal{S}_{2}\subset\mathbb{R}$.

  3. (iii)

$\mathbb{E}[\log(U_{i})]$ and $\mathbb{E}[\log(1-U_{i})]$ are finite.

######  Assumption 3 (Strong laws of large numbers).

Let $\pi_{0}\equiv P(H_{i}=0)>0$. For any $t\in(0,1]$, it holds that

(D.1) $m^{-1}\sum_{i=1}^{m}I(T_{i}^{*}\leq t)\longrightarrow G(t),\\\
m^{-1}\sum_{i=1}^{m}(1-H_{i})I(T_{i}^{*}\leq
t)\longrightarrow\pi_{0}G_{0}(t)\text{ and }\\\
m^{-1}\sum_{i=1}^{m}I(S_{i}^{*}\geq
1-c^{*}_{i}(t))\longrightarrow\bar{G}_{0}(t)$

almost surely, where $G(t)$, $G_{0}(t)$ and $\bar{G}_{0}(t)$ are positive
continuous functions in $t$. Moreover, $\bar{G}_{0}(t_{0})/G(t_{0})<\alpha$
for some $t_{0}>0$, and the limiting threshold

$t_{\alpha}^{\infty}\equiv\sup\left\\{t\in(0,1]:\frac{\bar{G}_{0}(t)}{G(t)}\leq\alpha\right\\}$

is such that $\max_{i\in\mathbb{N}}c_{i}^{*}(t_{\alpha}^{\infty})<\bar{s}$ for
$\bar{s}<1$, or equivalently,
$t_{\alpha}^{\infty}<\min_{i}{c_{i}^{*}}^{-1}(\bar{s})$. Note that the strong
laws above, as well as the marginal probability $\pi_{0}$, are with respect to
the joint law of $\\{H_{i},Z_{i},X_{i}\\}$.

Assumption 2 regulates the tail behaviors of the random variables $X_{i}$ and
$U_{i}$; in particular, $(ii)$ implies that conditional on $X_{i}$, the
density of $U_{i}$ can be unbounded at the two tails, which is natural for
multiple testing as it provides room for non-null tail behaviors. Assumption 3
states properties of the strong law limits involved. For technical reasons,
that $\sup_{i\in\mathbb{N}}c_{i}^{*}(t_{\alpha}^{\infty})$ is bounded away
from 1 ensures that our result won’t rely on the strong law limits for very
large values of $t$, which is hardly restrictive in practice: any sensible
multiple testing procedure should only consider rejecting $S_{i}^{*}$, which
is uniformly distributed under the null, if it is much less than the typically
small target FDR level $\alpha$. Similar assumptions have also appeared in the
works of Zhang and Chen (2020), Storey et al. (2004).

We remark that our current assumptions for Theorem 3.1 are no stronger than
those in Zhang and Chen (2020) in any essential way, and can conceivably be
further relaxed; for example, if Assumption 2$(i)$ is phrased as a
probabilistic bound, one can still likely establish a version of Theorem 3.1
which says that the FDR is less than $\alpha$ with probability approaching 1.
Moreover, we have assume, as stated in Section 2.1, that
$\\{H_{i},Z_{i},X_{i}\\}$ are independent across $i$, which can be further
relaxed to a generic weak dependence condition under which the strong laws in
Assumption 3 hold. In fact it is possible to prove an FDR bound in terms of
the conditional expectation $\mathbb{E}[\cdot|H_{1},\dots,H_{m}]$, treating
the hypotheses as fixed. These embellishments have not been pursued here for a
more streamlined presentation.

###  D.2. Technical lemmas

Under Assumption 1 and Assumption 2$(i)$, all $k_{ri}=k_{r}(X_{i})$,
$k_{li}=k_{l}(X_{i})$, $\pi_{ri}=\pi_{r}(X_{i})$, $\pi_{li}=\pi_{l}(X_{i})$
are bounded away from one and zero, i.e.
$k_{ri},k_{li}\in[\underline{k},\bar{k}]$,
$\pi_{ri},\pi_{li}\in[\underline{\pi},\bar{\pi}]$ for some compact intervals
$[\underline{k},\bar{k}],[\underline{\pi},\bar{\pi}]\subset(0,1)$, by the
compactness of $\Theta\times{\bf B}$. Moreover, by continuity of the beta
functions we can define

$B_{\max}:=\max_{k\in[\underline{k},\bar{k}]}(B(k,\gamma_{l})\vee(B(\gamma_{r},k))\text{
and
}B_{\min}:=\max_{k\in[\underline{k},\bar{k}]}(B(k,\gamma_{l})\vee(B(\gamma_{r},k)),$

which are both positive numbers.

The following “uniformity” properties will be heavily relied on later:

######  Lemma D.1 (Uniformity properties).

Under Assumption 1 and Assumption 2$(i)$, the following are true for any
$(\theta,\beta)\in\Theta\times\boldsymbol{B}$:

  1. (i)

For any $t_{0}>0$, there exists a $u_{0}=u_{0}(t_{0})>0$ not depending on
$(\theta,\beta)$ such that for all $t\geq t_{0}$,
$(\omega_{iL}(t),\omega_{iR}(t))\subset[u_{0},1-u_{0}]$ for all $i$.

  2. (ii)

For any $s_{0}>0$, there exists a $u_{0}=u_{0}(s_{0})>0$ not depending on
$(\theta,\beta)$ such that for all $s\geq s_{0}$,
$(\psi_{iL}(s),\psi_{iR}(s))\subset[u_{0},1-u_{0}]$ for all $i$.

  3. (iii)

There exists a small positive constant $u_{0}\in(0,0.5)$ not depending on
$(\theta,\beta)$ such that,

$\underline{u}_{i}\in[u_{0},1-u_{0}]$

for all $i$, where $\underline{u}_{i}$ is as in Lemma C.2.

  4. (iv)

$\lim_{\epsilon\rightarrow
0}\max_{i\in\mathbb{N},\theta\in\Theta,\beta\in\boldsymbol{B}}\lambda\Bigl{(}{b_{i}}^{-1}\bigl{(}[b_{i}(\underline{u}_{i}),b_{i}(\underline{u}_{i})+\epsilon\bigr{)}\Bigr{)}=0,$

where $\lambda$ is the Lebesgue measure, $\underline{u}_{i}$ is as in Lemma
C.2 and

${b_{i}}^{-1}(\mathcal{T})\equiv\\{u:{b_{i}}(u)\in\mathcal{T}\\}.$

for any interval $\mathcal{T}$ in $\mathbb{R}$.

###### Proof.

$(i)$: Note that

(D.2)
$b_{i}(u)\geq\max\left(\frac{\underline{\pi}u^{\bar{k}-1}(1-u)^{\gamma_{l}-1}}{(1-2\underline{\pi})B_{\max}},\frac{\underline{\pi}u^{\gamma_{r}-1}(1-u)^{\bar{k}-1}}{(1-2\underline{\pi})B_{\max}}\right)\text{
for all }i,$

which implies $\lim_{u\rightarrow 0}\min_{i}b_{i}(u)=\lim_{u\rightarrow
1}\min_{i}b_{i}(u)=\infty$, since the right hand side of (D.2) tends to
$\infty$ as $u$ tends to 0 or 1. Hence one must be able to find a small enough
$u_{0}>0$ such that $\min(b_{i}(u_{0}),b_{i}(1-u_{0}))>1/t_{0}$ for all $i$,
which implies that $(\omega_{iL}(t),\omega_{iR}(t))\subset[u_{0},1-u_{0}]$ for
all $i$ by Definition C.1. This proves $(i)$.

$(ii)$: Suppose towards a contradiction, such a $u_{0}$ doesn’t exist. Without
loss of generality, we assume there is a subsequence $\\{i_{1},i_{2},\dots\\}$
such that $\lim_{j\rightarrow\infty}\psi_{i_{j}L}(s_{0})=0$. As such,
$\lim_{j\rightarrow\infty}\psi_{i_{j}R}(s_{0})=1-s_{0}$ by the property stated
in Definition C.2, which implies

(D.3) $\limsup_{j}b_{i_{j}}(\psi_{i_{j}R}(s_{0}))\leq
1+\frac{\bar{\pi}(1-s_{0})^{\underline{k}-1}(s_{0})^{\gamma_{l}-1}}{(1-2\bar{\pi})B_{\min}}+\frac{\bar{\pi}(1-s_{0})^{\gamma_{r}-1}(s_{0})^{\underline{k}-1}}{(1-2\bar{\pi})B_{\min}}.$

On the other hand,

(D.4) $\lim_{j}b_{i_{j}}^{*}(\psi^{*}_{i_{j}L}(s_{0}))=\infty$

in consideration of (D.2) and
$\lim_{j\rightarrow\infty}\psi^{*}_{i_{j}L}(s_{0})=0$. (D.3) and (D.4)
together reach a contradiction since it must be that
$\lim_{j}b_{i_{j}}^{*}(\psi^{*}_{i_{j}L}(s_{0}))=\lim_{j}b_{i_{j}}^{*}(\psi^{*}_{i_{j}R}(s_{0}))$
as
$b_{i_{j}}^{*}(\psi^{*}_{i_{j}L}(s_{0}))=b_{i_{j}}^{*}(\psi^{*}_{i_{j}R}(s_{0}))$by
Definition C.2.

$(iii)$: By the fact that that $b_{i}^{\prime}(b_{i}(\underline{u}_{i}))=0$
for all $i$, it suffices to show that

(D.5) $\lim_{u\rightarrow
0}\max_{i,\theta,\beta}b_{i}^{\prime}(u)=-\infty\text{ and }\lim_{u\rightarrow
1}\min_{i,\theta,\beta}b^{\prime}_{i}(u)=\infty.$

Note that

$b_{i}^{\prime}(u)=\frac{\pi_{li}}{1-\pi_{ri}-\pi_{li}}h^{\prime}_{li}(u)+\frac{\pi_{ri}}{1-\pi_{ri}-\pi_{li}}h^{\prime}_{ri}(u),$

where $h^{\prime}_{li}$ has the form

$h_{li}^{\prime}(u)=B(k_{li},\gamma_{l})^{-1}[(k_{li}-1)u^{k_{li}-2}(1-u)^{\gamma_{l}-1}-({\gamma_{l}-1})(1-u)^{\gamma_{l}-2}u^{k_{li}-1}]$

as shown in the proof of Lemma C.1. Define, for $u\in(0,1)$, the functions

$\bar{h}^{\prime}_{l}(u)=B_{\max}^{-1}[(\bar{k}-1)u^{\bar{k}-2}(1-u)^{\gamma_{l}-1}-(\gamma_{l}-1)(1-u)^{\gamma_{l}-2}u^{\bar{k}-1}],$
$\underline{h}^{\prime}_{l}(u)=B_{\min}^{-1}[(\underline{k}-1)u^{\underline{k}-2}(1-u)^{\gamma_{l}-1}-(\gamma_{l}-1)(1-u)^{\gamma_{l}-2}u^{\underline{k}-1}]$

so that $\underline{h}^{\prime}_{l}(u)\leq
h^{\prime}_{li}(u)\leq\bar{h}^{\prime}_{l}(u)<0$ for all $i$. Note that

$\lim_{u\rightarrow 0}\bar{h}^{\prime}_{l}(u)=-\infty\text{ and
}\lim_{u\rightarrow 1}\underline{h}^{\prime}_{l}(u)=0.$

One can similarly define functions $\bar{h}_{r}^{\prime}$ and
$\underline{h}_{r}^{\prime}$ on $(0,1)$ such that
$0<\underline{h}^{\prime}_{r}(u)\leq
h^{\prime}_{ri}(u)\leq\bar{h}^{\prime}_{r}(u)$ for all $i$ and

$\lim_{u\rightarrow 0}\bar{h}^{\prime}_{r}(u)=0\text{ and }\lim_{u\rightarrow
1}\underline{h}^{\prime}_{r}(u)=\infty.$

The fact that
$b_{i}^{\prime}(u)\leq\frac{\bar{\pi}}{1-2\bar{\pi}}\bar{h}_{l}^{\prime}(u)+\frac{\bar{\pi}}{1-2\bar{\pi}}\bar{h}_{r}^{\prime}(u)$,
together with two of the limit results above, has shown the first limit in
(D.5). Similarly, that
$b_{i}^{\prime}(u)\geq\frac{\underline{\pi}}{1-2\underline{\pi}}\underline{h}_{l}^{\prime}(u)+\frac{\underline{\pi}}{1-2\underline{\pi}}\underline{h}_{r}^{\prime}(u)$,
together with the other two limit results above, has shown the second limit in
(D.5).

$(iv)$: By Lemma D.1$(iii)$, pick $u_{0}\in(0,0.5)$ such that
$\underline{u}_{i}\in[u_{0},1-u_{0}]$ for all $i$. Note that

$b_{i}^{\prime\prime}(u)=\frac{\pi_{li}}{1-\pi_{li}-\pi_{ri}}h_{li}^{\prime\prime}(u)+\frac{\pi_{ri}}{1-\pi_{li}-\pi_{ri}}h_{ri}^{\prime\prime}(u),$

where the dependence on $X_{i}$ has been suppressed in notations for brevity.
By Lemma C.1 , $b_{i}^{\prime\prime}$ is always positive, hence there exists a
universal positive number $c>0$ such that

$b_{i}^{\prime\prime}(u)>c$

for all $i$ and all $u\in[u_{0}/2,1-u_{0}/2]$, considering that $\Theta$, $B$
and $[u_{0}/2,1-u_{0}/2]$ are all compact. Now for each $i$ consider the
quadratic function

$f_{i}(u)=\frac{c(u-\underline{u}_{i})^{2}}{2}+b_{i}(\underline{u}_{i})$

defined on $[u_{0}/2,1-u_{0}/2]$. Then on the interval $[u_{0}/2,1-u_{0}/2]$ ,
$b_{i}\geq f_{i}$ since $g_{i}=b_{i}-f_{i}$ is strictly convex with
$g_{i}^{\prime}(\underline{u}_{i})=g_{i}(\underline{u}_{i})=0$. Then

$\lambda\left(b_{i}^{-1}([m_{i},m_{i}+\epsilon))\right)\leq\lambda\left(f_{i}^{-1}([m_{i},m_{i}+\epsilon))\right)=\sqrt{\frac{8\epsilon}{c}}$

where the right hand side obviously converges to zero as $\epsilon\rightarrow
0$.

∎

We will now state two crucial “event inclusion” lemmas that involve the most
delicate proofs in this paper, and may be skipped at first reading. To state
them, we conveniently define the long vectors
$\Pi=\Pi(\theta)=(\pi_{li},\pi_{ri})_{i=1}^{m}$ and
$K=K(\beta)=(k_{li},k_{ri})_{i=1}^{m}$ with $2m$ components. Note that they
implicitly depend on the unspecified parameters $\\{\theta,\beta\\}$. As such,
we can also define $\Pi^{*}=\Pi(\theta^{*})$ and $K^{*}=K(\beta^{*})$ to be
the versions evaluated at $\theta^{*}$ and $\beta^{*}$.

######  Lemma D.2 (First event inclusion lemma).

Suppose Assumptions 1 and 2$(i)$ are true and let
$(\theta,\beta)\in\Theta\times B$. For given $\underline{t}>0$ and
$\epsilon>0$, there exists a $\delta>0$ such that whenever
$\|\Pi-\Pi^{*}\|_{\infty}\vee\|K-K^{*}\|_{\infty}<\delta$,

$\left\\{S_{i}^{*}\leq c_{i}^{*}(t)-\epsilon\right\\}\subset\left\\{S_{i}\leq
c_{i}(t)\right\\}\subset\left\\{S_{i}^{*}\leq
c_{i}^{*}(t)+\epsilon\right\\}\text{ for all }i\text{ and all
}t>\underline{t}.$

###### Proof of Lemma D.2.

By Definition C.1, we will show, equivalently, that there exists $\delta>0$
such that whenever $\|\Pi-\Pi^{*}\|_{\infty}\vee\|K-K^{*}\|_{\infty}<\delta$,

(D.6)
$\\{U_{i}\in\underbrace{(\omega^{*}_{iL}(t_{i,+\epsilon}),\omega^{*}_{iR}(t_{i,+\epsilon}))}_{\text{length}=1-c_{i}^{*}(t)-\epsilon}\\}\subset\\{U_{i}\in\underbrace{(\omega_{iL}(t),\omega_{iR}(t))}_{\text{length}=1-c_{i}(t)}\\}\subset\\{U_{i}\in\underbrace{(\omega^{*}_{iL}(t_{i,-\epsilon}),\omega^{*}_{iR}(t_{i,-\epsilon}))}_{\text{length}=1-c_{i}^{*}(t)+\epsilon}\\},$

where we define $t_{i,-\epsilon}:={c_{i}^{*}}^{-1}(c_{i}^{*}(t)-\epsilon)$ and
$t_{i,+\epsilon}:={c_{i}^{*}}^{-1}(c_{i}^{*}(t)+\epsilon)$ that are
respectively less and greater than $t$. In particular, we will first focus on
showing the second inclusion in in (D.6) , which amounts to showing

(D.7) $b_{i}(\omega_{iL}^{*}(t_{i,-\epsilon}))\wedge
b_{i}(\omega_{iR}^{*}(t_{i,-\epsilon}))\geq 1/t$

whenever $\|\Pi-\Pi^{*}\|_{\infty}\vee\|K-K^{*}\|_{\infty}<\delta$, in light
of the fact that

(D.8)
$b_{i}^{*}(\omega_{iL}^{*}(t_{i,-\epsilon}))=b_{i}^{*}(\omega_{iR}^{*}(t_{i,-\epsilon}))=1/t_{i,-\epsilon},$

by the definition of $\omega^{*}_{iL}(\cdot)$ and $\omega^{*}_{iR}(\cdot)$ in
Definition C.1 and properties of $b_{i}$ from Lemma C.2.

Let $u_{0}>0$ be a small positive number such that
$(\omega^{*}_{iL}(\underline{t}),\omega^{*}_{iR}(\underline{t}))\subset[u_{0},1-u_{0}]$
by Lemma D.1(i), and consider the even larger compact interval
$[u_{0}/2,1-u_{0}/2]$. Consider each $b_{i}(u)=b(u;\pi_{i},k_{i})$ as a
function in $(u,\pi_{i},k_{i})$, and let

$\nabla_{\pi,k}b_{i}(u)\equiv\left({\frac{\partial}{\partial\pi_{i}}b(u;\pi_{i},k_{i})}^{T},{\frac{\partial}{\partial
k_{i}}b(u;\pi_{i},k_{i})}^{T}\right)^{T}$

be the gradient of $b_{i}$ with respect to $(\pi_{i},k_{i})$ evaluated at $u$.
Using the compactness of $[u_{0}/2,1-u_{0}/2]\times\Theta\times\mathcal{B}$
and Assumption 2$(i)$ again, one can find a universal constant $C(u_{0})>0$
such that the gradient bounds

(D.9) $\|\nabla_{\pi,k}b_{i}(u)\|_{1}<C(u_{0})\text{ for all }i,\text{ all
}u\in\left[\frac{u_{0}}{2},1-\frac{u_{0}}{2}\right],\text{ all
}(\theta,\beta)\in\Theta\times\boldsymbol{B}.$

On the other hand, without loss of generality, we will let

(D.10) $\epsilon<u_{0}/2$

and, with Lemma D.1$(iv)$, take $\tilde{\epsilon}>0$ be a small enough
constant such that

(D.11)
$\lambda({b^{*}_{i}}^{-1}([b^{*}_{i}(\underline{u}_{i}^{*}),b^{*}_{i}(\underline{u}_{i}^{*})+\tilde{\epsilon}))<\epsilon\text{
for all }i.$

By the mean-value theorem and the gradient bound (D.9), one can then find
$\delta>0$ such that when
$\|\Pi-\Pi^{*}\|_{\infty}\vee\|K-K^{*}\|_{\infty}<\delta$

(D.12) $|b_{i}(u)-b_{i}^{*}(u)|<\tilde{\epsilon}\text{ for all
}u\in\left[\frac{u_{0}}{2},1-\frac{u_{0}}{2}\right].$

By the construction of $\tilde{\epsilon}$ in (D.11) and convexity properties
from Lemma C.2, one must have for all $i$

$b_{i}^{*}(\omega_{iL}^{*}(t_{i,-\epsilon}))-b_{i}^{*}(\omega_{iL}^{*}(t))=b_{i}^{*}(\omega_{iR}^{*}(t_{i,-\epsilon}))-b_{i}^{*}(\omega_{iR}^{*}(t))\geq\tilde{\epsilon},$

which implies

(D.13) $1/t_{i,-\epsilon}\geq\tilde{\epsilon}+1/t.$

by the last property in Definition C.1. Since
$\\{S_{i}^{*}>c^{*}_{i}(t)-\epsilon\\}\supset\\{S_{i}^{*}>c^{*}_{i}(t)\\}$,
from the property (C.5) in Definition C.2 both

$\omega_{iL}^{*}(t_{i,-\epsilon})\in(\omega_{iL}^{*}(t)-\epsilon,\omega_{iL}^{*}(t))\text{
and
}\omega_{iR}^{*}(t_{i,-\epsilon})\in(\omega_{iR}^{*}(t),\omega_{iR}^{*}(t)+\epsilon)$

are true, which implies

$[\omega_{iL}^{*}(t_{i,-\epsilon}),\omega_{iR}^{*}(t_{i,-\epsilon})]\subset[u_{0}/2,1-u_{0}/2],$

considering (D.10) and
$\omega_{iL}^{*}(t),\omega_{iR}^{*}(t)\in[u_{0},1-u_{0}]$ (as
$t\geq\underline{t}$). Therefore by (D.12), we must have

(D.14) $b_{i}(\omega_{iL}^{*}(t_{i,-\epsilon}))\wedge
b_{i}(\omega_{iR}^{*}(t_{i,-\epsilon}))\geq
1/t_{i,-\epsilon}-\tilde{\epsilon}$

given (D.8). Combining (D.13) and (D.14) gives (D.7).

The proof for the first inclusion in (D.6) follows an analogous argument but
is with less resistance, since
$[\omega_{iL}^{*}(t_{i,+\epsilon}),\omega_{iR}^{*}(t_{i,+\epsilon})]\subset[u_{0},1-u_{0}]$
for all $i$. We leave it to the reader.

∎

######  Lemma D.3 (Second event inclusion lemma).

Suppose Assumptions 1 and 2$(i)$ are true and let
$(\theta,\beta)\in\Theta\times B$. For any fixed
$\bar{t}<\min_{i}{c_{i}^{*}}^{-1}(\bar{s})$ with $\bar{s}<1$ and any
${\epsilon}>0$, there exists a $\delta>0$ such that whenever
$\|\Pi-\Pi^{*}\|_{\infty}\vee\|K-K^{*}\|_{\infty}<\delta$,

$\\{S_{i}^{*}>1-c^{*}_{i}(t)-\epsilon\\}\supset\\{S_{i}>1-c_{i}(t)\\}\supset\\{S_{i}^{*}>1-c^{*}_{i}(t)+\epsilon\\},$

for all $i$ and all $t\leq\bar{t}$.

###### Proof of Lemma D.3.

Note that from (C.5) in Definition C.2 and Lemma D.2 we can conclude there
exists a $\delta_{1}>0$ such that whenever
$\|\Pi-\Pi^{*}\|_{\infty}\vee\|K-K^{*}\|_{\infty}<\delta_{1}$,

(D.15) $1-c^{*}_{i}(t)-\epsilon\leq 1-c_{i}(t)\leq 1-c^{*}_{i}(t)+\epsilon.$

Based on (D.15), it suffices to show that there exists a $\delta_{2}>0$ such
that whenever $\|\Pi-\Pi^{*}\|_{\infty}\vee\|K-K^{*}\|_{\infty}<\delta_{2}$,

(D.16)
$\\{S_{i}>1-c^{*}_{i}(t)-\epsilon\\}\subset\\{S_{i}^{*}>1-c^{*}_{i}(t)-2\epsilon\\},$

and

(D.17)
$\\{S_{i}>1-c^{*}_{i}(t)+\epsilon\\}\supset\\{S_{i}^{*}>1-c^{*}_{i}(t)+2\epsilon\\}$

which conclude the lemma by taking $\delta=\delta_{1}\vee\delta_{2}$ and
replacing $\epsilon$ with $\epsilon/2$. In fact, since $t\leq\bar{t}$ and
$c_{i}^{*}(\bar{t})$ are bounded away from 1, we will show the more general
statement: For a given $\underline{s}>0$, there exists
$\delta=\delta(\underline{s})>0$ such that whenever
$\|\Pi-\Pi^{*}\|_{\infty}\vee\|K-K^{*}\|_{\infty}<\delta$,

(D.18) $\\{S_{i}>s\\}\subset\\{S_{i}^{*}>s-\epsilon\\}$

and

(D.19) $\\{S_{i}>s\\}\supset\\{S_{i}^{*}>s+\epsilon\\}$

for all $s\geq\underline{s}$ and all $i$. This will necessitate (D.16) and
(D.17) for $t\leq\bar{t}$.

We will first show (D.18) which amounts to

(D.20)
$(\psi_{iL}(s),\psi_{iR}(s))\subset(\psi_{iL}^{*}(s-\epsilon),\psi_{iR}^{*}(s-\epsilon))$

in light of Definition C.2. In particular, it suffices to only consider the
case where $\epsilon<\underline{s}$, since if $s-\epsilon\leq 0$,
$\\{S_{i}^{*}>s-\epsilon\\}=\\{S_{i}^{*}\geq 0\\}$ becomes the whole
underlying probability space which makes (D.20) trivially true. Now for each
$i$, let

$s(i):=\sup\\{s^{\prime}\in(0,1):(\psi_{iL}(s),\psi_{iR}(s))\subset(\psi_{iL}^{*}(s^{\prime}),\psi_{iR}^{*}(s^{\prime}))\\}.$

By Definition C.2 it must be the case that

(D.21) $s(i)=1-\psi_{iR}^{*}(s(i))+\psi_{iL}^{*}(s(i))\leq
1-\psi_{iR}(s)+\psi_{iL}(s)=s,$

and only one of the following possibilities can be true:

  1. (i)

$\psi_{iL}^{*}(s(i))=\psi_{iL}(s)\text{ and
}\psi_{iR}(s)<\psi_{iR}^{*}(s(i))$,

  2. (ii)

$\psi_{iL}^{*}(s(i))<\psi_{iL}(s)\text{ and
}\psi_{iR}(s)=\psi_{iR}^{*}(s(i))$,

  3. (iii)

$\psi_{iL}^{*}(s(i))=\psi_{iL}(s)\text{ and
}\psi_{iR}(s)=\psi_{iR}^{*}(s(i))$.

In light of the monotone properties in Definition C.2, it suffices to show
that

(D.22) $s-s(i)<\epsilon,$

which will then imply (D.20). Obviously, if $(iii)$ is true then (D.22) must
be true in light of (D.21). We will focus on showing (D.22) in the case of
$(i)$ since the proof for the case of $(ii)$ follows a parallel argument.

By Lemma D.1$(ii)$, there exists a $u_{0}=u_{0}(\underline{s})>0$ such that

(D.23) $(\psi_{iL}(s),\psi_{iR}(s))\subset[u_{0},1-u_{0}]\text{ for all
}i\in\mathbb{N}\text{ and }s\geq\underline{s}.$

Consider each $b_{i}(u)=b(u;\pi_{i},k_{i})$ as a function in
$(u,\pi_{i},k_{i})$, and let

$\nabla_{\pi,k}b_{i}(u):=\left({\frac{\partial}{\partial\pi_{i}}b(u;\pi_{i},k_{i})}^{T},{\frac{\partial}{\partial
k_{i}}b(u;\pi_{i},k_{i})}^{T}\right)^{T}$

be the gradient of $b_{i}$ with respect to $(\pi_{i},k_{i})$ evaluated at $u$.
Using the compactness of $[u_{0},1-u_{0}]\times\Theta\times\mathcal{B}$ and
Assumption 2$(i)$ again, one can find a constant $C(u_{0})>0$ such that the
gradient bounds

(D.24) $\|\nabla_{\pi,k}b_{i}(u)\|_{1}<C(u_{0})\text{ for all }i\text{ and for
all }u\in\left[\frac{u_{0}}{2},1-\frac{u_{0}}{2}\right].$

On the other hand, without loss of generality, with Lemma D.1$(iv)$, let
$\tilde{\epsilon}>0$ be a small enough constant such that

(D.25)
$\lambda({b^{*}_{i}}^{-1}([m^{*}_{i},m^{*}_{i}+\tilde{\epsilon}))<\frac{\epsilon}{2}\text{
for all }i.$

By the mean-value theorem and the gradient bound (D.24), one can then find
$\delta>0$ such that when
$\|\Pi-\Pi^{*}\|_{\infty}\vee\|K-K^{*}\|_{\infty}<\delta$

(D.26) $|b_{i}(u)-b_{i}^{*}(u)|<\tilde{\epsilon}\text{ for all
}u\in\left[u_{0},1-u_{0}\right].$

Since $\psi_{iL}^{*}(s(i))=\psi_{iL}(s)\in[u_{0},1-u_{0}]$, (D.26) and the
last property in Definition C.2 suggest that

(D.27)
$b_{i}^{*}\Bigl{(}\psi_{iR}^{*}(s(i))\Bigr{)}=b_{i}^{*}\Bigl{(}\psi_{iL}^{*}(s(i))\Bigr{)}<b_{i}\Bigl{(}\psi_{iL}(s)\Bigr{)}+\tilde{\epsilon}=b_{i}\Bigl{(}\psi_{iR}(s)\Bigr{)}+\tilde{\epsilon}.$

But since $\psi_{iR}(s)$ is also in the interval $[u_{0},1-u_{0}]$, we must
have

(D.28)
$b_{i}^{*}\Bigl{(}\psi_{iR}(s)\Bigr{)}>b_{i}\Bigl{(}\psi_{iR}(s)\Bigr{)}-\tilde{\epsilon}.$

Combining (D.27) and (D.28), we get that
$2\tilde{\epsilon}>b_{i}^{*}\Bigl{(}\psi_{iR}^{*}(s(i))\Bigr{)}-b_{i}^{*}\Bigl{(}\psi_{iR}(s)\Bigr{)}$
which in light of the construction of $\tilde{\epsilon}$ in (D.25) and
convexity properties from Lemma C.2 gives that

$\psi_{iR}^{*}(s(i))-\psi_{iR}(s)<\epsilon,$

which in turn implies (D.22) by the property (C.5) in Definition C.2. The
proof of (D.19) is similar. It amounts to showing

$(\psi_{iL}^{*}(s+\epsilon),\psi_{iR}^{*}(s+\epsilon))\subset(\psi_{iL}(s),\psi_{iR}(s)).$

We will alternatively define

$\tilde{s}(i):=\inf\\{s^{\prime}\in(0,1):(\psi_{iL}^{*}(s^{\prime}),\psi_{iR}^{*}(s^{\prime}))\subset(\psi_{iL}(s),\psi_{iR}(s))\\}.$

then show $\tilde{s}(i)-s<\epsilon$. We leave the details to the reader. ∎

###  D.3. A Glivenko-Cantelli theorem

######  Lemma D.4 (Pre-Glivenko-Cantelli theorem).

Under Assumptions 1-3 , for any $\epsilon>0$ and positive numbers
$0<\underline{t}<\bar{t}<\min_{i}{c_{i}^{*}}^{-1}(\bar{s})$, there exists
$\delta=\delta(\epsilon)>0$ such that, for sufficiently large $m$,

(D.29)
$\displaystyle\sup_{\begin{subarray}{c}\max(\|K-K^{*}\|,\|\Pi-\Pi^{*}\|)<\delta\\\
\underline{t}\leq
t\leq\bar{t}\end{subarray}}\left|\frac{1}{m}\sum_{i=1}^{m}I(T_{i}\leq
t)-G(t)\right|\leq\epsilon$ (D.30)
$\displaystyle\sup_{\begin{subarray}{c}\max(\|K-K^{*}\|,\|\Pi-\Pi^{*}\|)<\delta\\\
\underline{t}\leq
t\leq\bar{t}\end{subarray}}\left|\frac{1}{m}\sum_{i=1}^{m}(1-H_{i})I(T_{i}\leq
t)-\pi_{0}{G}_{0}(t)\right|\leq\epsilon,$ (D.31)
$\displaystyle\sup_{\begin{subarray}{c}\max(\|K-K^{*}\|,\|\Pi-\Pi^{*}\|)<\delta\\\
\underline{t}\leq
t\leq\bar{t}\end{subarray}}\left|\frac{1}{m}\sum_{i=1}^{m}I(S_{i}\geq
1-c_{i}(t))-\bar{G}_{0}(t)\right|\leq\epsilon$

with probability 1.

###### Proof of Lemma D.4.

In this proof, for any function $F(\cdot)$, $F(t-)$ denotes the left limit at
the point $t$.

Proof of (D.29): Let

$G_{\theta,\beta}(t)=\frac{1}{m}\sum_{i=1}^{m}I(T_{i}\leq t),$

where the subscript emphasizes that the $T_{i}$’s are defined with an
unspecified $(\theta,\beta)$, to distinguish from $G$ in Assumption 3. Let $n$
be large enough such that $1/n<\epsilon/2$ and consider
$G^{\leftarrow}(1/n)\leq\dots\leq G^{\leftarrow}(n/n)$. If we define

$n^{\prime}:=\min\left\\{i:\underline{t}<G^{\leftarrow}(i/n)<\bar{t},i=1,\dots,n\right\\},$
$d^{\prime}:=|\left\\{i:\underline{t}<G^{\leftarrow}(i/n)<\bar{t},i=1,\dots,n\right\\}|.$

Define
$t_{1}:=G^{\leftarrow}(n^{\prime}/n),t_{2}:=G^{\leftarrow}((n^{\prime}+1)/n),\dots,t_{d^{\prime}}:=G^{\leftarrow}((n^{\prime}+d^{\prime}-1)/n)$,
as well as $t_{0}=\underline{t}$ and $t_{d}=\bar{t}$ with $d=d^{\prime}+1$.
Following the proof of the Glivenko-Cantelli theorem in Resnick (2019, p.224),
we have

(D.32) $\sup_{\underline{t}\leq
t\leq\bar{t}}\left|G_{\theta,\beta}(t)-G(t)\right|\leq\left(\bigvee_{v=0}^{d}|G_{\theta,\beta}(t_{v})-G(t_{v})|\vee|G_{\theta,\beta}(t_{v}-)-G(t_{v}-)|\right)+1/n$

We will first bound the terms of the form $|G_{\theta,\beta}(t_{v})-G(t_{v})|$
in (D.32). The strong law of large numbers for $G(\cdot)$ in Assumption 3
suggests that

(D.33)
$|G_{\theta,\beta}(t_{v})-G(t_{v})|\leq\left|m^{-1}\sum_{i=1}^{m}(I(T_{i}\leq
t_{v})-I(T_{i}^{*}\leq t_{v}))\right|+R_{v},$

where the remainder term $R_{v}\longrightarrow 0$ almost surely. Now,
realizing $\\{T_{i}^{*}\leq t_{v}\\}=\\{S_{i}^{*}\leq c_{i}^{*}(t_{v})\\}$, by
Lemma D.2 and $\underline{t}>0$, pick $\delta>0$ such that

(D.34) $\left|\frac{1}{m}\sum_{i=1}^{m}(I(T_{i}\leq t_{v})-I(T_{i}^{*}\leq
t_{v}))\right|\leq\\\
\underbrace{\left|m^{-1}\sum_{i=1}^{m}\left[I(S^{*}_{i}\leq
c_{i}^{*}(t_{v})+\epsilon)-I(S_{i}^{*}\leq
c_{i}^{*}(t_{v}))\right]\right|}_{(A)}\vee\underbrace{\left|m^{-1}\sum_{i=1}^{m}\left[I(S^{*}_{i}\leq
c_{i}^{*}(t_{v})-\epsilon)-I(S_{i}^{*}\leq
c_{i}^{*}(t_{v}))\right]\right|}_{(B)}.$

for $\max(\|\Pi-\Pi^{*}\|,\|K-K^{*}\|)<\delta$. This is because
$\sum_{i=1}^{m}\left[I(S^{*}_{i}\leq
c_{i}^{*}(t_{v})+\epsilon)-I(S_{i}^{*}\leq
c_{i}^{*}(t_{v}))\right]\geq\sum_{i=1}^{m}(I(T_{i}\leq t_{v})-I(T_{i}^{*}\leq
t_{v}))$ if the latter term is greater than $0$; likewise,
$\sum_{i=1}^{m}\left[I(S^{*}_{i}\leq
c_{i}^{*}(t_{v})-\epsilon)-I(S_{i}^{*}\leq
c_{i}^{*}(t_{v}))\right]\leq\sum_{i=1}^{m}(I(T_{i}\leq t_{v})-I(T_{i}^{*}\leq
t_{v}))$ if the latter is less than $0$.

We will first develop a bound for term $(A)$. One have

$\displaystyle(A)$
$\displaystyle\leq\mathbb{E}\Bigl{|}\frac{1}{m}\sum_{i=1}^{m}[P(S^{*}_{i}>c_{i}^{*}(t_{v})|X_{i})-P(S^{*}_{i}>c_{i}^{*}(t_{v})+\epsilon|X_{i})]\Bigr{|}+Q_{v}^{+}$
$\displaystyle=\mathbb{E}\Bigl{|}\frac{1}{m}\sum_{i=1}^{m}\bigl{[}P\bigl{(}U_{i}\in(\psi^{*}_{iL}(c_{i}^{*}(t_{v})),\psi^{*}_{iR}(c_{i}^{*}(t_{v}))|X_{i}\bigr{)}-$
$\displaystyle\hskip
85.35826ptP(U_{i}\in(\psi^{*}_{iL}(c_{i}^{*}(t_{v})+\epsilon),\psi^{*}_{iR}(c_{i}^{*}(t_{v})+\epsilon)|X_{i})\bigr{]}\Bigr{|}+Q_{v}^{+}$
(D.35) $\displaystyle\leq C(\underline{t})\epsilon+Q_{v}^{+}$

where $Q_{v}^{+}=o_{a.s.}(1)$ is a remainder term coming from the strong law
of $G(\cdot)$ in Assumption 3. The second equality comes from Definition C.2.
Note that the intervals

$(\psi^{*}_{iL}(c_{i}^{*}(t_{v})),\psi^{*}_{iR}(c_{i}^{*}(t_{v})),i=1,\dots,m$

can be equivalently represented as

$(\omega^{*}_{iL}(t_{v}),\omega^{*}_{iR}(t_{v})),i=1,\dots,m$

by Definition C.1, which all belong to a compact sub-interval in $(0,1)$ by
Lemma D.1$(i)$ and the fact that $\underline{t}>0$. As such, Assumption
2$(ii)$ can be applied to give the last inequality (D.35). By realizing, from
Definitions C.1 and C.2, that the event

$\\{S^{*}_{i}>c_{i}^{*}(t_{v})-\epsilon\\}$

is equivalent to $U_{i}$ belonging to an interval that is
$(\omega^{*}_{iL}(t_{v}),\omega^{*}_{iR}(t_{v}))$ expanded by a further
$\epsilon$ width, it is obvious that one can analogously develop the bound

(D.36) $(B)\leq C(\underline{t})\epsilon+Q_{v}^{-}$

for a constant $C(\underline{t})$ and $Q_{v}^{-}=o_{a.s.}(1)$. Combining
(D.33), (D.35) and (D.36) gives (by appropriately adjusting $\epsilon$)

(D.37) $|G_{\theta,\beta}(t_{v})-G(t_{v})|\leq\epsilon/2\text{ a.s.},$

for sufficiently large $m$. A similar bound

(D.38) $|G_{\theta,\beta}(t_{v}-)-G(t_{v}-)|\leq\epsilon/2\text{ a.s.},$

can be derived in much the same way with no difficulty by first writing

$G_{\theta,\beta}(t_{v}-)-G(t_{v}-)=m^{-1}\sum_{i=1}^{m}I(T_{i}<t_{v})-G(t_{v})$

using the continuity of $G(\cdot)$, and the proof is omitted for brevity.
Combining (D.32), (D.37) and (D.38) give (D.29). (D.30) can be proved in the
same way by first writing

$\sum_{i=1}^{m}(1-H_{i})I(T_{i}\leq t)=\sum_{H_{i}=0}I(T_{i}\leq t),$

and noting that $U_{i}$ is uniformly distributed given $H_{i}=0$ and $X_{i}$,
and thus omitted.

Proving (D.31) is similar. One proceed by developing the bound

$\displaystyle\left|\frac{1}{m}\sum_{i}[I(S_{i}>1-c_{i}(t_{v}))-P(S_{i}^{*}>1-c_{i}^{*}(t_{v}))]\right|$
$\displaystyle\leq\left|\frac{1}{m}\sum_{i}[I(S_{i}^{*}>1-c_{i}^{*}(t_{v})-\epsilon)-P(S_{i}^{*}>1-c_{i}^{*}(t_{v}))]\right|\vee$
$\displaystyle\hskip
113.81102pt\left|\frac{1}{m}\sum_{i}[I(S^{*}_{i}>1-c^{*}_{i}(t_{v})+\epsilon)-P(S_{i}^{*}>1-c_{i}^{*}(t_{v}))]\right|$

with Lemma D.3 for $t_{v}$ that are now quantiles of $\bar{G}_{0}$, for
sufficiently small $\delta>0$ and $\max(\|\Pi-\Pi^{*}\|,\|K-K^{*}\|)<\delta$.
From Definition C.2, the events $\\{S_{i}^{*}>1-c_{i}^{*}(t_{v})\\}$ have the
form

$\\{U_{i}\in(\psi^{*}_{iL}(1-c_{i}^{*}(t_{v})),\psi^{*}_{iR}(1-c_{i}^{*}(t_{v})))\\},$

and to show that the intervals
$(\psi^{*}_{iL}(1-c_{i}^{*}(t_{v})),\psi^{*}_{iR}(1-c_{i}^{*}(t_{v}))$ can be
placed in a compact sub-interval of $(0,1)$ by Lemma D.1$(ii)$ to apply
Assumption 2$(ii)$, one need to show that
$c_{1}^{*}(t_{v}),\dots,c_{m}^{*}(t_{v})$ are bounded away from 1. This is
true because $c_{i}^{*}(t_{v})\leq c_{i}^{*}(\bar{t})\leq\bar{s}<1$ for all
$i$ by definition. The same proof rundown goes through, again, by realizing
that
$(\psi^{*}_{iL}(1-c_{i}^{*}(t_{v})-\epsilon),\psi^{*}_{iR}(1-c_{i}^{*}(t_{v})-\epsilon)$
are just $\epsilon$-expansion of
$(\psi^{*}_{iL}(1-c_{i}^{*}(t_{v})),\psi^{*}_{iR}(1-c_{i}^{*}(t_{v}))$ from
Definition C.2. The rest of the proof goes thru with no resistance.

∎

######  Lemma D.5 (Glivenko-Cantelli Theorems).

Under Assumptions 1-3, for any
$0<\underline{t}<\bar{t}\leq\min_{i}{c_{i}^{*}}^{-1}(\bar{s})$, we have

(D.39) $\displaystyle\sup_{\underline{t}\leq
t\leq\bar{t}}\left|\frac{1}{m}\sum_{i=1}^{m}I(\hat{T}_{i}\leq
t)-G(t)\right|\longrightarrow 0,$ (D.40) $\displaystyle\sup_{\underline{t}\leq
t\leq\bar{t}}\left|\frac{1}{m}\sum_{i=1}^{m}(1-H_{i})I(\hat{T}_{i}\leq
t)-\pi_{0}G_{0}(t)\right|\longrightarrow 0,$ (D.41)
$\displaystyle\sup_{\underline{t}\leq
t\leq\bar{t}}\left|\frac{1}{m}\sum_{i=1}^{m}I(\hat{S}_{i}\geq
1-\hat{c}_{i}(t))-\bar{G}_{0}(t)\right|\longrightarrow 0,$

almost surely.

###### Proof of Lemma D.5.

Let
$\hat{\theta}=(\hat{\theta}_{l}^{T},\hat{\theta}_{r}^{T})^{T},\hat{\beta}=(\hat{\beta}_{l}^{T},\hat{\beta}_{r}^{T})^{T}$,
and define $\hat{\Pi}=\Pi(\hat{\theta})$ and $\hat{K}=K(\hat{\beta})$. We will
first show that

(D.42)
$\|\hat{\Pi}-\Pi^{*}\|_{\infty}\vee\|\hat{K}-K^{*}\|_{\infty}\longrightarrow
0\text{ a.s.},$

which is a consequence of

(D.43) $\hat{\theta}\longrightarrow\theta^{*}\text{ and
}\hat{\beta}\longrightarrow\beta^{*}\ \ a.s.$

by the mean value theorem, the compactness of $\Theta\times\boldsymbol{B}$ and
Assumption 2$(i)$. To show (D.43), it suffices to bound

(D.44)
$\left|\log[(1-\pi_{li}-\pi_{ri})+\pi_{li}h_{li}(U_{i})+\pi_{ri}h_{ri}(U_{i})]\right|=|\log
h_{X_{i}}(U_{i})|$

by an integrable function in $U_{i}$ that doesn’t depend on $(\theta,\beta)$
(White, 1981, Theorem 2.1). We first let $u_{0}>0$ be as in Lemma D.1$(iii)$.
By the compactness of $\Theta\times\boldsymbol{B}\times[u_{0},1-u_{0}]$ and
Assumption 2$(i)$, there exists universal constants $C,c>0$ such that
$\log(c)<0$, $\log(C)>0$, and

$C\geq h_{X_{i}}(U_{i})\geq c$

for all $(\theta,\beta,X_{i})$ whenever $U_{i}\in[u_{0},1-u_{0}]$. Note that
Lemma C.2 also implies that

$h_{X_{i}}(U_{i})\geq c\text{ for }U_{i}\in(0,u_{0})\cup(1-u_{0},1).$

Moreover, for $U_{i}\in(0,u_{0})$, $\pi_{ri}h_{ri}(U_{i})\leq C$ by Lemma C.1,