11footnotetext: Department of Statistics, Florida State University,
Tallahassee, FL 32306, USA.22footnotetext: School of Mathematics, Jilin
University, Changchun, Jilin 130012, China.

# A powerful empirical Bayes approach for high dimensional replicability
analysis

Haochen Lei1, Yan Li2 and Hongyuan Cao1

###### Abstract

Identifying replicable signals across different studies provides stronger
scientific evidence and more powerful inference. Existing literature on high
dimensional replicability analysis either imposes strong modeling assumptions
or has low power. We develop a powerful and robust empirical Bayes approach
for high dimensional replicability analysis. Our method effectively borrows
information from different features and studies while accounting for
heterogeneity. We show that the proposed method has better power than
competing methods while controlling the false discovery rate, both empirically
and theoretically. Analyzing datasets from the genome-wide association studies
reveals new biological insights that otherwise cannot be obtained by using
existing methods.

Keywords: False discovery rate; High dimension; Local false discovery rate;
Replicability analysis

##  1 Introduction

Replicability is a fundamental principle in scientific research. Consistent
findings among different times, places, and populations provide stronger
scientific evidence. Replicability attracts attention, scrutiny, and debate in
the popular press and the scientific community. An analysis of past studies
indicates that the cumulative prevalence of irreplicable preclinical research
exceeds $50\%$ (Ioannidis, 2005). Approximately $\$28$ billion per year is
spent on preclinical research that is not replicable in the United States
(Freedman et al., 2015). There are a myriad of reasons for irreplicable
results, such as study design, experimental protocols, and data analysis. We
focus on statistical analysis of data from high-throughput experiments. In
particular, we study conceptual replicability, where consistent results are
obtained using different populations that target the same scientific question.
Such replicable findings rule out spurious signals that unmeasured confounders
may cause in a single study, thus providing more robust scientific evidence.

In high-throughput experiments, many candidates are evaluated for their
association with a biological feature of interest, and those with significant
associations are identified for further analysis. We aim to simultaneously
identify replicable features from high-throughput experiments in multiple
studies. To identify significant features of a single study, a commonly used
metric is the false discovery rate (FDR) proposed by Benjamini and Hochberg
(1995) (BH), which is defined as the ratio of false rejections to total
rejections. An ad hoc approach to replicability analysis is to apply the BH
procedure separately in each study and then declare the overlapped results as
replicable findings. This approach does not control FDR and, moreover, has low
power as it does not borrow information from different studies (Bogomolov and
Heller, 2022). To have a valid FDR control, the MaxP method is proposed in
Benjamini et al. (2009), where the maximum of $p$-values in multiple studies
is used as the test statistic for each feature, and the BH procedure is
implemented afterward. This is overly conservative, as the composite null
feature of replicability is not considered. To improve power, Lyu et al.
(2023) enumerates the composition null by incorporating the proportion
estimation, but does not develop any theoretical results. To borrow
information from different studies, Bogomolov and Heller (2018) pre-screens
features in each study and develops a cross screening method to identify
replicable signals. Motivated by the agreement of top-ranked signals, Li et
al. (2011) proposes a reproducible discovery rate and graphical tools to
assess the replicability and control the error rate. A non-parametric
procedure based on the ranking of different features is proposed in Philtron
et al. (2018). Both Li et al. (2011) and Philtron et al. (2018) ignore the
heterogeneity of different studies and assume that a feature is either null or
non-null in all studies. An empirical Bayes approach for replicability
analysis is proposed in Heller and Yekutieli (2014), where the test statistic
is transformed into a $z$-score and parametric modeling is used. From a
Bayesian perspective, Zhao et al. (2020) proposes probabilistic models to
evaluate and control replicability for high-throughput experiments in which
the tuning parameters must be selected. The criteria for the evaluation of
replicability can be found in Hung and Fithian (2020) and Bogomolov and Heller
(2013) examine replicability in a two-stage analysis with a primary study and
a follow-up study. Replicability analysis shares the composite null structure
with mediation analysis but has a distinctive feature of information borrowing
across different studies. Related work on high dimensional mediation analysis
can be found in Liu et al. (2022) and James Y. Dai and LeBlanc (2022), among
others. A recent survey on replicability analysis can be found in Bogomolov
and Heller (2022).

In this paper, we develop a robust and powerful empirical Bayes approach for
high dimensional replicability analysis. We assume that the data are
summarized in $p$-values for each study. We use $p$-values mainly for
versatility. Without loss of generality, we use two studies to illustrate. To
account for the heterogeneity of different studies, we use four different
states to model hidden states of $p$-values from two studies. The composite
null of replicability analysis comprises three different states: zero effect
in both studies, zero effect in the first study, and non-zero effect in the
second study and vice versa. The empirical Bayes approach allows us to
enumerate different states of the composite null. Conditional on hidden
states, the distribution of paired $p$-values is modeled by a four-group
mixture model (Efron, 2012). We allow the density functions of $p$-values
under non-null to vary across different studies. Furthermore, instead of the
predominant parametric modeling of the $p$-value density function under the
non-null, we use a non-parametric density function estimation under the shape
constraint (Grenander, 1956). This has the flexibility of non-parametric
modeling and the convenience of no tuning parameters. The local false
discovery rate (Lfdr), defined as the posterior probability of being
replicability null, is used as a test statistic. We combine the EM algorithm
(Dempster et al., 1977) and the pool-adjacent-violator-algoirthm (PAVA)
(Robertson et al., 1988) to estimate unknown parameters and unknown functions
in the calculation of Lfdr. With tools from empirical processes (Van de Geer,
1993, 2000), we establish the convergence rate of the estimation through the
Hellinger distance. The proposed method can achieve asymptotic FDR control.
Extensive simulation studies support our theoretical predictions. By applying
our procedure to genome-wide association studies (GWAS) of type 2 diabetes
(T2D) of European and East Asian ancestry, we identify replicable single
nucleotide polymorphisms (SNPs) that could be missed by using existing
methods.

The rest of the paper proceeds as follows. Section 2 describes the proposed
empirical Bayesian method for replicability analysis, where estimates of
unknown parameters and unknown functions are obtained. Section 3 investigates
theoretical properties of the proposed estimators, and establishes the
asymptotic FDR control of our procedure based on estimated Lfdr. Sections 4
and 5 are devoted to simulation studies and data analysis, respectively.
Concluding remarks are given in Section 6.

##  2 Methodology

###  2.1 Notation and problem set up

Suppose we are interested in testing $m$ hypotheses simultaneously. We have
data from different studies summarized in the $p$-values for each hypothesis.
Without loss of generality, we use two studies to illustrate. Let
$\theta_{ji},j=1,2;i=1,\ldots,m$ denote the true but unobserved states in
study $j$ and hypothesis $i$. Conditional on the hidden states, assume the
p-values follow a mixture model (Efron, 2012):

$\displaystyle
p_{1i}|\theta_{1i}\sim\left(1-\theta_{1i}\right)f_{0}+\theta_{1i}f_{1},i=1,\ldots,m,$
(1) $\displaystyle
p_{2i}|\theta_{2i}\sim\left(1-\theta_{2i}\right)f_{0}+\theta_{2i}f_{2},i=1,\ldots,m,$

where $f_{0}$ is the density function of $p$-values under the null and $f_{1}$
and $f_{2}$ are the density functions of $p$-values under the non-null from
study 1 and study 2, respectively. This approach accounts for the
heterogeneity of hidden states ($\theta_{1i}$ and $\theta_{2i}$) and effect
sizes ($f_{1}$ and $f_{2}$) of two studies.

For $i=1,\ldots,m,$ denote

$P\left(\tau_{i}=(u,v)\right)=\xi_{u,v},\quad u,v=0,1,$

where

$\tau_{i}=\begin{cases}(0,0)&\text{ if }\theta_{1i}=\theta_{2i}=0,\\\
(0,1)&\text{ if }\theta_{1i}=0,\theta_{2i}=1,\\\ (1,0)&\text{ if
}\theta_{1i}=1,\theta_{2i}=0,\\\ (1,1)&\text{ if
}\theta_{1i}=\theta_{2i}=1,\quad i=1,\ldots,m.\end{cases}$

The null hypothesis for replicability is

$H_{0i}:\theta_{1i}\theta_{2i}=0,i=1,\cdots,m.$ (2)

This is a composite null hypothesis that consists of three different states. A
feature is replicable only if the true but unobserved states of both studies
take the value of $1.$ This effectively provides stronger evidence as a
feature in a confounded study may be falsely identified as a signal but will
not be picked up as replicable with the null specified in (2). Define the
local false discovery rate as the posterior probability of being null given
data (Efron, 2012). Specifically, we have

$\displaystyle\operatorname{Lfdr}(p_{1i}$
$\displaystyle,p_{2i}):=P\left(\theta_{1i}\theta_{2i}=0\mid
p_{1i},p_{2i}\right)$ (3)
$\displaystyle=1-P\left(\theta_{1i}=\theta_{2i}=1\mid p_{1i},p_{2i}\right)$
$\displaystyle=\frac{\xi_{00}f_{0}\left(p_{1i}\right)f_{0}\left(p_{2i}\right)+\xi_{01}f_{0}\left(p_{1i}\right)f_{2}\left(p_{2i}\right)+\xi_{10}f_{1}\left(p_{1i}\right)f_{0}\left(p_{2i}\right)}{\xi_{00}f_{0}\left(p_{1i}\right)f_{0}\left(p_{2i}\right)+\xi_{01}f_{0}\left(p_{1i}\right)f_{2}\left(p_{2i}\right)+\xi_{10}f_{1}\left(p_{1i}\right)f_{0}\left(p_{2i}\right)+\xi_{11}f_{1}\left(p_{1i}\right)f_{2}\left(p_{2i}\right)}.$

In (3), the three components of the composite null are enumerated as the four-
component mixture model treats each component equally. It is usually assumed
that $f_{0}$ follows a standard uniform distribution. For $f_{1}$ and $f_{2},$
non-parametric estimations are desirable for robust inference. We impose the
monotonicity condition (Sun and Cai, 2007; Cao et al., 2013) and assume that

$\displaystyle f_{1}/f_{0}\mbox{ and }f_{2}/f_{0}\mbox{ are non-increasing}.$
(4)

Under (4), for $x_{1}\leq x_{2}$ and $y_{1}\leq y_{2},$ we have
$\mbox{Lfdr}(x_{1},y_{1})\leq\mbox{Lfdr}(x_{2},y_{2}).$

We use the $\operatorname{Lfdr}$ defined in (3) as the test statistic for the
replicability null specified in (2). For the $i$th hypothesis, we can write
the rejection rule as

$I\left\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\right\\}.$ (5)

The total number of rejections is

$R_{m}(\lambda)=\sum\limits_{i=1}^{m}I\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\\}.$

In the rejections $R_{m}$, the number of false rejections is

$V_{m}(\lambda)=\sum\limits_{i=1}^{m}I\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\\}(1-\theta_{1i}\theta_{2i}).$

We aim to control the false discovery rate (FDR) at level $\alpha$, where

$\mbox{FDR}(\lambda)=E\left\\{\frac{V_{m}(\lambda)}{1\vee
R_{m}(\lambda)}\right\\}=E\left\\{\mbox{FDP}_{m}(\lambda)\right\\},$ (6)

where the false discovery proportion is defined as

$\mbox{FDP}_{m}(\lambda)=\frac{\sum_{i=1}^{m}I\left\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\right\\}(1-\theta_{1i}\theta_{2i})}{\sum_{i=1}^{m}I\left\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\right\\}}.$

To calculate the Lfdr and find critical value $\hat{\lambda}_{m}$, we estimate
the unknown parameters $\xi_{00},\xi_{01},\xi_{10},\xi_{11}$ and unknown
functions $f_{1},$ and $f_{2}$ by the maximum likelihood estimation.

###  2.2 Estimation

Assume that $f_{0}=1\\{0\leq x\leq 1\\}.$ Define unknown parameters and
unknown functions as

$w=(\xi_{00},\xi_{01},\xi_{10},\xi_{11},f_{1},f_{2}),$

which falls in the space $\Omega,$ where $\Omega$ is a semi-parametric space

$\Omega=\Delta\times\mathcal{G}_{1}\times\mathcal{G}_{2},$ (7)

where

$\Delta=\\{(\xi_{00},\xi_{10},\xi_{01},\xi_{11})\in
R^{4}_{+}:\xi_{00}+\xi_{10}+\xi_{01}+\xi_{11}=1\\}$
$\mathcal{G}_{1}=\\{f_{1}:f_{1}\in\mathcal{G}\\}\quad\text{and}\quad\mathcal{G}_{2}=\\{f_{2}:f_{2}\in\mathcal{G}\\},$
(8)

where

$\mathcal{G}=\\{f:\text{ non-increasing density function supported on
}(0,1)\\}.$

The probability density function of $p$-values from two studies can be written
as

$p_{w}(p_{1},p_{2})=\xi_{00}+\xi_{10}f_{1}(p_{1})+\xi_{01}f_{2}(p_{2})+\xi_{11}f_{1}(p_{1})f_{2}(p_{2}).$

The maximum likelihood estimator is

$\hat{w}_{m}=\operatornamewithlimits{argmax}_{w\in\Omega}\frac{1}{m}\sum\limits_{i=1}^{m}\log
p_{w}(p_{1i},p_{2i}).$ (9)

We use the EM algorithm (Dempster et al., 1977) and the pool-adjacent-violator
algorithm (PAVA)(de Leeuw et al., 2009) to solve the optimization problem (9).
Specifically, let $\mathbf{p}_{1}=\\{p_{1i}\\}_{i=1}^{m}$ and
$\mathbf{p}_{2}=\\{p_{2i}\\}_{i=1}^{m}$ be the observed $p$-values, and
$\bm{\theta}_{1}=\\{\theta_{1i}\\}_{i=1}^{m}$ and
$\bm{\theta}_{2}=\\{\theta_{2i}\\}_{i=1}^{m}$ be corresponding hidden states.
The joint log-likelihood function of
($\mathbf{p}_{1},\mathbf{p}_{2},\bm{\theta}_{1},\bm{\theta}_{2}$) can be
written as

$\displaystyle
l(\mathbf{p}_{1},\mathbf{p}_{2},\bm{\theta}_{1},\bm{\theta}_{2})$
$\displaystyle=\frac{1}{m}\sum\limits_{i=1}^{m}\\{(1-\theta_{1i})(1-\theta_{2i})\log\xi_{00}+\theta_{1i}(1-\theta_{2i})\log\xi_{10}+(1-\theta_{1i})\theta_{2i}\log\xi_{01}$
$\displaystyle\quad+\theta_{1i}\theta_{2i}\log\xi_{11}+\theta_{1i}\log
f_{1}(p_{1i})+\theta_{2i}\log f_{2}(p_{2i})\\}.$

We iteratively implement the following two steps.

E-step: Given current estimates of
$w^{(k)}=(\xi_{00}^{(k)},\xi_{10}^{(k)},\xi_{01}^{(k)},\xi_{11}^{(k)},f_{1}^{(k)},f_{2}^{(k)})$,
obtain the conditional expectation of the log-likelihood function as

$\displaystyle Q(w|w^{(k)})$
$\displaystyle=E_{\bm{\theta}_{1},\bm{\theta}_{2}\mid
w^{(k)}}\\{l(\mathbf{p}_{1},\mathbf{p}_{2},\bm{\theta}_{1},\bm{\theta}_{2})\\}$
$\displaystyle=\frac{1}{m}\sum\limits_{i=1}^{m}\\{\gamma_{i,00}^{(k)}\log\xi_{00}+\gamma_{i,10}^{(k)}\log\xi_{10}+\gamma_{i,01}^{(k)}\log\xi_{01}+\gamma_{i,11}^{(k)}\log\xi_{11}$
$\displaystyle\quad+(\gamma_{i,10}^{(k)}+\gamma_{i,11}^{(k)})\log
f_{1}(p_{1i})+(\gamma_{i,01}^{(k)}+\gamma_{i,11}^{(k)})\log f_{2}(p_{2i})\\},$

where the posterior probabilities of hidden states are calculated as follows.

$\displaystyle\gamma_{i,00}^{(k)}$
$\displaystyle=\frac{\xi_{00}^{(k)}}{\xi_{00}^{(k)}+\xi_{10}^{(k)}f_{1}^{(k)}(p_{1i})+\xi_{01}^{(k)}f_{2}^{(k)}(p_{2i})+\xi_{11}^{(k)}f_{1}^{(k)}(p_{1i})f_{2}^{(k)}(p_{2i})},$
$\displaystyle\gamma_{i,10}^{(k)}$
$\displaystyle=\frac{\xi_{10}^{(k)}f_{1}^{(k)}(p_{1i})}{\xi_{00}^{(k)}+\xi_{10}^{(k)}f_{1}^{(k)}(p_{1i})+\xi_{01}^{(k)}f_{2}^{(k)}(p_{2i})+\xi_{11}^{(k)}f_{1}^{(k)}(p_{1i})f_{2}^{(k)}(p_{2i})},$
$\displaystyle\gamma_{i,01}^{(k)}$
$\displaystyle=\frac{\xi_{01}^{(k)}f_{2}^{(k)}(p_{2i})}{\xi_{00}^{(k)}+\xi_{10}^{(k)}f_{1}^{(k)}(p_{1i})+\xi_{01}^{(k)}f_{2}^{(k)}(p_{2i})+\xi_{11}^{(k)}f_{1}^{(k)}(p_{1i})f_{2}^{(k)}(p_{2i})},$
$\displaystyle\gamma_{i,11}^{(k)}$
$\displaystyle=\frac{\xi_{11}^{(k)}f_{1}^{(k)}(p_{1i})f_{2}^{(k)}(p_{2i})}{\xi_{00}^{(k)}+\xi_{10}^{(k)}f_{1}^{(k)}(p_{1i})+\xi_{01}^{(k)}f_{2}^{(k)}(p_{2i})+\xi_{11}^{(k)}f_{1}^{(k)}(p_{1i})f_{2}^{(k)}(p_{2i})}.$

M-step: Update $w^{(k+1)}$ by

$w^{(k+1)}=\operatornamewithlimits{argmax}_{w\in\Omega}Q(w|w^{(k)}).$

subject to the constraint that $\xi_{00}+\xi_{01}+\xi_{10}+\xi_{11}=1$ and
$f_{1}$ and $f_{2}$ are non-increasing density functions. We have

$\xi_{00}^{(k+1)}=\frac{1}{m}\sum\limits_{i=1}^{m}\gamma_{i,00}^{(k)},\quad\xi_{10}^{(k+1)}=\frac{1}{m}\sum\limits_{i=1}^{m}\gamma_{i,10}^{(k)},\quad\xi_{01}^{(k+1)}=\frac{1}{m}\sum\limits_{i=1}^{m}\gamma_{i,01}^{(k)},\quad\xi_{11}^{(k+1)}=\frac{1}{m}\sum\limits_{i=1}^{m}\gamma_{i,11}^{(k)},$

and

$\displaystyle f_{1}^{(k+1)}$
$\displaystyle=\operatornamewithlimits{argmax}_{f_{1}\in\mathbb{H}}\sum\limits_{i=1}^{m}(\gamma_{i,10}^{(k)}+\gamma_{i,11}^{(k)})\log
f_{1}(p_{1i}),\quad$ (10) $\displaystyle f_{2}^{(k+1)}$
$\displaystyle=\operatornamewithlimits{argmax}_{f_{2}\in\mathbb{H}}\sum\limits_{i=1}^{m}(\gamma_{i,01}^{(k)}+\gamma_{i,11}^{(k)})\log
f_{2}(p_{2i}),$ (11)

where $\mathbb{H}$ is the space of a non-increasing density function. We
repeat the above E-step and M-step until the algorithm converges.

Next, we provide specific steps to solve (10). Let $0=p_{1,(0)}\leq\cdots\leq
p_{1,(i)}\leq\cdots\leq p_{1,(m)}\leq 1=p_{1,(m+1)}$ be the order statistics
of $\mathbf{p_{1}}.$ Denote
$t_{1i}=\gamma_{i,10}^{(k)}+\gamma_{i,11}^{(k)},i=1,\dots,m$ and let
$t_{1,(i)}$ corresponds to $p_{1,(i)}$. Denote $y_{i}=f_{1}(p_{1,(i)})$,
$\mathbf{y}=(y_{1},\cdots,y_{m})$ and
$\mathcal{Q}=\\{(y_{1},\cdots,y_{m}):y_{1}\geq y_{2}\geq\cdots\geq y_{m}\\}.$
The optimization problem (10) can be written as

$\hat{\mathbf{y}}=\operatornamewithlimits{argmax}_{\mathbf{y}\in\mathcal{Q}}\sum\limits_{i=1}^{m}t_{1,(i)}\log
y_{{i}},\text{ subject to
}\sum\limits_{i=1}^{m}\\{p_{1,(i)}-p_{1,(i-1)}\\}y_{i}=1.$

We first ignore the monotonic constraint $\mathcal{Q}$. By applying the
Lagrangian multiplier, the objective function to maximize is

$L(\mathbf{y},\lambda)=\sum\limits_{i=1}^{m}t_{1,(i)}\log
y_{i}+\lambda\left[\sum\limits_{i=1}^{m}\\{p_{1,(i)}-p_{1,(i-1)}\\}y_{i}-1\right].$

The maximizer of $\lambda$ and $y_{i}$ satisfy

$\displaystyle\frac{\partial L(\mathbf{y},\lambda)}{\partial y_{i}}$
$\displaystyle=\frac{t_{1,(i)}}{y_{i}}+\lambda\\{p_{1,(i)}-p_{1,(i-1)}\\}=0,$
$\displaystyle\frac{\partial L(\mathbf{y},\lambda)}{\partial\lambda}$
$\displaystyle=\sum\limits_{i=1}^{m}\\{p_{1,(i)}-p_{1,(i-1)}\\}y_{i}-1=0.$

We have

$\displaystyle\tilde{\lambda}=-\sum_{i=1}^{m}t_{1,(i)},$
$\displaystyle\tilde{y}_{i}=\frac{t_{1,(i)}}{\\{\sum_{i=1}^{m}t_{1,(i)}\\}\\{p_{1,(i)}-p_{1,(i-1)}\\}},i=1,\dots,m.$

To incorporate the monotonic constraint $\mathcal{Q}$, we have

$\mathbf{\hat{y}}=\operatornamewithlimits{argmin}_{\mathbf{y}\in\mathcal{Q}}\left\\{-L(\mathbf{y},\tilde{\lambda})\right\\}=\operatornamewithlimits{argmin}_{\mathbf{y}\in\mathcal{Q}}\sum_{i=1}^{m}t_{1,(i)}\left[-\log
y_{i}-\frac{-\\{\sum_{i=1}^{m}t_{1,(i)}\\}\\{p_{1,(i)}-p_{1,(i-1)}\\}}{t_{1,(i)}}y_{i}\right].$

Let $\mathbf{u}=(u_{1},\dots,u_{m})$ and

$\hat{\mathbf{u}}=\operatornamewithlimits{argmin}_{\mathbf{u}\in\mathcal{Q}}\sum_{i=1}^{m}t_{1,(i)}\left[u_{i}-\frac{-\\{\sum_{i=1}^{m}t_{1,(i)}\\}\\{p_{1,(i)}-p_{1,(i-1)}\\}}{t_{1,(i)}}\right]^{2}.$

The solution takes the max-min form

$\hat{u}_{i}=\max_{b\geq i}\min_{a\leq
i}\frac{-\\{\sum_{i=1}^{m}t_{1,(i)}\\}\sum_{k=a}^{b}\\{p_{1,(k)}-p_{1,(k-1)}\\}}{\sum_{k=a}^{b}t_{1,(k)}}.$

According to Theorem 3.1 of Barlow and Brunk (1972), the final estimates of
(10) are given by $\hat{f}_{1}(p_{1,(i)})=\hat{y}_{i}=-\frac{1}{\hat{u}_{i}}$.
This can be implemented by PAVA (de Leeuw et al., 2009).

The estimation of $f_{2}$ can be obtained by solving the optimization problem
(11) in the same way, and we omit the details.

##  3 Theory

###  3.1 Notations

Let $f_{1}$ be the probability density function of $p_{1i}$ when
$\theta_{1i}=1$ and $f_{2}$ be the probability density function of $p_{2i}$
when $\theta_{2i}=1,i=1,\ldots,n.$ Define a two-dimensional functional space

$\mathcal{F}=\left\\{p_{w}(x,y)=\xi_{00}+\xi_{10}f_{1}(x)+\xi_{01}f_{2}(y)+\xi_{11}f_{1}(x)f_{2}(y):w\in\Omega,x,y\in(0,1)\right\\}.$
(12)

Let $w=(\xi_{00},\xi_{10},\xi_{01},\xi_{11},f_{1},f_{2})$ and
$w^{*}=(\xi^{*}_{00},\xi^{*}_{10},\xi^{*}_{01},\xi^{*}_{11},f_{1}^{*},f_{2}^{*})$.
Define a metric $\tau$ on $\Omega$ by

$\tau(w,w^{*})=|\xi_{00}-\xi_{00}^{*}|+|\xi_{10}-\xi_{10}^{*}|+|\xi_{01}-\xi_{01}^{*}|+|\xi_{11}-\xi_{11}^{*}|+d(f_{1},f_{1}^{*})+d(f_{2},f_{2}^{*}),$

where $d(f_{1},f_{1}^{*})=\int_{0}^{1}|f_{1}(x)-f_{1}^{*}(x)|dx$ and
$d(f_{2},f_{2}^{*})=\int_{0}^{1}|f_{2}(x)-f_{2}^{*}(x)|dx.$

We define the two-dimensional Hellinger distance between $p_{w}$ and
$p_{w^{*}}$ as follows.

$h(p_{w},p_{w^{*}})=\left[\frac{1}{2}\int_{0}^{1}\int_{0}^{1}\left\\{p_{w}(x,y)^{1/2}-p_{w^{*}}(x,y)^{1/2}\right\\}^{2}dxdy\right]^{1/2}.$
(13)

###  3.2 Identifiability

We first want to make sure that

$p_{w}(x,y)=\xi_{00}+\xi_{01}f_{2}\left(y\right)+\xi_{10}f_{1}\left(x\right)+\xi_{11}f_{1}\left(x\right)f_{2}\left(y\right)$

is identifiable, which means that when $h(p_{w},p_{w^{*}})=0$,
$\tau(w,w^{*})=0$.

Suppose $h(p_{w},p_{w^{*}})=0$, we have $p_{w}(x,y)=p_{w^{*}}(x,y)$ almost
everywhere. Integrating out $y$ and $x$ from $p_{w}(x,y)$ respectively, we get

$g_{w,1}(x)=\int_{0}^{1}p_{w}(x,y)dy=(\xi_{00}+\xi_{01})+(\xi_{10}+\xi_{11})f_{1}(x),$
(14)

and

$g_{w,2}(y)=\int_{0}^{1}p_{w}(x,y)dx=(\xi_{00}+\xi_{10})+(\xi_{01}+\xi_{11})f_{2}(y).$
(15)

The marginal distributions of $p_{w}$ and $p_{w^{*}}$ should be the same
almost everywhere. Specifically,

$g_{w,1}(x)=g_{w^{*},1}(x),\quad\mbox{and}\quad
g_{w,2}(y)=g_{w^{*},2}(y)\quad\text{a.e.}$

In addition, we have

$\displaystyle p_{w}(x,y)$
$\displaystyle=\xi_{00}+\frac{\xi_{10}}{\xi_{10}+\xi_{11}}\\{g_{w,1}(x)-(\xi_{00}+\xi_{01})\\}+\frac{\xi_{01}}{\xi_{01}+\xi_{11}}\\{g_{w,2}(y)-(\xi_{00}+\xi_{10})\\}$
$\displaystyle\quad+\frac{\xi_{11}}{(\xi_{10}+\xi_{11})(\xi_{01}+\xi_{11})}\\{g_{w,1}(x)-(\xi_{00}+\xi_{01})\\}\\{g_{w,2}(y)-(\xi_{00}+\xi_{10})\\}$
$\displaystyle=\frac{\xi_{00}\xi_{11}-\xi_{10}\xi_{01}}{(\xi_{10}+\xi_{11})(\xi_{01}+\xi_{11})}\\{1-g_{w,1}(x)\\}\\{1-g_{w,2}(y)\\}+g_{w,1}(x)g_{w,2}(y).$

As $p_{w}(x,y)=p_{w^{*}}(x,y)$ almost everywhere, we have

$\displaystyle\frac{\xi_{00}\xi_{11}-\xi_{10}\xi_{01}}{(\xi_{10}+\xi_{11})(\xi_{01}+\xi_{11})}\\{1-g_{w,1}(x)\\}\\{1-g_{w,2}(y)\\}+g_{w,1}(x)g_{w,2}(y)$
$\displaystyle=$
$\displaystyle\frac{\xi^{*}_{00}\xi^{*}_{11}-\xi^{*}_{10}\xi^{*}_{01}}{(\xi^{*}_{10}+\xi^{*}_{11})(\xi^{*}_{01}+\xi^{*}_{11})}\\{1-g_{w^{*},1}(x)\\}\\{1-g_{w^{*},2}(y)\\}+g_{w^{*},1}(x)g_{w^{*},2}(y).$

Consequently $(\xi_{00},\xi_{10},\xi_{01},\xi_{11})$ and
$(\xi_{00}^{*},\xi_{10}^{*},\xi_{01}^{*},\xi_{11}^{*})$ have to satisfy the
following conditions.

  1. (i)

$\xi_{00}+\xi_{10}+\xi_{01}+\xi_{11}=\xi_{00}^{*}+\xi_{10}^{*}+\xi_{01}^{*}+\xi_{11}^{*}=1$

  2. (ii)

$\frac{\xi_{00}\xi_{11}-\xi_{10}\xi_{01}}{(\xi_{10}+\xi_{11})(\xi_{01}+\xi_{11})}=\frac{\xi^{*}_{00}\xi^{*}_{11}-\xi^{*}_{10}\xi^{*}_{01}}{(\xi^{*}_{10}+\xi^{*}_{11})(\xi^{*}_{01}+\xi^{*}_{11})}$
is a constant.

We summarize the identifiability in the following theorem.

######  Theorem 1.

Let $w=(\xi_{00},\xi_{01},\xi_{10},\xi_{11},f_{1},f_{2})$, and
$w^{*}=(\xi_{00}^{*},\xi_{01}^{*},\xi_{10}^{*},\xi_{11}^{*},f_{1}^{*},f_{2}^{*})$.
Assume

$\lim\limits_{x\to 1}f_{1}(x)=\lim\limits_{y\to
1}f_{2}(y)=0,\quad\mbox{and}\quad\lim\limits_{x\to
1}f^{*}_{1}(x)=\lim\limits_{y\to 1}f^{*}_{2}(y)=0.$

If $p_{w}(x,y)=p_{w^{*}}(x,y)$ almost everywhere, we have

$\xi_{00}=\xi_{00}^{*},\xi_{10}=\xi_{10}^{*},\xi_{01}=\xi_{01}^{*},\xi_{11}=\xi_{11}^{*},$

and

$f_{1}(x)=f_{1}^{*}(x)\quad a.e.,\quad f_{2}(x)=f_{2}^{*}(x)\quad a.e.$

Therefore, in order to identify individual parameters and functions, we need
the following condition.

(C1) Identifiability condition: for all
$w=(\xi_{00},\xi_{10},\xi_{01},\xi_{11},f_{1},f_{2})\in\Omega$,

$\lim\limits_{x\to 1}f_{1}(x)=\lim\limits_{y\to 1}f_{2}(y)=0.$

###  3.3 Consistency

This section uses tools from the empirical processes(Van de Geer, 2000; van
der Vaart and Wellner, 1996) to get asymptotic results of the estimated
parameters and functions. Recall $\mathcal{F}$ defined in (12). Let
$\alpha\geq 1.$ Define a metric $\|\cdot\|_{\alpha}$ on $\mathcal{F}$ as

$\|f\|_{\alpha}=\left(\int_{0}^{1}|f(x)|^{\alpha}dx\right)^{1/\alpha}.$

If $\|f_{i}^{U}-f_{i}^{L}\|_{\alpha}\leq\epsilon,\forall i=1,\cdots,n,$ and
for any $f\in\mathcal{F},\exists i$ such that $f^{L}_{i}\leq f\leq f^{U}_{i}$,
we call $[f^{L}_{1},f^{U}_{1}],\cdots,[f^{L}_{n},f^{U}_{n}]$ an $\epsilon$
bracket. The minimal number of $\epsilon$ brackets that cover $\mathcal{F}$ is
called the bracketing number, denoted as
$N_{[]}(\epsilon,\mathcal{F},\|\cdot\|_{\alpha}).$ The bracketing entropy is
$\log N_{[]}(\epsilon,\mathcal{F},\|\cdot\|_{\alpha}).$

Assume that the true parameter is $w_{0}\in\Omega.$ Denote

$\mathcal{\bar{F}}^{1/2}=\\{\sqrt{\frac{p_{w}+p_{w_{0}}}{2}},w\in\Omega\\}.$

To get the rate of convergence of the MLE based on the Hellinger distance
(13), we use the bracketing entropy of $\mathcal{\bar{F}}^{1/2}$ employed with
the metric $\|\cdot\|_{2}.$

To bound $\log N_{[]}(\epsilon,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2}),$ we
need the following assumption.

(C2) There exists an envelope function $F(x)>0,x\in(0,1)$ such that
$\int_{0}^{1}F(x)dx<\infty$ and $\int_{0}^{1}F(x)^{2}dx<C,$ where $C>0$ is a
constant. Moreover, for all

$w=(\xi_{00},\xi_{10},\xi_{01},\xi_{11},f_{1},f_{2})\in\Omega,$

we have $f_{1}(x)<F(x),f_{2}(y)<F(y)$ and $\xi_{00}\geq c,$ where $c>0$ is a
constant. The following proposition gives the bracketing entropy of
$\bar{\mathcal{F}}^{1/2}$ equipped with the $L_{2}$ norm.

######  Proposition 1.

Under (C2), there exist $C_{1}>0$ and $C_{2}>0$ such that

$C_{2}\epsilon^{-1}\leq\log
N_{[]}(\epsilon,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2})\leq
C_{1}\epsilon^{-1}.$

######  Remark 1.

For the space $\mathcal{G}_{1}$ defined in (8), we allow the envelop function
to diverge at 0. Without additional assumptions, the bracketing number of
$\mathcal{G}_{1}$ may be infinite (Gao and Wellner, 2009). We assume that the
envelop function of $\mathcal{G}_{1}$ has a finite $L_{2}$ norm specified in
(C2). Under (C2), we are able to get
$N_{[]}(\epsilon,\mathcal{F},\|\cdot\|_{2})$ without requiring that the
envelop function $F(x)$ is bounded.

Proposition 1 calculates the bracketing entropy, which allows us to get the
convergence rate of the Hellinger distance between $p_{\hat{w}_{n}}$ and
$p_{w_{0}}.$

Our main result is the following theorem.

######  Theorem 2.

Under (C2), suppose $\mathcal{F}$ is the space specified in (12). Denote
$p_{w_{0}}\in\mathcal{F}$ as the true density function, and $p_{\hat{w}_{n}}$
as the MLE. Then $\exists C_{3},C_{4},C_{5}>0$ such that

$P\left(h\left(p_{\hat{w}_{n}},p_{w_{0}}\right)>C_{3}n^{-1/3}\right)\leq
C_{4}\exp\left[-C_{5}n^{1/3}\right].$

In particular,

$h(p_{\hat{w}_{n}},p_{w_{0}})=O_{p}(n^{-1/3}).$

In Gao and Wellner (2007), they calculated the bracketing entropy with respect
to $\|\cdot\|_{2}$ for uniformly bounded monotone functions supported on
$[0,1]^{2}.$ The order of bracketing entropy is $\epsilon^{-2}(\log
1/\epsilon)^{2},$ resulting $n^{-\frac{1}{4}}\log n$ rate of convergence for
MLE in Hellinger distance. The functional space $\Omega$ (7) we are working
with is much smaller than that in Gao and Wellner (2007). Consequently, we get
a faster rate of convergence. Furthermore, for log-concave density functions,
the order of bracketing entropy is $\epsilon^{-1}\log^{3/2}(1/\epsilon)$
resulting in $n^{-1/3}\log^{1/2}n$ rate of convergence for MLE in Hellinger
distance (Kim and Samworth, 2016). Next, we show the consistency of
$\hat{w}_{n}.$

######  Theorem 3.

Under (C1) and (C2), we have $\tau(\hat{w}_{n},w_{0})\overset{p}{\to}0,$ where
$\hat{w}_{n}$ is the MLE, $w_{0}$ is the true parameter and $\overset{p}{\to}$
means convergence in probability.

###  3.4 False discovery rate control

With MLE
$\hat{w}_{n}=(\hat{\xi}_{00},\hat{\xi}_{01},\hat{\xi}_{10},\hat{\xi}_{11},\hat{f}_{1},\hat{f}_{2}),$
we have

$\widehat{\mbox{Lfdr}}(p_{1i},p_{2i})=\frac{\hat{\xi}_{00}+\hat{\xi}_{01}\hat{f}_{2}\left(p_{2i}\right)+\hat{\xi}_{10}\hat{f}_{1}\left(p_{1i}\right)}{\hat{\xi}_{00}+\hat{\xi}_{01}\hat{f}_{2}\left(p_{2i}\right)+\hat{\xi}_{10}\hat{f}_{1}\left(p_{1i}\right)+\hat{\xi}_{11}\hat{f}_{1}\left(p_{1i}\right)\hat{f}_{2}\left(p_{2i}\right)}.$
(16)

With (5) and (16), at threshold $\lambda,$ we write the estimated total number
of rejections

$\hat{R}_{m}(\lambda)=\sum_{i=1}^{m}\mathbf{I}\left\\{\widehat{\operatorname{Lfdr}}(p_{1i},p_{2i})\leq\lambda\right\\}.$

Note that

$E\left[\mbox{Lfdr}(p_{1i},p_{2i})I\\{\mbox{Lfdr}(p_{1i},p_{2i})\leq\lambda\\}\right]=E[I\\{\mbox{Lfdr}(p_{1i},p_{2i})\leq\lambda\\}(1-\theta_{1i}\theta_{2i})].$
(17)

we have the estimated number of false rejections

$\hat{V}_{m}(\lambda)=\sum_{i=1}^{m}\widehat{\operatorname{Lfdr}}(p_{1i},p_{2i})\mathbf{I}\left\\{\widehat{\operatorname{Lfdr}}(p_{1i},p_{2i})\leq\lambda\right\\},$

the estimated false discovery proportion

$\widehat{\mbox{FDP}}_{m}(\lambda)=\frac{\hat{V}_{m}(\lambda)}{\hat{R}_{m}(\lambda)},$

and the critical value at false discovery rate level $\alpha$

$\hat{\lambda}_{m}=\sup\\{\lambda\in[0,1]:\widehat{\mbox{FDP}}_{m}(\lambda)\leq\alpha\\}.$
(18)

We reject $H_{0i}$ when
$\widehat{\mbox{Lfdr}}(p_{1i},p_{2i})\leq\hat{\lambda}_{m}.$ In practice, we
use a step-up procedure (Sun and Cai, 2007). Specifically, we order
$\widehat{\mbox{Lfdr}}(p_{1i},p_{2i})$ from the smallest to the largest and
denote the order statistics
$\widehat{\mbox{Lfdr}}_{(1)}\leq\cdots\leq\widehat{\mbox{Lfdr}}_{(m)}$ and
corresponding ordered hypothesis $H^{(1)},\cdots,H^{(m)}.$ Reject
$H^{(1)},\cdots H^{(j)},$ where

$j=\mbox{max}\left\\{1\leq k\leq
m:\frac{1}{k}\sum_{i=1}^{k}\widehat{\mbox{Lfdr}}_{(i)}\leq\alpha\right\\}.$

We need additional conditions.

(C3) $f_{1}$ and $f_{2}$ are continuous non-increasing density functions on
$(0,1)$.

(C4) $\forall\lambda\in(0,1)$,

$\begin{gathered}\frac{1}{m}\sum_{i=1}^{m}\mathbf{I}\left\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\right\\}\overset{p}{\to}B_{1}(\lambda),\\\
\frac{1}{m}\sum_{i=1}^{m}\operatorname{Lfdr}(p_{1i},p_{2i})\mathbf{I}\left\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\right\\}\overset{p}{\to}B_{2}(\lambda),\end{gathered}$
$\frac{1}{m}\sum\limits_{i=1}^{m}I\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\\}(1-\theta_{1i}\theta_{2i})\overset{p}{\to}B_{2}(\lambda),$
(19)

where $B_{1}(\lambda)$ and $B_{2}(\lambda)$ are continuous functions on
$(0,1)$.

(C5) $\exists\lambda_{0}\in(0,1),$ such that $B_{1}(\lambda_{0})>0$, and
$\frac{B_{2}(\lambda_{0})}{B_{1}(\lambda_{0})}\leq\alpha.$

In contrast to the predominant parametric assumptions on the $p$-value density
functions under non-null, (C3) is a mild assumption. (C4) is similar to the
assumptions in Storey et al. (2004). In view of (17), (19) follows from the
law of large numbers. (C5) ensures the existence of critical value
$\hat{\lambda}_{m}$ in the population level.

######  Theorem 4.

Under (1), (C1)-(C5), we have

$\limsup\limits_{m\to\infty}\operatorname{FDR}(\hat{\lambda}_{m})\leq\alpha,$

where $\hat{\lambda}_{m}$ is chosen based on (18) and FDR is specified in (6).

######  Remark 2.

We estimate unknown parameters and unknown functions based on the MLE, where
we assume independence. In our simulation studies, we show that our method is
robust to dependent data. We conjecture that similar results would hold when
the data are weakly dependent and leave it for future research.

##  4 Simulation studies

###  4.1 Independent data

We perform simulations to assess the FDR control and power of the proposed
method. The data generation strategy is as follows. Given pre-specified prior
probabilities $(\xi_{00},\xi_{01},\xi_{10},\xi_{11})$, we generate the joint
hidden states from a multinomial distribution with probabilities
$P(\theta_{1i}=k,\theta_{2i}=l)=\xi_{kl}$ for $k,l\in\\{0,1\\}$. The test
statistics for study $j\ (j=1,2)$ are generated from
$X_{ji}\mid\theta_{ji}\sim\theta_{ji}N(0,\sigma_{j}^{2})+(1-\theta_{ji})N(\mu_{j},\sigma_{j}^{2}),\
i=1,\dots,m$, where $N(\mu,\sigma^{2})$ denotes the normal distribution with
mean $\mu$ and variance $\sigma^{2}$. The one-sided $p$-values are generated
from $p_{ji}=1-\Phi\left(\frac{X_{ji}}{\sigma_{j}}\right)$ for $i=1,\dots,m$
and $j=1,2$. We compare the proposed method with competing replicability
analysis methods, including ad hoc BH, MaxP (Benjamini et al., 2009), IDR (Li
et al., 2011), MaRR (Philtron et al., 2018), radjust (Bogomolov and Heller,
2018), and JUMP (Lyu et al., 2023). Details of the competing methods can be
found in the Supplementary Material.

Figure 1: FDR control of different methods for independent data. Figure 2:
Power comparison of different methods for independent data.

We set $m=10,000,\xi_{11}=0.01$ and $\xi_{01}=\xi_{10}$, such that with a
given value of $\xi_{00}$, we have
$\xi_{01}=\xi_{10}=(1-\xi_{00}-\xi_{11})/2$. We vary
$\mu_{1},\mu_{2},\sigma_{1},\sigma_{2}$ and $\xi_{00}$ across simulations to
evaluate performance. For each setting, the FDR and power of different methods
are calculated based on $100$ simulated datasets. We present the FDR control
and power comparison at the FDR level $0.05$ with different methods in Figure
1 and Figure 2. The horizontal dashed lines in Figure 1 indicate the target
FDR level. We observe that ad hoc BH and IDR do not control the FDR. The other
methods control the FDR across all settings, where MaxP is overly
conservative. MaRR works well when $\xi_{00}=0.85$, but loses power when
$\xi_{01}$ and $\xi_{10}$ are relatively large ($\xi_{00}=0.55,0.65$ or
$0.8$). JUMP and radjust have decent power, and our proposed method shows the
highest power across all settings.

###  4.2 Dependent data

We further simulate data with local dependence to examine the robustness of
the proposed method. Specifically, in each study, we divide the $10,000$
hypotheses into $100$ blocks and each block is further divided into two sub-
blocks of equal size. Within each sub-block, there is a constant positive
correlation (e.g., $0.2$) and a constant negative correlation (e.g., $-0.2$),
respectively. Data in different blocks are independent. In each block of study
$j$ ($j=1,2$), the corresponding hidden states are
$\bm{\theta}_{j}^{*}=(\theta_{jk}^{*})_{k=1}^{100}$, the test statistics are
generated from a multivariate normal distribution:
$\bm{X}_{j}\sim\mbox{MVN}(\bm{\mu},\Sigma)$, where the mean vector
$\bm{\mu}=\bm{\theta}^{*}_{j}\cdot\mu_{j},$ where $\mu_{j}$ is the same as
that in Section 4.1. We use a $4$-dimensional covariance matrix with
$\rho=0.2$ and $-0.2$ to illustrate the pattern of $\Sigma,$

$\Sigma=\begin{pmatrix}1&0.2&-0.2&-0.2\\\ 0.2&1&-0.2&-0.2\\\
-0.2&-0.2&1&0.2\\\ -0.2&-0.2&0.2&1\end{pmatrix}.$

Figure 3: FDR control of different methods for dependent data. Figure 4: Power
comparison of different methods for dependent data.

We present the FDR control and power comparison at the level $0.05$ with
different methods for dependent data in Figure 3 and Figure 4. Except ad hoc
BH and IDR, all methods control the FDR at the target level across all
settings. As to the power, we observe a similar pattern to the independent
data. Our proposed method has higher power than competing methods.

##  5 Data analysis

We illustrate our method by analyzing data from genome-wide association
studies (GWAS). Mahajan et al. (2022) assembled an ancestrally diverse
collection of 122 GWAS in $180,834$ type 2 diabetes (T2D) cases and
$1,159,055$ controls across five ancestry groups. We conduct replicability
analysis of the East Asian ancestry and European ancestry specific GWAS to
identify replicable T2D-associated SNPs across two populations. The datasets
can be downloaded from the DIAbetes Genetics Replication And Meta-analysis
Consortium at https://www.diagram-consortium.org/downloads.html. The East
Asian dataset includes summary statistics of $760,565$ SNPs, and the European
dataset includes summary statistics of $879,506$ SNPs in Chromosome 2. We
analyze the paired $p$-values of $760,565$ SNPs that are common in both
studies, where $p_{1i},i=1,\ldots,m$ denote $p$-values of European ancestry
and $p_{2i},i=1,\ldots,m$ denote $p$-values of East Asian ancestry.

Figure 5: Replicability analysis results of the T2D GWAS data.

We compare our method to competing replicability analysis methods for
comparison. The IDR method is not included, as it is not computationally
scalable and cannot complete the analysis in $10$ hours. The results of
different methods at the FDR level $0.01$ are summarized in Figure 5. MaxP is
the most conservative that identified $278$ replicable SNPs, $274$ of which
were also detected by other methods. Our method identified 622 replicable
SNPs, including 280 SNPs that were not detected by any other methods. This is
consistent with the simulation studies, where improved power of our method has
been observed.

Figure 6: Estimates of the non-null probability density functions in two T2D
studies

More significant SNPs are expected to be detected in the European population
than in the East Asian population, as the European ancestry dataset has a much
larger sample size. Our method demonstrated this by estimating the proportions
of four states as
$\hat{\xi}_{00}=0.8923,\hat{\xi}_{01}=0.0984,\hat{\xi}_{10}=0.0003$ and
$\hat{\xi}_{11}=0.0089$. The estimated non-null probability density functions
in the two studies, namely $\hat{f}_{1}$ and $\hat{f}_{2}$, are plotted in
Figure 6. Furthermore, we examine the $280$ replicable SNPs only identified by
our method in SNPnexus, which is a Web server-based annotation tool for
functional annotation of human genome sequence variation (Oscanoa et al.,
2020). To interpret our results, we map the identified SNPs to the genes
containing these SNPs or the genes closest to these SNPs. The $280$ SNPs are
mapped to $21$ genes, many of which have been reported as related to T2D in
previous literature. For instance, the GCKR gene is significantly related in
multiple ways to T2D and kidney-related phenotypes (Ho et al., 2022). BCL11A
deficiency is associated with decreased fetal hemoglobin (Xu et al., 2013),
which is significantly associated with a decreased risk of T2D. BCL11A gene
variants affect the insulin response to glucose and glucagon secretion, and
thus have been shown to increase the risk of T2D (Cauchi et al., 2012). The
Grb14 gene has been identified as an inhibitor of insulin signaling and action
and is involved in insulin resistance associated with T2D and obesity (Gondoin
et al., 2017).

##  6 Concluding remarks

In this paper, we propose an empirical Bayes approach for robust and powerful
inference of high dimensional replicability analysis. We account for the
heterogeneity of different studies by jointly modeling different states and
allowing non-null density functions to have different distributions.
Furthermore, we use a non-parametric approach for density function estimation
under monotonic constraint to achieve robust inference. By enumerating
different components of the composite null hypothesis of replicability
analysis, we are able to borrow information from different studies and
different features to improve power. Extensive simulation studies demonstrate
favorable performance over competing methods. An analysis of GWAS data sets
reveals novel biological findings that competing methods may not obtain.

We assume that conditional on the joint hidden states, $p$-values are
independent. This can be relaxed to allow weak dependence, which was
demonstrated in our simulation studies. The theoretical development is more
involved, and we leave it to future research. For strong dependences, such as
those induced by the hidden Markov model, the proposed method is not
applicable, and new tools are needed.

We use significance to assess replicability. We are aware of the controversy
over the use of $p$-values. Although acknowledging the drawbacks of
$p$-values, we use them due to their versatility and easy availability. Due to
privacy concerns and/or logistics, it is difficult to obtain raw data, yet
summary-level statistics are readily available. In addition, $p$-values can be
obtained in various ways, including asymptotic calibration, permutations, or
finite sample calculation. Our method is platform-independent and is
applicable to various study designs, and only needs $p$-values as input.

The empirical Bayes approach enumerates composite nulls, which makes the
method not scalable. In principle, our method can handle more than two
studies. In practice, the replicability analysis of $n$ studies requires the
estimation of $2^{n}$ mixtures, which can be computationally prohibitive.

## References

  * Barlow and Brunk (1972) Barlow, R. E. and H. D. Brunk (1972).  The isotonic regression problem and its dual.  Journal of the American Statistical Association 67(337), 140–147.
  * Benjamini et al. (2009) Benjamini, Y., R. Heller, and D. Yekutieli (2009).  Selective inference in complex research.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 367(1906), 4255–4271.
  * Benjamini and Hochberg (1995) Benjamini, Y. and Y. Hochberg (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.
  * Blei et al. (2007) Blei, R., F. Gao, and W. Li (2007).  Metric entropy of high dimensional distributions.  Proceedings of the American Mathematical Society 135(12), 4009–4018.
  * Bogomolov and Heller (2013) Bogomolov, M. and R. Heller (2013).  Discovering findings that replicate from a primary study of high dimension to a follow-up study.  Journal of the American Statistical Association 108(504), 1480–1492.
  * Bogomolov and Heller (2018) Bogomolov, M. and R. Heller (2018).  Assessing replicability of findings across two studies of multiple features.  Biometrika 105(3), 505–516.
  * Bogomolov and Heller (2022) Bogomolov, M. and R. Heller (2022).  Replicability across multiple studies.  arXiv preprint arXiv:2210.00522.
  * Cao et al. (2013) Cao, H., W. Sun, and M. R. Kosorok (2013).  The optimal power puzzle: Scrutiny of the monotone likelihood ratio assumption in multiple testing.  Biometrika 100(2), 495–502.
  * Cauchi et al. (2012) Cauchi, S., I. Ezzidi, Y. El Achhab, N. Mtiraoui, L. Chaieb, D. Salah, C. Nejjari, Y. Labrune, L. Yengo, D. Beury, et al. (2012).  European genetic variants associated with type 2 diabetes in north african arabs.  Diabetes & Metabolism 38(4), 316–323.
  * de Leeuw et al. (2009) de Leeuw, J., K. Hornik, and P. Mair (2009).  Isotone optimization in r: Pool-adjacent-violators algorithm (pava) and active set methods.  Journal of Statistical Software 32(5), 1–24.
  * Dempster et al. (1977) Dempster, A. P., N. M. Laird, and D. B. Rubin (1977).  Maximum likelihood from incomplete data via the em algorithm.  Journal of the Royal Statistical Society: Series B (Methodological) 39(1), 1–22.
  * Efron (2012) Efron, B. (2012).  Large-scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction.  Cambridge University Press.
  * Freedman et al. (2015) Freedman, L. P., I. M. Cockburn, and T. S. Simcoe (2015).  The economics of reproducibility in preclinical research.  PLoS Biology 13(6), e1002165.
  * Gao and Wellner (2007) Gao, F. and J. A. Wellner (2007).  Entropy estimate for high-dimensional monotonic functions.  Journal of Multivariate Analysis 98(9), 1751–1764.
  * Gao and Wellner (2009) Gao, F. and J. A. Wellner (2009).  On the rate of convergence of the maximum likelihood estimator of a k-monotone density.  Science in China Series A: Mathematics 52(7), 1525–1538.
  * Gondoin et al. (2017) Gondoin, A., C. Hampe, R. Eudes, C. Fayolle, C. Pierre-Eugène, M. Miteva, B. O. Villoutreix, F. Charnay-Pouget, D. J. Aitken, T. Issad, et al. (2017).  Identification of insulin-sensitizing molecules acting by disrupting the interaction between the insulin receptor and grb14.  Scientific Reports 7(1), 16901.
  * Grenander (1956) Grenander, U. (1956).  On the theory of mortality measurement: part ii.  Scandinavian Actuarial Journal 1956(2), 125–153.
  * Hanche-Olsen and Holden (2010) Hanche-Olsen, H. and H. Holden (2010).  The kolmogorov–riesz compactness theorem.  Expositiones Mathematicae 28(4), 385–394.
  * Heller and Yekutieli (2014) Heller, R. and D. Yekutieli (2014).  Replicability analysis for genome-wide association studies.  The Annals of Applied Statistics 8(1), 481–498.
  * Ho et al. (2022) Ho, L.-J., C.-H. Lu, R.-Y. Su, F.-H. Lin, S.-C. Su, F.-C. Kuo, N.-F. Chu, Y.-J. Hung, J.-S. Liu, and C.-H. Hsieh (2022).  Association between glucokinase regulator gene polymorphisms and serum uric acid levels in taiwanese adolescents.  Scientific Reports 12(1), 5519.
  * Hung and Fithian (2020) Hung, K. and W. Fithian (2020).  Statistical methods for replicability assessment.  The Annals of Applied Statistics 14(3), 1063–1087.
  * Ioannidis (2005) Ioannidis, J. P. (2005).  Why most published research findings are false.  PLoS Medicine 2(8), e124.
  * James Y. Dai and LeBlanc (2022) James Y. Dai, J. L. S. and M. LeBlanc (2022).  A multiple-testing procedure for high-dimensional mediation hypotheses.  Journal of the American Statistical Association 117(537), 198–213.
  * Kim and Samworth (2016) Kim, A. K. and R. J. Samworth (2016).  Global rates of convergence in log-concave density estimation.  The Annals of Statistics 44(6), 2756–2779.
  * Li et al. (2011) Li, Q., J. B. Brown, H. Huang, and P. J. Bickel (2011).  Measuring reproducibility of high-throughput experiments.  Annals of Applied Statistics 5(3), 1752–1779.
  * Liu et al. (2022) Liu, Z., J. Shen, R. Barfield, J. Schwartz, A. A. Baccarelli, and X. Lin (2022).  Large-scale hypothesis testing for causal mediation effects with applications in genome-wide epigenetic studies.  Journal of the American Statistical Association 117(537), 67–81.
  * Lyu et al. (2023) Lyu, P., Y. Li, X. Wen, and H. Cao (2023).  Jump: Replicability analysis of high-throughput experiments with applications to spatial transcriptomic studies.  Bioinformatics 39(6), 366.
  * Mahajan et al. (2022) Mahajan, A., C. N. Spracklen, W. Zhang, M. C. Ng, L. E. Petty, H. Kitajima, G. Z. Yu, S. Rüeger, L. Speidel, Y. J. Kim, et al. (2022).  Multi-ancestry genetic study of type 2 diabetes highlights the power of diverse populations for discovery and translation.  Nature Genetics 54(5), 560–572.
  * Oscanoa et al. (2020) Oscanoa, J., L. Sivapalan, E. Gadaleta, A. Z. Dayem Ullah, N. R. Lemoine, and C. Chelala (2020).  Snpnexus: a web server for functional annotation of human genome sequence variation (2020 update).  Nucleic Acids Research 48(W1), W185–W192.
  * Philtron et al. (2018) Philtron, D., Y. Lyu, Q. Li, and D. Ghosh (2018).  Maximum rank reproducibility: a nonparametric approach to assessing reproducibility in replicate experiments.  Journal of the American Statistical Association 113(523), 1028–1039.
  * Resnick (2005) Resnick, S. (2005).  A Probability Path.  Springer.
  * Robertson et al. (1988) Robertson, T., R. L. Dykstra, and F. T. Wright (1988).  Order restricted statistical inference.  In Wiley Series in Probability and Mathematical Statistics. John Wiley and Sons.
  * Stein and Shakarchi (2009) Stein, E. M. and R. Shakarchi (2009).  Real Analysis: Measure Theory, Integration, and Hilbert spaces.  Princeton University Press.
  * 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. (2004) Storey, J. D., J. E. Taylor, and D. Siegmund (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.
  * Storey and Tibshirani (2003) Storey, J. D. and R. Tibshirani (2003).  Statistical significance for genomewide studies.  Proceedings of the National Academy of Sciences 100(16), 9440–9445.
  * Sun and Cai (2007) Sun, W. and T. T. Cai (2007).  Oracle and adaptive compound decision rules for false discovery rate control.  Journal of the American Statistical Association 102(479), 901–912.
  * Van de Geer (1993) Van de Geer, S. (1993).  Hellinger-consistency of certain nonparametric maximum likelihood estimator.  The Annals of Statistics 21(1), 14–44.
  * Van de Geer (2000) Van de Geer, S. (2000).  Empirical Processes in M-estimation.  Cambridge University Press.
  * van der Vaart and Wellner (1996) van der Vaart, A. W. and J. A. Wellner (1996).  Weak Convergence and Empirical Processes: with Applications to Statistics.  Springer New York, NY.
  * Xu et al. (2013) Xu, J., D. E. Bauer, M. A. Kerenyi, T. D. Vo, S. Hou, Y.-J. Hsu, H. Yao, J. J. Trowbridge, G. Mandel, and S. H. Orkin (2013).  Corepressor-dependent silencing of fetal hemoglobin expression by bcl11a.  Proceedings of the National Academy of Sciences 110(16), 6518–6523.
  * Zhao et al. (2020) Zhao, Y., M. G. Sampson, and X. Wen (2020).  Quantify and control reproducibility in high-throughput experiments.  Nature Methods 17(12), 1207–1213.

## Supplementary material

##  Appendix A Proofs

###  A.1 Proof of Theorem 1

###### Proof.

Assume that $p_{w}(x,y)=p_{w^{*}}(x,y)$ almost everywhere, where

$w=(\xi_{00},\xi_{01},\xi_{10},\xi_{11},f_{1},f_{2}),\quad\text{and}\quad
w^{*}=(\xi_{00}^{*},\xi_{01}^{*},\xi_{10}^{*},\xi_{11}^{*},f_{1}^{*},f_{2}^{*}).$

Notice that when $x\to 1$,

$\lim\limits_{x\to 1}f_{1}(x)=\lim\limits_{x\to
1}f_{1}^{*}(x)=\lim\limits_{x\to 1}f_{2}(x)=\lim\limits_{x\to
1}f_{2}^{*}(x)=0.$

By (14) and (15), when $x\to 1$ we have

$\displaystyle\xi_{00}+\xi_{10}$
$\displaystyle=\xi^{*}_{00}+\xi^{*}_{10}=\lim\limits_{x\to 1}g_{w,1}(x)$ (A.1)
$\displaystyle\xi_{00}+\xi_{01}$
$\displaystyle=\xi^{*}_{00}+\xi^{*}_{01}=\lim\limits_{x\to 1}g_{w,2}(x)$ (A.2)

By the analysis in section 3.2, the relationship between $w$ and $w^{*}$ is

$\displaystyle\frac{\xi_{00}\xi_{11}-\xi_{10}\xi_{01}}{(\xi_{10}+\xi_{11})(\xi_{01}+\xi_{11})}$
$\displaystyle=\frac{\xi^{*}_{00}\xi^{*}_{11}-\xi^{*}_{10}\xi_{11}^{*}}{(\xi_{10}^{*}+\xi_{11}^{*})(\xi_{01}^{*}+\xi_{11}^{*})}$
(A.3) $\displaystyle\xi_{00}+\xi_{10}+\xi_{01}+\xi_{11}$
$\displaystyle=\xi_{00}^{*}+\xi^{*}_{10}+\xi^{*}_{01}+\xi^{*}_{11}=1.$ (A.4)

By (A.1),(A.2),(A.4), we get

$\xi_{01}+\xi_{11}=\xi^{*}_{01}+\xi^{*}_{11},\quad\xi_{10}+\xi_{11}=\xi^{*}_{10}+\xi^{*}_{11}.$
(A.5)

Plugging into (A.3), we have

$\xi_{00}\xi_{11}-\xi_{10}\xi_{01}=\xi^{*}_{00}\xi^{*}_{11}-\xi^{*}_{10}\xi^{*}_{01}.$

Note that

$\xi_{00}\xi_{11}-\xi_{10}\xi_{01}=\xi_{00}\xi_{11}+\xi_{00}\xi_{01}-\xi_{00}\xi_{01}-\xi_{10}\xi_{01}=\xi_{00}(\xi_{01}+\xi_{11})-\xi_{01}(\xi_{00}+\xi_{10}).$

We have

$\xi_{00}(\xi_{01}+\xi_{11})-\xi_{01}(\xi_{00}+\xi_{10})=\xi^{*}_{00}(\xi^{*}_{01}+\xi^{*}_{11})-\xi^{*}_{01}(\xi^{*}_{00}+\xi^{*}_{10}).$
(A.6)

By(A.1), (A.5), and (A.6), we obtain

$(\xi_{00}-\xi^{*}_{00})(\xi_{01}+\xi_{11})-(\xi_{01}-\xi^{*}_{01})(\xi_{00}+\xi_{10})=0.$
(A.7)

By (A.2) and (A.7), we have $\xi_{00}-\xi^{*}_{00}=-(\xi_{01}-\xi^{*}_{01})$.
Therefore,

$(\xi_{01}-\xi_{01}^{*})(\xi_{00}+\xi_{10}+\xi_{01}+\xi_{11})=\xi_{01}-\xi^{*}_{01}=0.$
(A.8)

Plugging (A.8) into (A.1), (A.2), and (A.4), we have

$\xi_{00}=\xi^{*}_{00},\xi_{10}=\xi^{*}_{10},\xi_{01}=\xi^{*}_{01},\xi_{11}=\xi^{*}_{11}.$

By (14) and (15), we obtain

$f_{1}(x)=f_{1}^{*}(x),\quad f_{2}(y)=f_{2}^{*}(y)\quad a.e.$

∎

###  A.2 Proof of Proposition 1

Under conditions (C1) and (C2), we rewrite $\mathcal{F}$ as follows.

$\mathcal{F}=\\{p_{w}=\xi_{00}+\xi_{10}f_{1}(x)+\xi_{01}f_{2}(y)+\xi_{11}f_{1}(x)f_{2}(y):w=(\xi_{00},\xi_{01},\xi_{10},\xi_{11},f_{1}(x),f_{2}(y))\in\Omega\\},$

where

$\Omega=\Delta\times\mathcal{G}_{1}\times\mathcal{G}_{2},$

where

$\displaystyle\Delta$
$\displaystyle=\\{(\xi_{00},\xi_{10},\xi_{01},\xi_{11})\in
R^{4}_{+}:\xi_{00}+\xi_{10}+\xi_{01}+\xi_{11}=1,\xi_{00}\geq c\\},$
$\displaystyle\mathcal{G}_{1}$
$\displaystyle=\\{f_{1}:f_{1}\in\mathcal{G},f_{1}\leq F,\lim\limits_{x\to
1}f_{1}(x)=0\\},$ $\displaystyle\mathcal{G}_{2}$
$\displaystyle=\\{f_{2}:f_{2}\in\mathcal{G},f_{2}\leq F,\lim\limits_{y\to
1}f_{2}(y)=0\\}.$

Without loss of generality, we further assume that $F(x)$ is non-increasing
and large enough such that $F(x)>3$, and $c$ is small enough such that
$c<1/5$.

We first show that the bracketing number of $\mathcal{\bar{F}}^{1/2}$ can be
bounded by the bracketing number of $\mathcal{F}$, summarized in Proposition
A.1. Then we calculate the upper bound of
$N_{[]}(\epsilon,\mathcal{F},\|\cdot\|_{2})$, and the lower bound of
$N_{[]}(\epsilon,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2}),$ respectively.

######  Proposition A.1.

Under (C2), we have

$N_{[]}(\epsilon,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2})\leq
N_{[]}(2\sqrt{2c}\epsilon,\mathcal{F},\|\cdot\|_{2}).$

####  A.2.1 Proof of Proposition A.1

Let $[p^{L}_{1},p^{U}_{1}],\cdots,[p^{L}_{n},p^{U}_{n}]$ be the smallest
$\epsilon$-brackets of $\mathcal{F}.$ Without loss of generality, we further
assume that $p^{U}_{i},p^{L}_{i}\geq 0,i=1,\ldots,n,$ because every function
in the functional space $\mathcal{F}$ is bounded away from $c>0$. We have

$\left|\sqrt{\frac{p^{U}_{i}+p_{w_{0}}}{2}}-\sqrt{\frac{p^{L}_{i}+p_{w_{0}}}{2}}\right|=\frac{1}{\sqrt{2}}\frac{|p_{i}^{U}-p_{i}^{L}|}{\sqrt{p^{U}_{i}+p_{w_{0}}}+\sqrt{p^{L}_{i}+p_{w_{0}}}}\leq\frac{1}{2\sqrt{2c}}|p^{U}_{i}-p^{L}_{i}|,$

where we used $p_{w_{0}}\geq c$ specified in (C2) in the last inequality.

Then

$\left\|\sqrt{\frac{p^{U}_{i}+p_{w_{0}}}{2}}-\sqrt{\frac{p^{L}_{i}+p_{w_{0}}}{2}}\right\|_{2}\leq\frac{1}{2\sqrt{2c}}\|p^{U}_{i}-p^{L}_{i}\|_{2}\leq\frac{1}{2\sqrt{2c}}\epsilon.$

Therefore,

$N_{[]}(\frac{1}{2\sqrt{2c}}\epsilon,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2})\leq
N_{[]}(\epsilon,\mathcal{F},\|\cdot\|_{2}).$

####  A.2.2 Calculating the upper bound

Step 1: Calculate the bracketing entropy of one-dimensional non-increasing
density functions supported on $(0,1)$ with envelop function $F\in L^{2}$.

######  Lemma 1.

Let $\mathcal{M}$ be the space of the one-dimensional non-increasing density
function supported on $(0,1)$ with envelop function $F\in L^{2}$. Then we have

$\log N_{[]}(\delta,\mathcal{M},\|\cdot\|_{2})\leq C\frac{1}{\delta}.$

###### Proof.

Let

$g(x)=\int_{0}^{x}F^{2}(t)dt.$

Then $g(x)$ is continuous, and
$g(0)=0,g(1)\geq\int_{0}^{1}f^{2}(t)dt\geq(\int_{0}^{1}f(t)dt)^{2}=1$, where
$f\in\mathcal{M}$. For any small $\delta>0$, there exists $s$ such that
$g(s)=\delta^{2}/2$. For any $f\in\mathcal{M}$, when $x\geq s$, $f(x)\leq
f(s)\leq F(s)$. So that the space

$\mathcal{M}^{\prime}=\\{f(x)\cdot 1\\{x\geq s\\}:f(x)\in M\\}$

is uniformly bounded by $F(s)$.

It has been shown in (Van de Geer, 2000) that

$\log
N_{[]}(\frac{\delta}{\sqrt{2}},\mathcal{M}^{\prime},\|\cdot\|_{2})\leq\frac{C}{\delta}.$

We can find the $\frac{\delta}{\sqrt{2}}$ brackets
$[f_{1}^{L},f_{1}^{U}],\cdots,[f_{n}^{L},f_{n}^{U}]$ with
$n=\lceil\exp(\frac{C}{\delta})\rceil$, where $f_{i}^{L}$ and $f_{i}^{U}$ are
supported on $[s,1]$, and $\lceil x\rceil$ is the smallest interger that is no
smaller than $x$.

Let $g_{i}^{L}(x)=f_{i}^{L}(x)$, $g_{i}^{U}=F(x)\cdot
1\\{0<x<s\\}+f_{i}^{U}(x)$. For any $f\in\mathcal{M}$, $\exists i$ such that
$f_{i}^{L}(x)\leq f\cdot 1\\{x>s\\}\leq f_{i}^{U}(x)$. So that
$g_{i}^{L}(x)\leq f(x)\leq g_{i}^{U}(x)$.

In addition,

$\displaystyle\int_{0}^{1}|g_{i}^{U}(t)-g_{i}^{L}(t)|^{2}dt$
$\displaystyle=\int_{0}^{s}|g_{i}^{U}(t)-g_{i}^{L}(t)|^{2}dt+\int_{s}^{1}|g_{i}^{U}(t)-g_{i}^{L}(t)|^{2}dt$
$\displaystyle=\int_{0}^{s}F^{2}(t)dt+\int_{s}^{1}|f_{i}^{U}(t)-f_{i}^{L}(t)|^{2}dt$
$\displaystyle\leq\frac{\delta^{2}}{2}+(\frac{\delta}{\sqrt{2}})^{2}$
$\displaystyle=\delta^{2}.$

So $\|g_{i}^{L}-g_{i}^{U}\|_{2}\leq\delta$, and
$[g_{1}^{L},g_{1}^{U}],\cdots,[g_{n}^{L},g_{n}^{U}]$ are $\delta$-brackets of
$\mathcal{M}$.

Therefore we have

$\log N_{[]}(\delta,\mathcal{M},\|\cdot\|_{2})\leq\frac{C}{\delta}.$

∎

Step 2: Bound $\|p_{w^{U}}-p_{w^{L}}\|_{2}.$

Write
$p_{w^{U}}(p_{1},p_{2})=\xi_{00}^{U}+\xi_{01}^{U}f_{2}(p_{2})+\xi_{10}^{U}f_{1}(p_{1})+\xi_{11}^{U}f_{1}^{U}(p_{1})f_{2}^{U}(p_{2})$,
and
$p_{w^{L}}(p_{1},p_{2})=\xi_{00}^{L}+\xi_{01}^{L}f_{2}(p_{2})+\xi_{10}^{L}f_{1}(p_{1})+\xi_{11}^{L}f_{1}^{L}(p_{1})f_{2}^{L}(p_{2})$,
where $[\xi_{ij}^{U},\xi_{ij}^{L}],i,j=0,1$ is one of the $\delta$ bracket of
$\Delta$; $[f_{1}^{U},f_{1}^{L}]$ is one of the $\delta$ bracket of
$\mathcal{G}_{1}$; and $[f_{2}^{U},f_{2}^{L}]$ is one of the $\delta$ bracket
of $\mathcal{G}_{2}$, where $\Delta$, $\mathcal{G}_{1}$, and $\mathcal{G}_{2}$
are defined in (7). So we have

$\delta\geq\max\\{|\xi_{ij}^{U}-\xi_{ij}^{L}|,\|f_{1}^{U}-f_{1}^{L}\|_{2},\|f_{2}^{U}-f_{2}^{L}\|_{2}\\},i,j=0,1.$

Using the fact that

$\int A(x)B(y)dxdy=\int A(x)dx\int B(y)dy,$

we obtain

$\displaystyle\|\xi_{11}^{U}f_{1}^{U}(p_{1})f_{2}^{U}(p_{2})$
$\displaystyle-\xi_{11}^{L}f_{1}^{L}(p_{1})f_{2}^{L}(p_{2})\|_{2}$
$\displaystyle\leq\xi_{11}^{U}\|f_{1}^{U}(p_{1})f_{2}^{U}(p_{2})-f_{1}^{L}(p_{1})f_{2}^{L}(p_{2})\|_{2}+\|f_{1}^{L}(p_{1})f_{2}^{L}(p_{2})\|_{2}\cdot|\xi^{U}_{11}-\xi^{L}_{11}|$
$\displaystyle\leq\|f_{1}^{U}(p_{1})f_{2}^{U}(p_{2})-f_{1}^{L}(p_{1})f_{2}^{L}(p_{2})\|_{2}+C\delta$
$\displaystyle\leq
C\delta+\|f_{1}^{U}(p_{1})(f_{2}^{U}(p_{2})-f_{2}^{L}(p_{2}))\|_{2}+\|f_{2}^{L}(p_{2})(f_{1}^{U}(p_{1})-f_{1}^{L}(p_{1}))\|_{2}$
$\displaystyle\leq
C\delta+\sqrt{C}\|f_{2}^{U}-f_{2}^{L}\|_{2}+\sqrt{C}\|f_{1}^{U}-f_{1}^{L}\|_{2}$
$\displaystyle\leq(C+2\sqrt{C})\delta.$

In the same way, we can further show that

$\displaystyle|\xi_{00}^{U}-\xi_{00}^{L}|$ $\displaystyle\leq\delta,$
$\displaystyle\|\xi_{10}^{U}f_{1}^{U}-\xi_{10}^{L}f_{1}^{L}\|_{2}$
$\displaystyle\leq(\sqrt{C}+1)\delta,$
$\displaystyle\|\xi_{01}^{U}f_{2}^{U}-\xi_{01}^{L}f_{2}^{L}\|_{2}$
$\displaystyle\leq(\sqrt{C}+1)\delta,$
$\displaystyle\|\xi_{11}^{U}f_{1}^{U}f_{2}^{U}-\xi_{11}^{L}f_{1}^{L}f_{2}^{L}\|_{2}$
$\displaystyle\leq(C+2\sqrt{C})\delta.$

Then

$\displaystyle\|p_{w^{U}}-p_{w^{L}}\|_{2}$
$\displaystyle\leq|\xi_{00}^{U}-\xi_{00}^{L}|+\|\xi_{10}^{U}f_{1}^{U}-\xi_{10}^{L}f_{1}^{L}\|_{2}+\|\xi_{01}^{U}f_{2}^{U}-\xi_{01}^{L}f_{2}^{L}\|_{2}+\|\xi_{11}^{U}f_{1}^{U}f_{2}^{U}-\xi_{11}^{L}f_{1}^{L}f_{2}^{L}\|_{2}$
$\displaystyle\leq(C+4\sqrt{C}+3)\cdot\delta.$

Step 3: Construct the pairs of functions with bracketing entropy. We first
construct the pairs in $\Delta,\mathcal{G}_{1},\mathcal{G}_{2}$ separately and
combine them together.

Let $\delta>0$, and $n_{1}=\lceil\frac{1}{\delta}\rceil\leq C_{5}\delta^{-1}$.
For $\xi_{11}$, write
$\xi_{11,i}^{L}=\frac{i-1}{n_{1}},\xi_{11,i}^{U}=\frac{i}{n_{1}}$,
$i=1,\cdots,n_{1}$. Then for any $\xi_{11}\in[0,1]$, there exists $i$ such
that $\xi_{11,i}^{L}\leq\xi_{11}\leq\xi_{11,i}^{U}$ and
$|\xi_{11,i}^{L}-\xi_{11,i}^{U}|\leq\delta$. For $\xi_{00},\xi_{10},\xi_{01}$,
we do the same thing.

By Lemma 1, we can choose the correponding bracketing function
$(f_{1,i}^{L},f_{1,i}^{U})$, $i=1,\cdots,n_{2}$, where
$n_{2}\leq\exp(C_{6}/\delta)$. Then for any $f_{1}\in\mathcal{G}_{1}$ there
exists $i$ such that $f_{1,i}^{L}\leq f_{1}\leq f_{1,i}^{U}$ with
$\|f_{1,i}^{L}-f_{1,i}^{U}\|_{2}\leq\delta$. The bracketing number of
$\mathcal{G}_{2}$ is the same as the bracketing number of $\mathcal{G}_{1}$.

Consequently, we obtain the bracketing entropy of $\Delta,\mathcal{G}_{1}$,and
$\mathcal{G}_{2}$, with

$\log n_{1}\leq C_{5}\log\delta^{-1},$

and

$\log n_{2}\leq C_{6}/\delta.$

Therefore, for any

$p_{w}(x,y)=\xi_{00}+\xi_{10}f_{1}(x)+\xi_{01}f_{2}(y)+\xi_{11}f_{1}(x)f_{2}(y),$

there exist $i_{00},i_{01},i_{10},i_{11},j_{1}$, and $j_{2}$ such that

$\displaystyle\xi_{00,i_{00}}^{L}$
$\displaystyle\leq\xi_{00}\leq\xi_{00,i_{00}}^{U},\quad\xi_{01,i_{01}}^{L}\leq\xi_{01}\leq\xi_{01,i_{01}}^{U},$
$\displaystyle\xi_{10,i_{10}}^{L}$
$\displaystyle\leq\xi_{10}\leq\xi_{00,i_{10}}^{U},\quad\xi_{11,i_{11}}^{L}\leq\xi_{11}\leq\xi_{11,i_{11}}^{U},$
$\displaystyle f_{1,j_{1}}^{L}(x)$ $\displaystyle\leq f_{1}(x)\leq
f_{1,j_{1}}^{U}(x),\quad\mbox{and}\quad f_{2,j_{2}}^{L}(y)\leq f_{2}(y)\leq
f_{2,j_{2}}^{U}(y).$

Let $I=\\{(i_{00},i_{01},i_{10},i_{01},j_{1},j_{2})\\}$, where
$i_{\cdot\cdot}=1,\cdots n_{1}$, and $j_{\cdot}=1,\cdots,n_{2}$, then

$\log N_{[]}(\delta,\mathcal{F},\|\cdot\|_{2})\lesssim\log|I|=\log
n_{1}^{4}n_{2}^{2}\leq C_{7}\delta^{-1}.$

To sum up, we have

$\log N_{[]}(\delta,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2})\lesssim\log
N_{[]}(\delta,\mathcal{F},\|\cdot\|_{2})\lesssim C_{7}\delta^{-1}.$

####  A.2.3 Calculating the lower bound

The lower bound calculation is motivated from Blei et al. (2007) and Gao and
Wellner (2007). These papers focus on the bounded monotone functional space
while we focus on the density functional space. Besides, the functional space
is $\bar{\mathcal{F}}^{1/2}$ instead of $\mathcal{F}$.

To calculate the lower bound of
$N_{[]}(\delta,\bar{\mathcal{F}}^{1/2},\|\cdot\|_{2})$, we need to determine
the packing number. The $\epsilon$ packing number of functional space
$\mathcal{G}$ is defined to be the minimal number of disjoint
$g_{i}\in\mathcal{G}$ such that when $i\neq j$, $d(g_{i},g_{j})>\epsilon$. To
generate $\delta$ packing of $\bar{\mathcal{F}}^{1/2}$, we divide the support
of $f\in\bar{\mathcal{F}}^{1/2}$, $(0,1)$, into three parts. The first part of
$f$ is the constant $3$, the second part of $f$ is a step function constructed
based on Blei et al. (2007) and Gao and Wellner (2007), and the third part of
$f$ is another constant so that $f$ is a density function. The values of the
three parts are set to ensure that $f$ is non-increasing.

Consider the following sequence

$A=\\{\mathbf{a}:\mathbf{a}=(a_{1},a_{2},\cdots,a_{i},\cdots,a_{n})\text{
where }a_{i}\in\\{0,1\\},i=1,\cdots,n\\},$

where $n=\lfloor\frac{1}{\sqrt{6912}\delta}\rfloor$, where $\lfloor x\rfloor$
is the biggest interger that is no greater than $x$. Assume that $\delta$ is
small enough so that $n$ is large enough. Define

$g_{\mathbf{a}}(x)=\begin{cases}3&0\leq x<\frac{1}{12}\\\
\frac{4}{3}\\{\frac{n-k+1}{n}+\frac{a_{k}}{n}\\}+\frac{1}{2}&\frac{k-1}{2n}+\frac{1}{12}\leq
x<\frac{k}{2n}+\frac{1}{12}\\\ t_{\mathbf{a}}&\frac{7}{12}\leq
x<1.\end{cases}$

where $1\leq k\leq n$, and $t_{\mathbf{a}}$ is choosen such that
$g_{\mathbf{a}}(x)$ is a density function. In other words,
$\int_{0}^{1}g_{\mathbf{a}}(x)dx=1$. It can be calculated by

$\displaystyle t_{\mathbf{a}}$
$\displaystyle=(1-\frac{7}{12})^{-1}\left\\{1-\frac{3}{12}-\int_{\frac{1}{12}}^{\frac{7}{12}}g_{\mathbf{a}}(x)dx\right\\}$
$\displaystyle=\frac{12}{5}\\{\frac{3}{4}-\frac{1}{4}-\frac{1}{n}\sum\limits_{k=1}^{n}\frac{2}{3}(\frac{n-k+1}{n}+\frac{a_{k}}{n})\\}$
$\displaystyle=\frac{2}{5}-\frac{4}{5n}-\frac{8}{5n^{2}}\sum\limits_{k=1}^{n}a_{k}.$

Note that
$g_{\mathbf{a}}(\frac{1}{12}+)=\frac{4}{3}(1+\frac{a_{1}}{n})+\frac{1}{2}<3=g_{\mathbf{a}}(\frac{1}{12}-)$
and
$g_{\mathbf{a}}(\frac{7}{12}+)<\frac{2}{5}<\frac{1}{2}<g_{\mathbf{a}}(\frac{7}{12}-)$,
and $g_{\mathbf{a}}(x)$ is nonincreasing in $(\frac{1}{12},\frac{7}{12})$.
This implies that $g_{\mathbf{a}}(x)$ is non-increasing in $(0,1)$.

For $\mathbf{a}\neq\mathbf{a^{\prime}}$, $\mathbf{a},\mathbf{a^{\prime}}\in
A$, define the Hamming distance

$d_{h}(\mathbf{a},\mathbf{a^{\prime}})=\mbox{card}\\{k:a_{k}\neq
a^{\prime}_{k},k=1,\cdots,n\\},$

where card means cardinality. For each $\mathbf{a}\in A$, consider the ball

$B(\mathbf{a})=\\{\mathbf{a^{\prime}}\in
A,d_{h}(\mathbf{a},\mathbf{a^{\prime}})\leq\lfloor n/16\rfloor\\}.$

Using the fact that

$\begin{pmatrix}n\\\
i\end{pmatrix}\leq(\frac{ne}{i})^{i}\quad\text{and}\quad(16e)^{1/16}<\sqrt{2},$

we have

$|B(\mathbf{a})|\leq\sum\limits_{i=0}^{\lfloor n/16\rfloor}\begin{pmatrix}n\\\
i\end{pmatrix}<(1+\frac{n}{16})\begin{pmatrix}n\\\ \lfloor
n/16\rfloor\end{pmatrix}\leq(1+\frac{n}{16})(\frac{ne}{n/16})^{n/16}=(1+\frac{n}{16})\\{(16e)^{1/16}\\}^{n}.$

So that when $n$ is large enough, we have

$|B(\mathbf{a})|\leq(1+\frac{n}{16})\\{(16e)^{1/16}\\}^{n}\leq(\sqrt{2})^{n}=2^{n/2}.$

Each ball contains less than $2^{n/2}$ functions. Each of the $|A|=2^{n}$
functions must be contained in at least one ball so that there are at least
$2^{n/2}$ disjoint balls.

Let $A^{\prime}\subset A$ such that $D=\\{B(\mathbf{a}),\mathbf{a}\in
A^{\prime}\\}$ is a set of disjoint balls. Then we have $|D|\geq 2^{n/2}.$
Randomly choose two $B(\mathbf{a}),B(\mathbf{a}^{\prime})\in D$ with
$\mathbf{a}\neq\mathbf{a}^{\prime}$, then
$d_{h}(\mathbf{a},\mathbf{a}^{\prime})>\lfloor n/16\rfloor$. Let
$f_{\mathbf{a}}(x,y)=g_{\mathbf{a}}(x),f_{\mathbf{a}^{\prime}}(x,y)=g_{\mathbf{a^{\prime}}}(x)$.
Notice that for any non-increasing density function $f(x,y)$ supported in
$[0,1]^{2}$, for any $0<a,b<1,$ we have

$1=\int_{0}^{1}\int_{0}^{1}f(x,y)dxdy\geq\int_{0}^{a}\int_{0}^{b}f(x,y)dxdy\geq\int_{0}^{a}\int_{0}^{b}f(a,b)dxdy=abf(a,b),$

then $f(a,b)\leq\frac{1}{ab}$.

We have

$\displaystyle\int_{0}^{1}$
$\displaystyle\int_{0}^{1}\left|\sqrt{\frac{f_{w_{0}}(x,y)+f_{\mathbf{a}}(x,y)}{2}}-\sqrt{\frac{f_{w_{0}}(x,y)+f_{\mathbf{a}^{\prime}}(x,y)}{2}}\right|^{2}dxdy$
$\displaystyle=\frac{1}{2}\int_{0}^{1}\int_{0}^{1}\left|\frac{|g_{\mathbf{a}}(x)-g_{\mathbf{a}^{\prime}}(x)|}{\sqrt{f_{w_{0}}(x,y)+f_{\mathbf{a}}(x,y)}+\sqrt{f_{w_{0}}(x,y)+f_{a^{\prime}}(x,y)}}\right|^{2}dxdy$
$\displaystyle\geq\frac{1}{2}\int_{0}^{1}\int_{0}^{1}\left|\frac{|g_{\mathbf{a}}(x)-g_{\mathbf{a}^{\prime}}(x)|}{2\sqrt{\frac{2}{xy}}}\right|^{2}dxdy$
$\displaystyle=\frac{1}{16}\int_{0}^{1}y\int_{0}^{1}x|g_{\mathbf{a}}(x)-g_{\mathbf{a}^{\prime}}(x)|^{2}dxdy$
$\displaystyle\geq\frac{1}{32}\int_{\frac{1}{12}}^{\frac{7}{12}}x|g_{\mathbf{a}}(x)-g_{\mathbf{a}^{\prime}}(x)|^{2}dx$
$\displaystyle\geq\frac{1}{384}\int_{\frac{1}{12}}^{\frac{7}{12}}|g_{\mathbf{a}}(x)-g_{\mathbf{a}^{\prime}}(x)|^{2}dx$
$\displaystyle\geq\frac{1}{384}\sum\limits_{k=1}^{n}I\\{a_{k}\neq
a^{\prime}_{k}\\}\frac{1}{2n}(\frac{4}{3n})^{2}$
$\displaystyle\geq\frac{1}{6912}n^{-2}.$

Recall that
$n=\lfloor\frac{1}{\sqrt{6912}\delta}\rfloor\leq\frac{1}{\sqrt{6912}\delta}$,
we obtain

$\int_{0}^{1}\int_{0}^{1}\left|\sqrt{\frac{f_{0}(x,y)+f_{a}(x,y)}{2}}-\sqrt{\frac{f_{0}(x,y)+f_{a^{\prime}}(x,y)}{2}}\right|^{2}dxdy\geq\delta^{2}.$

The number of $\delta$-packing is greater than $2^{n/2}$, or
$2^{1/\sqrt{27648}\delta}$, and the bracketing entropy satisfies

$\log N_{[]}(\epsilon,\bar{\mathcal{F}}^{1/2},\|\cdot\|_{2})\geq
C_{8}/\epsilon.$

###  A.3 Proof of Theorem 2

We use the following Proposition adapted from Theorem 7.4 of Van de Geer
(2000) to show the consistency of $p_{\hat{w}_{n}}$ based on Hellinger
distance.

######  Proposition A.2.

Let

$J_{B}(\delta,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2})=\int_{0}^{\delta}\sqrt{\log
N_{[]}(\epsilon,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2})}d\epsilon,$

where

$\bar{\mathcal{F}}^{1/2}=\\{\sqrt{\frac{p_{w}+p_{w_{0}}}{2}}:w\in\Omega\\}.$

Take $\Psi(\delta)\geq
J_{B}\left(\delta,\bar{\mathcal{F}}^{1/2},\|\cdot\|_{2}\right)$ in such a way
that $\Psi(\delta)/\delta^{2}$ is a non-increasing function of $\delta$. Then
for a universal constant $c$ and for

$\sqrt{n}\delta_{n}^{2}\geq c\Psi(\delta_{n}),$

we have for all $\delta\geq\delta_{n}$

$\mathbf{P}\left(h\left(p_{\hat{w}_{n}},p_{w_{0}}\right)>\delta\right)\leq
c\exp\left(-\frac{n\delta^{2}}{c^{2}}\right),$

where $\hat{w}_{n}$ is defined in (9).

From Proposition 1, we have

$\log N_{[]}(\epsilon,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2})\leq
C_{1}/\epsilon.$

Then, we obtain

$\displaystyle J_{B}(\delta,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2})$
$\displaystyle=\int_{0}^{\delta}\sqrt{\log
N_{[]}(\epsilon,\mathcal{\bar{F}}^{1/2},\|\cdot\|_{2})}d\epsilon$
$\displaystyle\leq\int_{0}^{\delta}\sqrt{C_{1}/\epsilon}d\epsilon$
$\displaystyle=2\sqrt{C_{1}}\delta^{1/2}.$

Define $\Psi(\delta)=2\sqrt{C_{1}}\delta^{1/2}$, and
$\delta_{n}=C_{2}n^{-1/3}$, then we have

$\sqrt{n}\delta_{n}^{2}=C_{2}^{2}n^{-1/6}=\frac{C_{2}^{3/2}}{2\sqrt{C_{1}}}\\{2\sqrt{C_{1}}(C_{2}n^{-1/3})^{1/2}\\}=c\Psi(\delta_{n}),$

where $c=\frac{C_{2}^{3/2}}{2\sqrt{C_{1}}}$. By applying Proposition A.2,
$\exists C_{2},C_{3},C_{4}>0$ such that

$\mathbf{P}\left(h\left(p_{\hat{w}_{n}},p_{w_{0}}\right)>C_{2}n^{-1/3}\right)\leq
C_{3}\exp\left[-C_{4}n^{1/3}\right].$

###  A.4 Proof of Theorem 3

We adopt Lemma 5.2 in (Van de Geer, 1993) to show the consistency of $w$.

######  Proposition A.3.

Let $\Omega\subset\Omega^{*}$, where $\left(\Omega^{*},\tau\right)$ is a
compact metric space. Suppose $w\rightarrow p_{w},w\in\Omega^{*}$ is
$\mu$-almost everywhere continuous and that $w_{0}$ is identifiable.

Then, if $h\left(f_{w_{n}},f_{w_{0}}\right)\rightarrow 0$ for some sequence
$\left\\{w_{n}\right\\}$, where $h$ is the Hellinger distance, we have
$\tau\left(w_{n},w_{0}\right)\rightarrow 0$.

According to Theorem 1, we have shown the identifiability of $w_{0}$. To use
Proposition A.3, we need to find a compact space $\Omega^{*}$, and show that
$w\to p_{w}$ is a continuous mapping.

Define a bigger parameter space
$\Omega^{*}=\Delta\times\mathcal{G}^{*}_{1}\times\mathcal{G}_{2}^{*}$, where

$\mathcal{G}_{1}^{*}=\\{f:f\text{ is non-increasing on }(0,1],\|f\|_{1}\leq
1,0\leq f\leq F\text{ a.e.}\\},$

and

$\mathcal{G}_{2}^{*}=\\{f:f\text{ is non-increasing on }(0,1],\|f\|_{1}\leq
1,0\leq f\leq F\text{ a.e.}\\}.$

The corresponding product metric $\tau$ is

$\displaystyle\tau(w,\hat{w})$
$\displaystyle=|\xi_{00}-\hat{\xi}_{00}|+|\xi_{01}-\hat{\xi}_{01}|+|\xi_{10}-\hat{\xi}_{10}|+|\xi_{11}-\hat{\xi}_{11}|$
$\displaystyle\qquad+\int_{0}^{1}|f_{1}(x)-\hat{f}_{1}(x)|dx+\int_{0}^{1}|f_{2}(y)-\hat{f}_{2}(y)|dy.$

The following two lemmas will show that $w\to p_{w}$ is continuous, and
$\Omega^{*}$ is compact.

######  Lemma 2.

Let $w=(\xi_{00},\xi_{01},\xi_{10},\xi_{11},f_{1}(x),f_{2}(y))$. The mapping
$p:\Omega^{*}\to\mathcal{F}$, where

$p(w)(x,y)=p_{w}(x,y)=\xi_{00}+\xi_{10}f_{1}(x)+\xi_{01}f_{2}(y)+\xi_{11}f_{1}(x)f_{2}(y),$

is a continuous mapping, where we use the product metric $\tau$ in
$\Omega^{*}$ and the $L^{1}$ metric in $\mathcal{F}$.

######  Lemma 3.

$\Omega^{*}$ is a compact space.

Identifiability Theorem 1, continuous mapping Lemma 2, and compact space Lemma
3 provide all conditions needed in Proposition A.3. So we have
$\tau(\hat{w}_{n},w_{0})\overset{p}{\to}0$.

####  A.4.1 Proof of Lemma 2

Now we show that $w\to p_{w},w\in\Omega^{*}$ is continuous. Assume that
$\tau(w^{*},w)\leq\delta$, using the fact that

$\int A(x)B(y)dxdy=\int A(x)dx\int B(y)dy,$

we obtain

$\displaystyle\|\xi_{11}f_{1}(p_{1})f_{2}(p_{2})$
$\displaystyle-\xi_{11}^{*}f_{1}^{*}(p_{1})f_{2}^{*}(p_{2})\|_{1}$
$\displaystyle\leq\xi_{11}\|f_{1}(p_{1})f_{2}(p_{2})-f_{1}^{*}(p_{1})f_{2}^{*}(p_{2})\|_{1}+\|f_{1}^{*}(p_{1})f_{2}^{*}(p_{2})\|_{1}\cdot|\xi_{11}-\xi_{11}^{*}|$
$\displaystyle\leq\|f_{1}(p_{1})f_{2}(p_{2})-f_{1}^{*}(p_{1})f_{2}^{*}(p_{2})\|_{1}+\delta$
$\displaystyle\leq\delta+\|f_{1}(p_{1})(f_{2}(p_{2})-f_{2}^{*}(p_{2}))\|_{1}+\|f_{2}^{*}(p_{2})(f_{1}(p_{1})-f_{1}^{*}(p_{1}))\|_{1}$
$\displaystyle\leq\delta+\|f_{2}-f_{2}^{*}\|_{1}+\|f_{1}-f_{1}^{*}\|_{1}$
$\displaystyle\leq 3\delta,$

In the same way, we can show that

$\displaystyle|\xi_{00}-\xi_{00}^{*}|$ $\displaystyle\leq\delta$
$\displaystyle\|\xi_{10}f_{1}-\xi_{10}^{*}f_{1}^{*}\|_{1}$ $\displaystyle\leq
2\delta$ $\displaystyle\|\xi_{01}f_{2}-\xi_{01}^{*}f_{2}^{*}\|_{1}$
$\displaystyle\leq 2\delta$
$\displaystyle\|\xi_{11}f_{1}(p_{1})f_{2}(p_{2})-\xi_{11}^{*}f_{1}^{*}(p_{1})f_{2}^{*}(p_{2})\|_{1}$
$\displaystyle\leq 3\delta$ (A.9)

Therefore

$\displaystyle\|p_{w}-p_{w^{*}}\|_{1}$
$\displaystyle\leq|\xi_{00}-\xi_{00}^{*}|+\|\xi_{10}f_{1}-\xi_{10}^{*}f_{1}^{*}\|_{1}+\|\xi_{01}f_{2}-\xi_{01}^{*}f_{2}^{*}\|_{1}$
$\displaystyle\qquad+\|\xi_{11}f_{1}(p_{1})f_{2}(p_{2})-\xi_{11}^{*}f_{1}^{*}(p_{1})f_{2}^{*}(p_{2})\|$
$\displaystyle\leq 8\delta,$

implying that $w\to p_{w}$ is continuous.

####  A.4.2 Proof of Lemma 3

$\mathcal{G}_{1}^{*}$ is totally bounded when $\mathcal{G}_{1}^{*}$ can be
covered by a finite number of $\epsilon$-balls. $\mathcal{G}_{1}^{*}$ is
complete when every Cauchy sequence in it converges.

Step 1: $\mathcal{G}_{1}^{*}$ is totally bounded.

The proof is an application of Kolmogorov?Riesz Theorem Hanche-Olsen and
Holden (2010). The Theorem says that $G_{1}^{*}$ is totally bounded with
metric $L^{1}$ if and only if

  1. 1.

$\mathcal{G}_{1}^{*}$ is bounded,

  2. 2.

for any $\epsilon>0$, there exists $R,$ and $\delta$, such that for any
$g\in\mathcal{G}^{*}_{1}$, $|y|<\delta$, we have

$\int_{|x|>R}f(x)<\epsilon,\quad\text{and}\quad\int_{-\infty}^{\infty}|f(x+y)-f(x)|dx<\epsilon.$

######  Remark 3.

$\mathcal{G}_{1}^{*}$ being bounded means that for any
$g,g^{\prime}\in\mathcal{G}_{1}^{*}$, $\int|g(x)-g^{\prime}(x)|dx<C$, where
$C$ is a constant. It is not the same as $\forall g\in\mathcal{G}_{1}^{*},g$
is bounded.

The remaining part is to verify the above conditions.

  1. 1.

$\mathcal{G}_{1}^{*}$ is bounded because for any $g\in\mathcal{G}_{1}^{*}$, we
have $\|g-0\|_{1}\leq 1$, so that for any
$g,g^{\prime}\in\mathcal{G}_{1}^{*}$,
$\|g-g^{\prime}\|_{1}\leq\|g-0\|_{1}+\|g^{\prime}-0\|_{1}\leq 2$.

  2. 2.

Taking $R=1$ will guarantee that $\int_{|x|>R}g(x)dx=0<\epsilon$.

  3. 3.

Take $\delta$ such that $2\int_{0}^{\delta}F(x)\leq\epsilon$. It is achievable
because $F\in L^{1}$. Without loss generality, we can assume that $y>0$, so
that

$\displaystyle\int_{-\infty}^{\infty}|g(x+y)-g(x)|dx$
$\displaystyle=\int_{-y}^{0}|g(x+y)-g(x)|dx+\int_{0}^{\infty}|g(x+y)-g(x)|dx$
$\displaystyle=\int_{-y}^{0}g(x+y)dx+\int_{0}^{\infty}\\{g(x)-g(x+y)\\}dx$
$\displaystyle=2\int_{0}^{y}g(x)dx$ $\displaystyle\leq 2\int_{0}^{y}F(x)dx$
$\displaystyle\leq 2\int_{0}^{\delta}F(x)dx$ $\displaystyle\leq\epsilon.$

After verifying above conditions, by Kolmogorov?Riesz Theorem,
$\mathcal{G}_{1}^{*}$ is totally bouned.

Step 2: $\mathcal{G}_{1}^{*}$ is complete.

If $\\{f_{n}\\}$ is a Cauchy sequence in $\mathcal{G}_{1}^{*}\subset L^{1}$,
and $f_{n}\in\mathcal{G}_{1}^{*}$, then there exists $f\in L^{1}$ such that
$\|f_{n}-f\|_{1}\to 0$, because of the completeness of $L^{1}$ space. By
Corollary 2.3 in (Stein and Shakarchi, 2009), there exists a subsequence of
$f_{n_{k}}$ such that $f_{n_{k}}(x)\to f(x)$ for almost everywhere. As $0\leq
f_{n}\leq F$, $f_{n}$ is non-increasing and $\|f_{n}\|\leq 1$, we have $f\leq
F$, $f$ is non-increasing and $0\leq f\leq F$ almost everywhere, which implies
that $f\in\mathcal{G}_{1}^{*}$. So that $\mathcal{G}_{1}^{*}$ is complete.

By showing $\mathcal{G}_{1}^{*}$ is totally bounded and complete, we prove
that $\mathcal{G}_{1}^{*}$ is compact. The proof of compactness of
$\mathcal{G}_{2}^{*}$ is exactly the same as that of $\mathcal{G}_{1}^{*},$ so
we omit the details.

###  A.5 Proof of Theorem 4

We first prove that when using the MLE, $\hat{w}_{n}$, to estimate the local
false discovery rate (Lfdr), the average absolute difference between estimated
Lfdr and true Lfdr converges to $0$ based on the following proposition.

######  Proposition A.4.

Under (C1) and (C2), if
$\hat{w}_{n}=(\hat{\xi}_{00},\hat{\xi}_{01},\hat{\xi}_{10},\hat{\xi}_{11},\hat{f}_{1},\hat{f}_{2})$
is the MLE, and $w_{0}=(\xi_{00},\xi_{01},\xi_{10},\xi_{11},f_{1},f_{2})$ is
the true parameter, we have

$\frac{1}{n}\sum\limits_{i=1}^{n}\left|\widehat{\mbox{Lfdr}}(p_{1i},p_{2i})-\mbox{Lfdr}(p_{1i},p_{2i})\right|\overset{p}{\to}0,$

where

$\widehat{\mbox{Lfdr}}(p_{1i},p_{2i})=\frac{\hat{\xi}_{00}+\hat{\xi}_{01}\hat{f}_{2}\left(p_{2i}\right)+\hat{\xi}_{10}\hat{f}_{1}\left(p_{1i}\right)}{\hat{\xi}_{00}+\hat{\xi}_{01}\hat{f}_{2}\left(p_{2i}\right)+\hat{\xi}_{10}\hat{f}_{1}\left(p_{1i}\right)+\hat{\xi}_{11}\hat{f}_{1}\left(p_{1i}\right)\hat{f}_{2}\left(p_{2i}\right)},$

and

$\mbox{Lfdr}(p_{1i},p_{2i})=\frac{\xi_{00}+\xi_{01}f_{2}\left(p_{2i}\right)+\xi_{10}f_{1}\left(p_{1i}\right)}{\xi_{00}+\xi_{01}f_{2}\left(p_{2i}\right)+\xi_{10}f_{1}\left(p_{1i}\right)+\xi_{11}f_{1}\left(p_{1i}\right)f_{2}\left(p_{2i}\right)}.$

####  A.5.1 Proof of Proposition A.4

Let
$\mathcal{H}=\\{\frac{\xi_{00}+\xi_{10}f_{1}(p_{1})+\xi_{01}f_{2}p_{2})}{\xi_{00}+\xi_{10}f_{1}(p_{1})+\xi_{01}f_{2}(p_{2})+\xi_{11}f_{1}(p_{1})f_{2}(p_{2})}:w\in\Omega\\}$.
For $w,w^{*}\in\Omega$ and $\tau(w,w^{*})<\delta$, notice that $\xi_{00}\geq
c,\xi_{00}^{*}\geq c,p_{w}\geq c,$ and $p_{w^{*}}\geq c$. We have

$\displaystyle\left\|\frac{p_{w}-\xi_{11}f_{1}f_{2}}{p_{w}}-\frac{p_{w^{*}}-\xi^{*}_{11}f^{*}_{1}f^{*}_{2}}{p_{w^{*}}}\right\|_{1}$
$\displaystyle=\left\|\frac{\xi_{11}f_{1}f_{2}p_{w^{*}}-\xi^{*}_{11}f_{1}^{*}f_{2}^{*}p_{w}}{p_{w}p_{w^{*}}}\right\|_{1}$
$\displaystyle\leq\left\|\frac{\xi_{11}f_{1}f_{2}|p_{w}-p_{w^{*}}|+p_{w^{*}}|\xi_{11}f_{1}f_{2}-\xi_{11}^{*}f_{1}^{*}f_{2}^{*}|}{p_{w}p_{w^{*}}}\right\|_{1}$
$\displaystyle\leq\left\|\frac{\xi_{11}f_{1}f_{2}}{p_{w}}\frac{|p_{w}-p_{w^{*}}|}{p_{w^{*}}}\right\|_{1}+\left\|\frac{1}{p_{w}}|\xi_{11}f_{1}f_{2}-\xi_{11}^{*}f_{1}^{*}f_{2}^{*}\right\|_{1}$
$\displaystyle\leq\frac{\|p_{w}-p_{w^{*}}\|_{1}+\|\xi_{11}f_{1}f_{2}-\xi_{11}^{*}f_{1}^{*}f_{2}^{*}\|_{1}}{c}$
$\displaystyle\leq\frac{11}{c}\delta.$

The last inequality holds because of (A.9). Then
$\mbox{Lfdr}(w):\Omega\to\mathcal{H}$ is a continuous function. By continuous
mapping theorem (with corresponding random element), we have

$\|\widehat{\mbox{Lfdr}}-\mbox{Lfdr}\|_{1}\overset{p}{\to}0.$

In addition, by Holder’s inequality, we have

$\displaystyle\int|\widehat{\mbox{Lfdr}}(p_{1i},p_{2i})$
$\displaystyle-\mbox{Lfdr}(p_{1i},p_{2i})|f_{w_{0}}(p_{1i},p_{2i})dp_{1i}dp_{2i}$
$\displaystyle\leq\|f_{w_{0}}\|_{2}\left\\{\int\left|\widehat{\mbox{Lfdr}}(p_{1i},p_{2i})-\mbox{Lfdr}(p_{1i},p_{2i})\right|^{2}dp_{1i}dp_{2i}\right\\}^{1/2}$
$\displaystyle\leq\|f_{w_{0}}\|_{2}\left\\{\int
2\left|\widehat{\mbox{Lfdr}}(p_{1i},p_{2i})-\mbox{Lfdr}(p_{1i},p_{2i})\right|dp_{1i}dp_{2i}\right\\}^{1/2}$
$\displaystyle\overset{p}{\to}0.$

So that $\mathbb{E}|\widehat{\mbox{Lfdr}}-\mbox{Lfdr}|\overset{p}{\to}0$.

Notice that

$\frac{\xi_{00}+\xi_{10}f_{1}(p_{1})+\xi_{01}f_{2}(p_{2})}{\xi_{00}+\xi_{10}f_{1}(p_{1})+\xi_{01}f_{2}(p_{2})+\xi_{11}f_{1}(p_{1})f_{2}(p_{2})}=1-\frac{\xi_{11}}{\frac{\xi_{00}}{f_{1}(p_{1})f_{2}(p_{2})}+\frac{\xi_{10}}{f_{2}(p_{2})}+\frac{\xi_{01}}{f_{1}(p_{1})}+\xi_{11}}$

is a monotone increasing function in $p_{1},p_{2}$, and bounded by 1. By Gao
and Wellner (2007) the bracketing number of $\mathcal{H}$ is bounded, and we
have

$N_{[]}(\mathcal{H},\epsilon,\|\cdot\|_{2})<\infty.$

Define $\mathcal{H}^{\prime}=\\{|h-\mbox{Lfdr}(w_{0})|:h\in\mathcal{H}\\}$. We
aim to show that

$N_{[]}(\mathcal{H}^{\prime},\epsilon,\|\cdot\|_{2})<\infty,$

which can be achieved by the following lemma.

######  Lemma 4.

Let $\mathcal{F}$ be a functional space, and let $\alpha\geq 1$. Define

$|\mathcal{F}|=\\{|f|,f\in\mathcal{F}\\}.$

If

$N_{[]}(\mathcal{F},\epsilon,\|\cdot\|_{\alpha})<\infty,$

we have

$N_{[]}(|\mathcal{F}|,\epsilon,\|\cdot\|_{\alpha})\leq
N_{[]}(\mathcal{F},\epsilon,\|\cdot\|_{\alpha}).$

###### Proof.

Let $[f_{i}^{L},f_{i}^{U}],i=1,\cdots,n$ be $\epsilon$ bracketing of
$\mathcal{F}$. Define $g^{U}_{i},g^{L}_{i}$ in the following way:

  * •

When $f_{i}^{U}(x)\geq f_{i}^{L}(x)>0$, $g^{U}_{i}(x)=f^{U}_{i}(x)$,
$g^{L}_{i}(x)=f^{L}_{i}(x)$;

  * •

When $f_{i}^{U}(x)\geq 0\geq
f_{i}^{L}(x),g^{U}_{i}(x)=\max\\{|f^{U}_{i}(x)|,|f^{L}_{i}(x)|\\}$,
$g^{L}_{i}(x)=0$;

  * •

When $0>f_{i}^{U}(x)\geq f_{i}^{L}(x)$, $g^{U}_{i}(x)=-f^{L}_{i}(x)$,
$g^{L}_{i}(x)=-f^{U}_{i}(x)$.

In all the three cases, we have
$|g_{i}^{U}(x)-g_{i}^{L}(x)|\leq|f_{i}^{U}(x)-f_{i}^{L}(x)|$. When
$f_{i}^{L}\leq f\leq f_{i}^{U}$, we also have

$g_{i}^{L}\leq|f|\leq g_{i}^{U}.$

This implies that $[g_{i}^{L},g_{i}^{U}],i=1,\cdots,n$ is also $\epsilon$
bracketing of functional space $|\mathcal{F}|$. Therefore we have

$N_{[]}(|\mathcal{F}|,\epsilon,\|\cdot\|_{\alpha})\leq
N_{[]}(\mathcal{F},\epsilon,\|\cdot\|_{\alpha}).$

∎

By Holder’s inequality, we have

$\|g\|_{L_{1}(P_{0})}=\|gf_{w_{0}}\|_{1}\leq\|g\|_{2}\|f_{w_{0}}\|_{2},$

so that the both the space $\mathcal{H}$, and $\mathcal{H}^{\prime}$ have
bounded bracketing number in $L_{1}(P_{0})$ metric. This satisfies the
condition of the Uniformly Law of Large Numbers (Lemma 3.1 in (Van de Geer,
2000)). So that

$P_{n}|\widehat{\mbox{Lfdr}}-\mbox{Lfdr}|\overset{p}{\to}E|\widehat{\mbox{Lfdr}}-\mbox{Lfdr}|\overset{p}{\to}0.$

####  A.5.2 Proof of Theorem 4

We present the sketch of the proof as follows. Define

$\displaystyle R_{m}(\lambda)$
$\displaystyle=\sum_{i=1}^{m}\mathbf{I}\left\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\right\\},$
$\displaystyle V_{m}(\lambda)$
$\displaystyle=\sum\limits_{i=1}^{m}\mathbf{I}\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\\}(1-\theta_{1i}\theta_{2i}),$
$\displaystyle D_{m}(\lambda)$
$\displaystyle=\sum_{i=1}^{m}\operatorname{Lfdr}(p_{1i},p_{2i})\mathbf{I}\left\\{\operatorname{Lfdr}(p_{1i},p_{2i})\leq\lambda\right\\},$

where $R_{m}$ is the number of rejections, $V_{m}$ is the number of false
discoveries and $D_{m}$ is the estimated number of false discoveries. If we
use the estimated local false discovery rate, the corresponding values are
$\hat{R}_{m}(\lambda),\hat{V}_{m}(\lambda),$ and $\hat{D}_{m}(\lambda)$, where

$\displaystyle\hat{R}_{m}(\lambda)$
$\displaystyle=\sum_{i=1}^{m}\mathbf{I}\left\\{\widehat{\operatorname{Lfdr}}(p_{1i},p_{2i})\leq\lambda\right\\},$
$\displaystyle\hat{V}_{m}(\lambda)$
$\displaystyle=\sum\limits_{i=1}^{m}\mathbf{I}\\{\widehat{\operatorname{Lfdr}}(p_{1i},p_{2i})\leq\lambda\\}(1-\theta_{1i}\theta_{2i}),$
$\displaystyle\hat{D}_{m}(\lambda)$
$\displaystyle=\sum_{i=1}^{m}\widehat{\operatorname{Lfdr}}(p_{1i},p_{2i})\mathbf{I}\left\\{\widehat{\operatorname{Lfdr}}(p_{1i},p_{2i})\leq\lambda\right\\}.$

With a good estimate of the local false discovery rate as that in Proposition
A.4, that satisfies

$\frac{1}{m}\sum_{i=1}^{m}\left|\widehat{\operatorname{Lfdr}}(p_{1i},p_{2i})-\operatorname{Lfdr}(p_{1i},p_{2i})\right|\overset{p}{\to}0,$

we have the following uniform consistency result:

$\displaystyle\sup\limits_{\lambda\in[\lambda_{0},1]}|\hat{R}_{m}(\lambda)-R_{m}(\lambda)|$
$\displaystyle\overset{p}{\to}0,$
$\displaystyle\sup\limits_{\lambda\in[\lambda_{0},1]}|\hat{D}_{m}(\lambda)-D_{m}(\lambda)|$
$\displaystyle\overset{p}{\to}0,$
$\displaystyle\sup\limits_{\lambda\in[\lambda_{0},1]}|\hat{V}_{m}(\lambda)-V_{m}(\lambda)|$
$\displaystyle\overset{p}{\to}0,$

where $\lambda_{0}$ is defined in (C5).

Using the similar way for proving “Glivenko?Cantelli Theorem” ((Resnick,
2005)), we show that $\hat{R}_{m}(\lambda)/m,$ $\hat{D}_{m}(\lambda)/m,$ and
$\hat{V}_{m}(\lambda)/m$ are uniformly consistent. In other words,

$\displaystyle\sup\limits_{\lambda\in[\lambda_{0},1]}|\hat{R}_{m}(\lambda)/m-B_{1}(\lambda)|$
$\displaystyle\overset{p}{\to}0,$
$\displaystyle\sup\limits_{\lambda\in[\lambda_{0},1]}|\hat{D}_{m}(\lambda)/m-B_{2}(\lambda)|$
$\displaystyle\overset{p}{\to}0,$
$\displaystyle\sup\limits_{\lambda\in[\lambda_{0},1]}|\hat{V}_{m}(\lambda)/m-B_{2}(\lambda)|$
$\displaystyle\overset{p}{\to}0.$

Besides, the estimated false discovery propotion
$\frac{\hat{D}_{m}(\lambda)}{\hat{R}_{m}(\lambda)}$, and oracle false
discovery proportion $\frac{V_{m}(\lambda)}{R_{m}(\lambda)}$ uniformly
converge to $\frac{B_{2}(\lambda)}{B_{1}(\lambda)}$, when
$\lambda\in[\lambda_{0},1]$, and

$\displaystyle\sup\limits_{\lambda\in[\lambda_{0},1]}\left|\frac{\hat{D}_{m}(\lambda)}{\hat{R}_{m}(\lambda)}-\frac{B_{2}(\lambda)}{B_{1}(\lambda)}\right|$
$\displaystyle\overset{p}{\to}0,$
$\displaystyle\sup\limits_{\lambda\in[\lambda_{0},1]}\left|\frac{V_{m}(\lambda)}{R_{m}(\lambda)}-\frac{B_{2}(\lambda)}{B_{1}(\lambda)}\right|$
$\displaystyle\overset{p}{\to}0.$

So when $m$ is sufficiently large we have

$P\left\\{\left|\frac{\hat{D}_{m}(\lambda_{0})}{\hat{R}_{m}(\lambda_{0})}-\frac{B_{2}(\lambda_{0})}{B_{1}(\lambda_{0})}\right|\leq\frac{\alpha-
B_{2}(\lambda_{0})/B_{1}(\lambda_{0})}{2}\right\\}\to 1.$

So that

$P\left\\{\frac{\hat{D}_{m}(\lambda_{0})}{\hat{R}_{m}(\lambda_{0})}<\alpha\right\\}\to
1,\quad{and}\quad P\\{\hat{\lambda}_{m}>\lambda_{0}\\}\to 1.$

Besides $P(R_{m}(\hat{\lambda}_{m})>0)\to 1$, so the false discovery rate can
be bounded by

$\displaystyle\mbox{FDR}(\hat{\lambda}_{m})$
$\displaystyle=\mathbb{E}\left(\frac{V_{m}(\hat{\lambda}_{m})}{R_{m}(\hat{\lambda}_{m})\vee
1}\right)$
$\displaystyle=\mathbb{E}\left(\frac{V_{m}(\hat{\lambda}_{m})}{R_{m}(\hat{\lambda}_{m})}\right)$
$\displaystyle=\mathbb{E}\left\\{\frac{\hat{D}_{m}(\hat{\lambda}_{m})}{\hat{R}_{m}(\hat{\lambda}_{m})}+\left(\frac{V_{m}(\hat{\lambda}_{m})}{R_{m}(\hat{\lambda}_{m})}-\frac{B_{2}(\hat{\lambda}_{m})}{B_{1}(\hat{\lambda}_{m})}\right)+\left(\frac{B_{2}(\hat{\lambda}_{m})}{B_{1}(\hat{\lambda}_{m})}-\frac{\hat{D}_{m}(\hat{\lambda}_{m})}{\hat{R}_{m}(\hat{\lambda}_{m})}\right)\right\\}$
$\displaystyle\leq\alpha+\mathbb{E}\left|\frac{V_{m}(\hat{\lambda}_{m})}{R_{m}(\hat{\lambda}_{m})}-\frac{B_{2}(\hat{\lambda}_{m})}{B_{1}(\hat{\lambda}_{m})}\right|+\mathbb{E}\left|\frac{\hat{D}_{m}(\hat{\lambda}_{m})}{\hat{R}_{m}(\hat{\lambda}_{m})}-\frac{B_{2}(\hat{\lambda}_{m})}{B_{1}(\hat{\lambda}_{m})}\right|$
$\displaystyle\leq\alpha+\mathbb{E}\sup\limits_{\lambda\in[\lambda_{0},1]}\left|\frac{V_{m}(\lambda)}{R_{m}(\lambda)}-\frac{B_{2}(\lambda)}{B_{1}(\lambda)}\right|+\mathbb{E}\sup\limits_{\lambda\in[\lambda_{0},1]}\left|\frac{\hat{D}_{m}(\lambda)}{\hat{R}_{m}(\lambda)}-\frac{B_{2}(\lambda)}{B_{1}(\lambda)}\right|$

Based on the following lemma, we prove the result that

$\limsup\limits_{m\to\infty}\mbox{FDR}(\hat{\lambda}_{m})\leq\alpha.$

######  Lemma 5.

If $|X_{m}|\leq C$, and $X_{m}\overset{p}{\to}0$, then we have

$\limsup\limits_{m\to\infty}\mathbb{E}[X_{m}]\leq 0.$

###### Proof.

We take a subsequence of $X_{m}$ such that
$\limsup_{m}\mathbb{E}[X_{m}]=\lim_{k}\mathbb{E}[X_{m_{k}}]$. As
$X_{m}\overset{p}{\to}0$, we can further select a subsequence such that
$X_{m_{k_{l}}}\overset{a.s.}{\to}0$ when $l\to\infty$. (see p. 172 of Resnick
(2005))

Then by Fatou’s Lemma, because $|X_{m}|\leq C$, we have

$\limsup_{m}\mathbb{E}[X_{m}]=\limsup_{k}\mathbb{E}[X_{m_{k}}]\leq\mathbb{E}[\limsup_{l}X_{m_{k_{l}}}]=0.$

∎

##  Appendix B Competing methods

In the simulation studies, we compared the proposed method to several other
methods for replicability analysis, including ad hoc BH, MaxP (Benjamini et
al., 2009), IDR (Li et al., 2011), MaRR (Philtron et al., 2018), radjust
(Bogomolov and Heller, 2018) and JUMP (Lyu et al., 2023). Let
$(p_{1i},p_{2i}),i=1,\dots,m$ denote the paired $p$-values from two studies.
We review these competing methods as follows.

###  B.1 The ad hoc BH method

BH (Benjamini and Hochberg, 1995) is the most popular multiple testing
procedure that conservatively controls the FDR for $m$ independent or
positively correlated tests. In study $j,j=1,2$, the BH procedure proceeds as
below:

  * •

Step 1. Let $p_{j(1)}\leq p_{j(2)}\leq\dots\leq p_{j(m)}$ be the ordered
$p$-values in study $j$, and denote by $H^{(j)}_{(i)}$ the null hypothesis
corresponding to $p_{j(i)}$;

  * •

Step 2. Find the largest $i$ such that $p_{j(i)}\leq\frac{i}{m}\alpha$, i.e.,
$\hat{k}=\mbox{max}\\{i\geq 1:p_{j(i)}\leq\frac{i}{m}\alpha\\}$, and
$\hat{k}=0$ if the set is empty;

  * •

Step 3. Reject all $H^{(j)}_{(i)}$ for $i=1,\dots,\hat{k}$.

The ad hoc BH method for replicability analysis identifies features rejected
by both studies as replicable signals.

###  B.2 The MaxP method

Define the maximum $p$-values as

$q_{i}=\mbox{max}\\{p_{1i},p_{2i}\\},i=1,\dots,m.$

The naïve MaxP method directly applies BH (Benjamini and Hochberg, 1995) to
$q_{i},i=1,\dots,m$ for FDR control of replicability analysis.

###  B.3 The IDR procedure

The IDR procedure Li et al. (2011) deals with high throughput experimental
data from two studies. For feature $i,$ we have the bivariate observations
$(x_{1i},x_{2i}),i=1,\dots,m$. It is assumed that
$(x_{1i},x_{2i}),i=1,\dots,m$ consist of genuine signals (replicable signals
across two studies) and spurious signals (non-replicable signals). Let $K_{i}$
denote whether the $i$th feature is a replicable signal ($K_{i}=1$) or not
($K_{i}=0$). It is assumed that $K_{i},i=1,\ldots,m$ are independent and
follow the Bernoulli distribution. Denote $\pi_{1}=P(K_{i}=1)$. To induce
dependence between $x_{1i}$ and $x_{2i}$, we use a copula model. Specifically,
we assume that the observed data $(x_{1i},x_{2i})$ are generated from latent
variables $(z_{1i},z_{2i})$. The latent variables

$\displaystyle\begin{pmatrix}z_{1i}\\\ z_{2i}\end{pmatrix}\Bigg{|}K_{i}=k$
$\displaystyle\sim N\begin{pmatrix}\begin{pmatrix}\mu_{k}\\\
\mu_{k}\end{pmatrix},\sigma_{k}^{2}\begin{pmatrix}1&\rho_{k}\\\
\rho_{k}&1\end{pmatrix}\end{pmatrix},\text{ }k=0,1,$

where $\mu_{0}=0,\mu_{1}>0,\sigma_{0}^{2}=1,\sigma_{1}^{2}>0,\rho_{0}=0,$ and
$0<\rho_{1}\leq 1$. The cdf of $z_{ji}$ is

$G(x)=P(z_{ji}\leq
x)=\pi_{1}\Phi\left(\frac{x-\mu_{1}}{\sigma_{1}}\right)+(1-\pi_{1})\Phi(x).$

Denote the marginal distribution function of $x_{ji},i=1,\ldots,m;j=1,2,$ as
$F_{j}.$ Generate

$x_{ji}=F_{j}^{-1}(G(z_{ji})),i=1,\dots,m;j=1,2.$

In this way, dependence across two studies is produced. To control the false
discovery rate, we use the local irreproducible discovery rate (idr) as the
test statistic, which is defined as the posterior probability of $K_{i}=0$
given $(x_{1i},x_{2i}).$ Specifically,

$\displaystyle idr(x_{1i}$ $\displaystyle,x_{2i}):=P(K_{i}=0\mid
x_{1i},x_{2i})$
$\displaystyle=\frac{(1-\pi_{1})h_{0}[G^{-1}\\{F_{1}(x_{1i})\\},G^{-1}\\{F_{2}(x_{2i})\\}]}{(1-\pi_{1})h_{0}[G^{-1}\\{F_{1}(x_{1i})\\},G^{-1}\\{F_{2}(x_{2i})\\}]+\pi_{1}h_{1}[G^{-1}\\{F_{1}(x_{1i})\\},G^{-1}\\{F_{2}(x_{2i})\\}]}.$

where

$h_{k}\sim N\begin{pmatrix}\begin{pmatrix}\mu_{k}\\\
\mu_{k}\end{pmatrix},\sigma_{k}^{2}\begin{pmatrix}1&\rho_{k}\\\
\rho_{k}&1\end{pmatrix}\end{pmatrix},\ k=0,1.$

The estimation of $(\pi_{1},\mu_{1},\sigma_{1}^{2},\rho_{1})$ and
$(F_{1},F_{2})$ is through the EM algorithm (Dempster et al., 1977). The step-
up procedure based on ordered idr can be used for FDR control (Sun and Cai,
2007). Specifically, let $idr_{(1)}\leq\dots\leq idr_{(m)}$ be the ranked
$idr$ values, and denote $H_{(1)},\dots,H_{(m)}$ as the corresponding
hypotheses. Find $l=\mbox{max}\\{i:i^{-1}\sum_{j=1}^{i}idr_{j}\leq\alpha\\}$,
and reject all $H_{(i)}$ with $i=1,\dots,l$.

###  B.4 The MaRR procedure

The MaRR procedure (Philtron et al., 2018) uses the maximum rank of each
feature. The null hypothesis is that $H_{i0}:p_{1i}\text{ and }p_{2i}$ are
irreproducible. Denote $(R_{1i},R_{2i})$ as the ranks of
$(p_{1i},p_{2i}),i=1,\dots,m$ within each study. Define

$M_{i}=\mbox{max}\\{R_{1i},R_{2i}\\},i=1,\dots,m.$

Let $\pi_{1}$ denote the proportion of replicable signals. Under the
assumptions:
(I1) if gene $g$ is reproducible and gene $h$ is irreproducible

$R_{1g}<R_{1h},\quad R_{2g}<R_{2h};$

(I2) the correlation between the ranks of the reproducible gene is non-
negative;
(I3) the two ranks of the irreproducible gene are independent,
irreproducible ranks $R_{1i}$ and $R_{2i}$ are uniformly distributed between
$\lfloor m\pi_{1}\rfloor+1$ and $m$. Denote the conditional null survival
function of $M_{i}/m$ as

$\displaystyle S_{m,\pi_{1}}(x)=P(M_{i}/m>x\mid\text{ gene }i\text{ is
irreproducible})$ $\displaystyle=$ $\displaystyle 1-P\left(R_{1i}/m\leq
x,R_{2i}/m\leq x\mid\text{ gene }i\text{ is irreproducible}\right)$
$\displaystyle=$ $\displaystyle 1-\prod_{j=1}^{2}P\left(R_{ji}/m\leq
x\mid\text{ gene }i\text{ is irreproducible}\right)$ $\displaystyle=$
$\displaystyle\begin{cases}1,\qquad&x<\pi_{1},\\\
1-\frac{(i_{x}-j_{\pi_{1}})^{2}}{(m-j_{\pi_{1}})^{2}},\quad&\pi_{1}\leq x\leq
1,\end{cases}$

where $i_{x}=\lfloor mx\rfloor$ and $j_{\pi_{1}}=\lfloor m\pi_{1}\rfloor$. The
limiting conditional survival function of $M_{i}/m$ under the null is

$\displaystyle S_{m,\pi_{1}}(x)\to S_{\pi_{1}}(x)=\begin{cases}1&x<\pi_{1}\\\
1-\frac{(x-\pi_{1})^{2}}{(1-\pi_{1})^{2}}&\pi_{1}\leq x\leq 1\\\
0&1<x\end{cases}.$

The empirical survival function can be estimated by
$\hat{S}_{m}(x)=\frac{1}{m}\sum_{i=1}^{m}I(M_{i}/m\geq x)$, $x\in(0,1)$. By
strong law of large numbers and Bayesian formula,

$\displaystyle\hat{S}_{m}(x)\rightarrow$ $\displaystyle P(M_{i}/m\geq x)$
$\displaystyle=$ $\displaystyle(1-\pi_{1})P(M_{i}/m\geq x\mid\text{ gene
}i\text{ is irreproducible})+\pi_{1}\times 0$ $\displaystyle=$
$\displaystyle(1-\pi_{1})S_{\pi_{1}}(x)\text{ for }x\in(\pi_{1},1).$

If we estimate $\pi_{1}$ by $i/m$, we can define the mean square error (MSE)
as follows.

$\mbox{MSE}(i/m)=(m-i)^{-1}\sum_{j=i}^{m}\left(\hat{S}_{m}(j/m)-(1-i/m)S_{i/m}(j/m)\right)^{2}.$

$\hat{k}$ is chosen to minimize the MSE in the range between $0$ and $\lfloor
0.9m\rfloor$.

$\displaystyle\hat{k}=\mathop{\arg\min}\limits_{i=0,1,\dots,\left\lfloor
0.9m\right\rfloor}\left\\{\mbox{MSE}(i/m)\right\\}.$

Thus $\hat{k}/m$ serves as a good estimation of $\pi_{1}$. To control the FDR
at level $\alpha$, the MaRR generates the rejection threshold as follows

$\displaystyle\text{Define }\hat{N}=\max_{\hat{k}<i\leq
n}\left\\{i:m\widehat{\mbox{FDR}}(i)=\frac{(i-\hat{k})^{2}}{Q(i)(m-\hat{k})}\leq\alpha\right\\},$

where $Q(i)=\sum_{j=1}^{m}I(M_{j}\leq i)$. Reject all features associated with
$M_{i}\leq\hat{N}$. Philtron et al. (2018) relax assumption (I1) to (R1):
$P(R_{1g}<R_{1h})>1/2$ and $P(R_{2g}<R_{2h})>1/2$, which is more plausible in
practice.

###  B.5 The radjust procedure

The radjust procedure (Bogomolov and Heller, 2018) works as follows,

  * •

Step 1. For a pre-specified FDR level $\alpha$, compute

$R=\mbox{max}\left[r:\sum_{i\in\mathcal{S}_{1}\cap\mathcal{S}_{2}}I\left\\{(p_{1i},p_{2i})\leq\left(\frac{r\alpha}{2|\mathcal{S}_{2}|},\frac{r\alpha}{2|\mathcal{S}_{1}|}\right)\right\\}=r\right],$

where $\mathcal{S}_{j}$ is the set of features pre-selected in study $j$ for
$j=1,2$. By default, it selects features with $p$-values less than or equal to
$\alpha/2$.

  * •

Step 2. Declare replicable features as those with indices in the set

$\mathcal{R}=\left\\{i:(p_{1i},p_{2i})\leq\left(\frac{R\alpha}{2|\mathcal{S}_{2}|},\frac{R\alpha}{2|S_{1}|}\right),i\in\mathcal{S}_{1}\cap\mathcal{S}_{2}\right\\}.$

In this paper, we implement the adaptive version of the radjust procedure
provided in Bogomolov and Heller (2018) for comparison, which first estimates
the fractions of true null hypotheses among the pre-selected features. The
fractions in the two studies are estimated as follows.

$\hat{\pi}_{0}^{(1)}=\frac{1+\sum_{i\in\mathcal{S}_{2,\alpha}}I(p_{1i}>\alpha)}{|\mathcal{S}_{2,\alpha}|(1-\alpha)},\quad\hat{\pi}_{0}^{(2)}=\frac{1+\sum_{i\in\mathcal{S}_{1,\alpha}}I(p_{2i}>\alpha)}{|\mathcal{S}_{1,\alpha}|(1-\alpha)},$
(B.1)

where $\mathcal{S}_{j,\alpha}=\mathcal{S}_{j}\cap\\{1\leq i\leq
m:p_{ji}\leq\alpha\\},\ j=1,2$. The adaptive procedure with a nominal FDR
level $\alpha$ works as follows.

  * •

Step 1. Compute $\hat{\pi}_{0}^{(1)}$ and $\hat{\pi}_{0}^{(2)}$ using (B.1).
Let

$R=\mbox{max}\left[r:\sum_{i\in\mathcal{S}_{1,\alpha}\cap\mathcal{S}_{2,\alpha}}I\left\\{(p_{1i},p_{2i})\leq\left(\frac{r\alpha}{2|\mathcal{S}_{2,\alpha}|\hat{\pi}_{0}^{(1)}},\frac{r\alpha}{2|\mathcal{S}_{1,\alpha}|\hat{\pi}_{0}^{(2)}}\right)\right\\}=r\right],$

  * •

Step 2. Reject features with indices in the set

$\mathcal{R}=\left\\{i:(p_{1i},p_{2i})\leq\left(\frac{R\alpha}{2|\mathcal{S}_{2,\alpha}|\hat{\pi}_{0}^{(1)}},\frac{R\alpha}{2|S_{1,\alpha}|\hat{\pi}_{0}^{(2)}}\right),i\in\mathcal{S}_{1,\alpha}\cap\mathcal{S}_{2,\alpha}\right\\}.$

###  B.6 The JUMP procedure

The JUMP procedure (Lyu et al., 2023) works on the maximum of $p$-values
across two studies. Define

$p_{i}^{\max}=\mbox{max}\\{p_{1i},p_{2i}\\},i=1,\dots,m.$

Let $s_{i}=(\theta_{1i},\theta_{2i}),\ i=1,\dots,m$ denote the joint hidden
states across two studies. Then $s_{i}\in\\{(0,0),(0,1),(1,0),(1,1)\\}$ with
$\mathbb{P}(s_{i}=(k,l))=\xi_{kl}$ for $k,l=0,1$ and $\sum_{k,l}\xi_{kl}=1$.
It can be shown that

$\displaystyle\mathbb{P}\left(p_{i}^{\max}\leq t\mid H_{i0}\text{ is
true}\right)$ $\displaystyle=$
$\displaystyle\frac{\xi_{00}\mathbb{P}(p_{i}^{\max}\leq t\mid
s_{i}=(0,0))}{\xi_{00}+\xi_{01}+\xi_{10}}+\frac{\xi_{01}\mathbb{P}(p_{i}^{\max}\leq
t\mid
s_{i}=(0,1))}{\xi_{00}+\xi_{01}+\xi_{10}}+\frac{\xi_{10}\mathbb{P}(p_{i}^{\max}\leq
t\mid s_{i}=(1,0))}{\xi_{00}+\xi_{01}+\xi_{10}}$ $\displaystyle\leq$
$\displaystyle\frac{\xi_{00}t^{2}+(\xi_{01}+\xi_{10})t}{\xi_{00}+\xi_{01}+\xi_{10}}\leq
t,$

which means that $p_{i}^{\max}$ follows a super-uniform distribution under the
replicability null. Denote

$G(t)=\frac{\xi_{00}t^{2}+(\xi_{01}+\xi_{10})t}{\xi_{00}+\xi_{01}+\xi_{10}}.$

For a given threshold $t\in(0,1)$, a conservative estimate of the FDR is
obtained by

$\mbox{FDR}^{*}(t)=\frac{m(\xi_{00}+\xi_{01}+\xi_{10})G(t)}{\sum_{i=1}^{m}I\\{p_{i}^{\max}\leq
t\\}\vee 1}.$

Following Storey (2002); Storey et al. (2004), the proportion of null
hypotheses in study $j$ can be estimated by

$\hat{\pi}_{0}^{(j)}(\lambda_{j})=\frac{\sum_{i=1}^{m}I\\{p_{ji}\geq\lambda_{j}\\}}{m(1-\lambda_{j})},\quad
j=1,2.$

Similarly, $\xi_{00}$ is estimated by

$\hat{\xi}_{00}(\lambda_{3})=\frac{\sum_{i=1}^{m}I\\{p_{1i}\geq\lambda_{3},p_{2i}\geq\lambda_{3}\\}}{m(1-\lambda_{3})^{2}},$

where $\lambda_{1},\lambda_{2}$ and $\lambda_{3}$ are tuning parameters that
can be selected by using the smoothing method provided in Storey and
Tibshirani (2003). Then we have

$\hat{\xi}_{01}=\hat{\pi}_{0}^{(1)}-\hat{\xi}_{00},\quad\hat{\xi}_{10}=\hat{\pi}_{0}^{(2)}-\hat{\xi}_{00}.$

With these estimates, we have a plug-in estimate of FDR,

$\widehat{\mbox{FDR}}^{*}(t)=\frac{m(\hat{\xi}_{00}t^{2}+\hat{\xi}_{01}t+\hat{\xi}_{10}t)}{\sum_{i=1}^{m}I\\{p_{i}^{\max}\leq
t\\}\vee 1}.$

The JUMP procedure works as follows.

  * •

Step 1. Let $p_{(1)}^{\max}\leq\dots\leq p_{(m)}^{\max}$ be the ordered
maximum of $p$-values and denote by $H_{(i)}$ the corresponding hypothesis;

  * •

Step 2. Find the largest $k$ such that the estimated FDR is controlled, i.e.,

$\hat{k}=\max\\{1\leq k\leq
m:\widehat{\mbox{FDR}}^{*}(p_{(k)}^{\max})\leq\alpha\\};$

  * •

Step 3. Reject $H_{(i)},$ $i=1,\dots,\hat{k}$.

###  B.7 Adjustment to our method

In the simulations, we estimate $(\xi_{00},\xi_{01},\xi_{10},\xi_{11})$ with
the conservative estimates proposed in Lyu et al. (2023) due to better finite
sample performance. Specifically, we estimate the proportion of null
hypotheses in study $j$, $\pi^{(j)}_{0}$, as

$\displaystyle\hat{\pi}_{0}^{(j)}(\lambda_{j})=\frac{\sum_{i=1}^{m}I\\{p_{ji}\geq\lambda_{j}\\}}{m(1-\lambda_{j})},j=1,2,$

and estimate $\xi_{00}$ as

$\displaystyle\hat{\xi}_{00}(\lambda_{3})=\frac{\sum_{i=1}^{m}I\\{p_{1i}\geq\lambda_{3},p_{2i}\geq\lambda_{3}\\}}{m(1-\lambda_{3})^{2}},$

where $\lambda_{1},\lambda_{2}$ and $\lambda_{3}$ are tuning parameters. We
use the smoothing method provided in Storey and Tibshirani (2003) to select
the tuning parameters. We have
$\hat{\xi}_{01}=\hat{\pi}_{0}^{(1)}-\hat{\xi}_{00},\hat{\xi}_{10}=\hat{\pi}_{0}^{(2)}-\hat{\xi}_{00}$
and $\hat{\xi}_{11}=1-\hat{\xi}_{00}-\hat{\xi}_{01}-\hat{\xi}_{10}.$