# A Unified View of False Discovery Rate Control: Reconciliation of Bayesian
and Frequentist Approaches

Xiaoquan Wen<EMAIL_ADDRESS>University of Michigan

###### Abstract

This paper explores the intrinsic connections between the Bayesian false
discovery rate (FDR) control procedures and their counterpart of frequentist
procedures. We attempt to offer a unified view of FDR control within and
beyond the setting of testing exchangeable hypotheses. Under the standard two-
groups model and the oracle condition, we show that the Bayesian and the
frequentist methods can achieve asymptotically equivalent FDR control at
arbitrary levels. Built on this result, we further illustrate that rigorous
post-fitting model diagnosis is necessary and effective to ensure robust FDR
controls for parametric Bayesian approaches. Additionally, we show that the
Bayesian FDR control approaches are coherent and naturally extended to the
setting beyond testing exchangeable hypotheses. Particularly, we illustrate
that $p$-values are no longer the natural statistical instruments for optimal
frequentist FDR control in testing non-exchangeable hypotheses. Finally, we
illustrate that simple numerical recipes motivated by our theoretical results
can be effective in examining some key model assumptions commonly assumed in
both Bayesian and frequentist procedures (e.g., zero assumption).

##  1 Introduction

Ever since the seminal work by Benjamini and Hochberg (Benjamini and Hochberg,
1995), false discovery rate (FDR) control has become increasingly popular in
the scientific practice of multiple hypothesis testing. Comparing to the
alternative statistical strategies, such as the control of family-wise error
rate (FWER), the FDR control is much more appealing in large-scale data
analyses due to its cost-effectiveness considerations. In the last decade,
both the theory and the computational algorithms for FDR control have matured
in the setting of testing multiple exchangeable hypotheses, thanks to the
important works of Benjamini and Hochberg (1995), Efron et al. (2001), Storey
(2002, 2003), Storey et al. (2004), Genovese and Wasserman (2004), Newton et
al. (2004), Müller et al. (2004), Ghosh et al. (2006), Sun and Cai (2007),
Whittemore (2007), Storey (2007), Efron (2008), Muralidharan (2010), Stephens
(2016), just name a few. In many fields, e.g., genetics, genomics, and
molecular biology, FDR control has emerged as the method of choice to quantify
uncertainty and safeguard potential type I errors in scientific discovery
processes.

The existing FDR control approaches can be classified into two distinct
categories. The first category is represented by the frequentist FDR control
approaches, e.g., the Benjamini-Hochberg (B-H) procedure and the Storey’s
$q$-value method, which are characterized by their use of $p$-values to
estimate and control FDR. The methods in the second category require explicit
computation of posterior probabilities, also commonly known as local fdr’s
(Efron et al., 2001), and are considered as Bayesian FDR control approaches.
The two categories of the methods are fundamentally different because the
underlying quantities that measure the expected false discovery proportion are
conceptually distinct (the details are explained in Section 2). Both
categories of approaches have their own strength and weakness. For example,
the frequentist FDR control methods, which requires no explicit specification
of the alternative models, are known to be robust to the misspecification of
the alternative models; whereas the Bayesian methods, especially the
parametric approaches, have great flexibility to incorporate ancillary
data/prior information that can significantly improve the power of hypothesis
testing.

There have been many efforts to bridge the gap between the two types of
approaches in the literature: Storey (Storey, 2002) shows that the Bayesian
interpretation of the $q$-value method and Efron (Efron, 2008) demonstrates
the connections between Bayesian FDR/local fdr and frequentist FDR under some
specific power functions. In this article, we aim to explore a more in-depth
connection between the two types of approaches. In particular, we are
interested in investigating if the Bayesian and frequentist FDR control
procedures can yield concordant results for analyzing the same data set. If
the answer is yes, we would like to identify the corresponding necessary and
sufficient conditions. Answers to these questions have important implications
for the best practice of FDR controls for both frequentist and Bayesian
approaches. Beyond the standard setting characterized by the two-groups model
(Efron, 2008), we further explore the extensibility an optimality of the two
types of FDR control methods in some more complex, but increasingly more
realistic, scenarios.

The paper is structured as the follows. We first provide necessary background
on frequentist and Bayesian FDR controls, then proceed to present our main
theoretical results. Subsequently, we provide numerical illustrations of our
theoretical results. Finally, we end the paper by the summarizing the main
conclusions and a brief discussion on future directions.

##  2 Background

###  2.1 Model and notation

We follow Genovese and Wasserman (2004), Efron (2008) to adopt a mixture model
description of testing $m$ exchangeable hypotheses as our starting point. Let
the latent binary indicator $\gamma_{i}=1$ denote that the observed data of of
the $i$-th hypothesis are generated from the alternative scenario, and 0
otherwise. We use $z_{i}$ to represent the observed data/summary statistic
from the $i$-th test, which can be computed independent of data from all other
tests. Note that $z_{i}$ can be a $p$-value in a frequentist approach or a
Bayes factor (i.e., a marginal likelihood ratio statistic) in a Bayesian
approach; it can also be a vector (Stephens, 2016) instead of a single number.
We denote the collection of the observed data and the corresponding latent
indicators by $\mbox{\boldmath$z$}=(z_{1},\dots,z_{m})$ and
$\mbox{\boldmath$\gamma$}=(\gamma_{1},\dots,\gamma_{m})$, respectively.

The assumed data generation mechanism can be described by the following
probabilistic model,

$\displaystyle\gamma_{i},\dots,\gamma_{m}\mid\pi_{0}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}{\rm Bernoulli}\,(1-\pi_{0}),$ (1) $\displaystyle
z_{i}\mid\gamma_{i}=0\,\sim f_{0},$ $\displaystyle z_{i}\mid\gamma_{i}=1\,\sim
f_{1}.$

In particular, we assume the independence of $\gamma_{i}$’s conditional on the
hyper-parameter $\pi_{0}$. The model also implies that $z_{i}$’s are
exchangeable, and we will relax this particular assumption later in the paper.
This particular mixture model is known as the _two-groups model_ (Efron,
2008). It is also a special case of conditional independent hierarchical model
(CIHM) discussed in Kass and Steffey (1989). Unless otherwise specified, we
assume that the theoretical null, $f_{0}$, is known and correctly specified
while both $f_{1}$ and $\pi_{0}$ are unknown. Finally, we use $f_{c}$ to
denote the mixture distribution $\pi_{0}f_{0}+(1-\pi_{0})f_{1}$.

###  2.2 FDR estimation and control

Let $\delta_{i}(\cdot)$ denote a binary decision function with respect to the
$i$-th hypothesis, i.e., $\delta_{i}=1$ indicates the rejection of the $i$th
hypothesis, and 0 otherwise. Under the two-groups model and given a decision
function, we define the false discovery proportion (FDP), a random variable,
by

$\displaystyle{\rm FDP}$ $\displaystyle=\frac{\mbox{number of falsely rejected
null hypothese}}{\max(1,\mbox{number of total rejections})}$ (2)
$\displaystyle=\frac{\sum_{i=1}^{m}\delta_{i}\,(1-\gamma_{i})}{\left(\sum_{i=1}^{m}\delta_{i}\right)\,\vee\,1}.$

Both Bayesian and frerquentist FDR control procedures aim to bound some type
of expectation of FDP.

####  2.2.1 Bayesian FDR and the optimal control procedure

The Bayesian FDR is defined as the conditional expectation of FDP given the
observed data. Because of the explicit conditioning, $\delta_{i}$’s are
considered deterministic. Müller et al. (2006) shows that the optimal decision
rule to control Bayesian FDR is given by

$\delta_{B,i}(t_{b})={\bf 1}_{\\{u_{i}\,\leq\,t_{b}\\}},$ (3)

where $u_{i}$ denotes the posterior probability
$\Pr(\gamma_{i}=0\mid\mbox{\boldmath$z$})$ and $t_{b}$ represents a pre-
defined threshold on $u_{i}$’s. The optimal Bayesian FDR control procedure
rejects the null hypotheses for which the individual posterior probabilities
of false discovery, i.e., $u_{i}$’s, are small. It straightforwardly follows
from the definition that

${\rm BFDR}(t_{b})={\rm E}({\rm
FDP}\mid\mbox{\boldmath$z$})=\frac{\sum_{i=1}^{m}u_{i}\,\delta_{B,i}(t_{b})}{\left[\sum_{i=1}^{m}\delta_{B,i}(t_{b})\right]\,\vee\,1}.$
(4)

In the literature, the posterior probability of false discovery is also
commonly referred to as local fdr, a terminology coined by Efron.

####  2.2.2 Frequentist FDR control

The frequentist FDR is defined as the unconditional expectation of FDP, i.e.,

${\rm FDR}={\rm E}({\rm FDP})={\rm E}({\rm BFDR}),$ (5)

where, in comparison to Bayesian FDR, an additional expectation is taken with
respect to the hypothetically re-sampled data, $z$. Widely used Benjamini-
Hochberg (B-H) and Storey’s $q$-value algorithms are designed to control
frequentist FDR. The decision rules embedded in those algorithms have the
common form,

$\delta_{F,i}(t_{f})={\bf 1}_{\\{p_{i}\,\leq\,t_{f}\\}},$ (6)

where $p_{i}$ and $t_{f}$ represent the $p$-value of test $i$ and a pre-
defined threshold, respectively. While the exact evaluation of frequentist FDR
(5) is difficult for a given $t_{f}$, Genovese and Wasserman (2004) shows an
accurate asymptotic approximation,

${\rm
FDR}(t_{f})=\frac{m\,\pi_{0}\,t_{f}}{\left[\sum_{i=1}^{m}\delta_{F,i}(t_{f})\right]\,\vee\,1}+o\left(\frac{1}{m}\right).$
(7)

The approximation,
$\frac{m\,\pi_{0}\,t_{f}}{\left[\sum_{i=1}^{m}\delta_{F,i}(t_{f})\right]\,\vee\,1}$,
is also referred to as an estimate of Fdr by Efron (2012). In Sun and Cai
(2007) the same quantity is refers to as marginal FDR (mFDR).

The B-H and the $q$-value procedures are two most commonly applied frequentist
FDR control methods. The difference between the two lies in their treatment of
unknown parameter $\pi_{0}$. The $q$-value procedure plugs in a non-parametric
estimator of $\pi_{0}$ that is derived based on sample quantile information.
In comparison, the B-H procedure simply sets $\pi_{0}=1$ in (7). As a result,
its estimate essentially becomes an upper-bound of desired FDR.

###  2.3 Rejection path and control of FDR

In this paper, we propose a new statistical device to perform comparisons
between FDR control procedures. We formally define a rejection path as
follows.

######  DEFINITION 1.

Given a decision rule of a FDR control procedure, $\delta_{i}={\bf
1}_{\\{s_{i}\leq t\\}}$, the rejection path for a data set is the sequence of
estimated (Bayesian or frequentist) false discovery rates sequentially
evaluated at the threshold values determined by the order statistics,
$s_{(1)},\dots,s_{(m)}$, i.e.,

${\rm FDR}(t=s_{(1)}),{\rm FDR}(t=s_{(2)}),\dots,{\rm FDR}(t=s_{(m)}).$

The rejection path contains the complete information required to control FDR
at any pre-defined level $\alpha$. Specifically, we find

$l=\arg\max_{k}\\{k:{\rm FDR}(s_{(k)})\leq\alpha\\},$

and reject the tests corresponding to $s_{(1)},...,s_{(l)}$.

Different FDR control procedures typically yield distinct rejection paths when
analyzing a given data set. If the rejection paths from two procedures are
always identical, the following two conditions are necessarily satisfied.
First, the two procedures rank individual tests in the same order. In
addition, the two procedures yield the same FDR estimates at all threshold
values in the rejection path. The ranking of individual tests is typically
determined by the choice of test statistics and can be examined by the method
like receiver operating characteristic (ROC) curves in simulation settings.
The evaluation of FDR is the one of the main interests of this paper.

Direct comparison of two FDR control procedures yields quantitative assessment
of relative conservativeness of the different FDR estimation approaches.
Figure 1 illustrates the pair-wise comparison of the rejection paths from the
B-H procedure and the $q$-value procedure for the Hedenfalk data distributed
in the R package q-value.

Figure 1: pairwise comparison of rejection paths of the $q$-value method and
the B-H procedure. We use the Hedenfalk data distributed with the R package
$q$-value, which contain 3,170 $p$-values. The comparison of the rejection
paths indicates that the B-H procedure always overestimate the FDR than the
$q$-value procedure by a linear factor.

##  3 Results

###  3.1 Frequentist property of Bayesian FDR control

Our first new result concerns the quantitative connection between the Bayesian
and frequentist FDR control procedures in a setting where $m$ is large. We
consider an ideal scenario where both $\pi_{0}$ and $f_{1}$ are correctly
specified, which corresponds to an oracle condition (Sun and Cai, 2007). Under
such setting, the posterior false discovery probability for each test,
$u_{i}^{*}=\Pr(\gamma_{i}=0\mid\pi_{0},z_{i})$, can be analytically computed
by

$\displaystyle u_{i}^{*}$
$\displaystyle=\frac{\pi_{0}f_{0}(z_{i})}{\pi_{0}f_{0}(z_{i})+(1-\pi_{0})f_{1}(z_{i})}$
(8) $\displaystyle=\frac{\pi_{0}}{\pi_{0}+(1-\pi_{0}){\rm BF}^{*}(z_{i})},$

where ${\rm BF}^{*}(z_{i}):=f_{1}(z_{i})/f_{0}(z_{i})$ denotes the marginal
likelihood ratio statistic/Bayes factor.

We show that under specific conditions, Bayesian and frequentist FDR control
procedures are asymptotically equivalent. The main result is summarized in
Theorem 1.

######  THEOREM 1.

Assume that

  1. 1.

$\pi_{0}$, $f_{0}$ and $f_{1}$ are correctly specified;

  2. 2.

There exists a monotone mapping from the marginal likelihood ratio statistics
to the $p$-values.

It follows that

${\rm BFDR}\left(t_{b,i}=u^{*}_{(i)}\right)\xrightarrow{a.s.}{\rm
FDR}\left(t_{f,i}=p_{bf,(i)}\right)\mbox{ for all }i.$

###### Proof.

See Appendix A. ∎

Remark 1. The existence of monotone mapping from the Bayes factors to the
$p$-values implies that the monotonic relationship between the two statistical
instruments holds for any observed data sets. The most straightforward way to
guarantee this condition is to derive $p$-values directly from the
corresponding likelihood ratio statistics. In this scenario, Bayesian and
frequentist analyses make exactly the same modeling assumptions and extract
exactly the same information from data, and the difference in results should
solely reflect their procedural difference in FDR controls. Nevertheless, it
should be noted that the monotonic relationship between Bayes factors and
$p$-values can hold in much more general settings. For example, Bayesian
hypothesis testing and inference based on $z$-statistics (or
$z^{2}$-statistics) has been widely studied and applied in a wide range of
applications (e.g., genetic association studies). It is commonly assumed that
$z_{i}\sim{\rm N}(0,1)$ under the null and $z_{i}\sim{\rm N}(0,1+k)$ (for some
$k>0$) under the alternative. (Note that the alternative model is sometimes
interpreted as a random effect assumption.) Under this setting, both Cox and
Hinkley (1979) and Wakefield (2009) have shown that

${\rm
BF}(z_{i})=\frac{1}{\sqrt{1+k}}\exp\left(\frac{1}{2}\,\frac{k}{1+k}\,z_{i}^{2}\right),$
(9)

which is a monotonic transformation of $z_{i}^{2}$. As a consequence, the
$p$-value of the Bayes factor (9) is identical to the two-sided $p$-value of
the corresponding $z$-statistic. More generally, the above relationship holds
for modeling of $z$-statistics under alternatives using mixtures of normal
and/or uniform distributions under the uni-modal alternative (UA) assumption
that is recently proposed by Stephens (2016).

Remark 2. Both Müller et al. (2006) and Sun and Cai (2007) have shown that the
oracle Bayesian procedure is optimal to minimize false non-discovery rate
(FNR, a quantity measuring type II errors) while controlling FDR. Our result
essentially extends the optimality results to the frequentist FDR control
procedures due to the asymptotic equivalence.

Remark 3. The key assumption on the correctness of model specification is
critical to ensure the convergence result. Although it is unrealistic to
expect this particular assumption to hold _exactly_ in practice, the high-
dimensional setting of practical applications typically yields a large amount
of informative data and offers an opportunity for accurate model estimation.
Additionally, it underscores the necessity and importance of careful model
specification and thorough diagnosis in the practice of Bayesian FDR control.
This topic will be extensively discussed in the subsequent sections.

Theorem 1 goes beyond the previous results by Efron (2012) and Storey (2003),
both of which focused on the Bayesian representation/interpretation of the
frequentist FDR estimation, i.e., the functional form of equation (7). It
highlights the point-wise convergence of rejection paths of the Bayesian and
frequentist procedures under the oracle condition. Thus, in a high-dimensional
setting and under correct model specifications, the Bayesian and frequentist
FDR control procedures are expected to be highly concordant for all pre-
defined FDR control levels. Figure 2 provides a numerical illustration of
Theorem 1 using a set of simulated data, where good numerical concordance
between the Bayesian and frequentist procedures can be observed when $m\sim
10^{3}$.

|  |
---|---|---
Figure 2: pairwise comparison of rejection paths from the Bayesian and
frequentist FDR control procedures under the oracle condtion. In all three
panels, $z$-scores are drawn from the mixutre distribution,
$\pi_{0}N(0,1)+(1-\pi_{0})N(0,1+k)$, where we set $\pi_{0}=0.5$ and $k=10$.
Each panel represents a simulation with $m$ ($=200,2000$ and 20,000,
respectively) independent tests.

###  3.2 Approximate inference of Bayesian FDR

Standard Bayesian inference requires parametric specification of $f_{1}$ by
$f_{\theta}$, where $\theta$ denotes the necessary hyper-parameters. Ideally,
$f_{\theta}$ should be flexible enough such that it can accurately approximate
the true $f_{1}$. The examples of such kind are illustrated by Stephens
(2016). To enable Bayesian computation, the priors for the set of all the
hyper-parameters $\lambda:=(\pi_{0},\theta)$ also need to be specified.

One of the factors that hinder the usage of the standard Bayesian inference
for FDR control in practical settings is its expensive computational cost due
to lacking of general analytic solutions for computing local fdr’s in a
general CIHM. Markov chain Monte Carlo (MCMC) algorithms are theoretically
possible, however, the lack of scalability and/or effective convergence
diagnosis often makes applying Bayesian procedure in high-dimensional settings
less appealing. Here, we advocate the use of parametric empirical Bayes (PEB)
procedure to fit CIHM and argue that PEB generally provides adequate
approximations to exact Bayesian inference in high-dimensional hypothesis
testing settings. Comparing to the full Bayesian inference, the key
distinction in the PEB procedure is to replace the desired local fdr,
$\Pr(\gamma_{i}=0\mid\mbox{\boldmath$z$})$, with an empirical Bayes estimate,
namely $\Pr(\gamma_{i}=0\mid\hat{\lambda},\,z_{i})$, where
$\hat{}\mbox{\boldmath$\lambda$}$ denotes the MLE of $\lambda$ and is a
function of all observed data $z$. This procedure is justified by the result
of Kass and Steffey (1989), which shows that the PEB approach can be regarded
as a _computational approximation_ of the full Bayesian computation in the
setting of CIHM (for reasonably specified, yet arbitrary, priors on
$\lambda$). More specifically, Result 1 of Kass and Steffey (1989) implies
that, with a pre-specified alternative model $f_{\theta}$,

$\Pr(\gamma_{i}=0\mid\mbox{\boldmath$z$})=\Pr(\gamma_{i}=0\mid
z_{i},\hat{\lambda})\cdot\left[1+o(\frac{1}{m})\right],$ (10)

which is obtained by applying Laplace’s method and re-normalizing the
resulting approximation into a valid probability distribution.

The PEB approximation is particularly attractive for at least two reasons.
First, its result has a _relative_ error of $o(\frac{1}{m})$ which is
appealing for approximating small to modest posterior probabilities. In
comparison, traditional Monte Carlo based approach only ensures the control of
_additive_ error at the same scale. This point is critically important because
the accurate evaluation of small to modest posterior probabilities is the key
to ensuring accurate Bayesian FDR control. Second, the EB procedure, which can
be carried out by some well-established optimization algorithms (e.g., EM), is
much more computationally efficient than the MCMC fitting of CIHM and is
particularly suitable in a big data setting. It should be clear that the
overall effectiveness of EB-based approximate inference procedure is dictated
by the informativeness of the data with respect to the key model parameters.
In this sense, the high-dimensionality of the simultaneous tests creates a
nearly ideal asymptotic setting for accurate parameter estimation.

The PEB also produces an estimate of the mixture density $f_{c}$, i.e.,

$\hat{f}_{c}=\hat{\pi}_{0}f_{0}+(1-\hat{\pi}_{0})f_{\hat{\theta}}.$ (11)

Interestingly, by the same argument of Kass and Steffey (1989), it can be
shown that $\hat{f}_{c}$ is an approximation of the density of the posterior
predictive distribution, $f(z_{\rm new}\mid\mbox{\boldmath$z$})$, i.e.,

$f(z_{\rm new}\mid\mbox{\boldmath$z$})=\int f(z_{\rm
new}\mid\lambda)f(\lambda\mid\mbox{\boldmath$z$})d\,\lambda=\hat{f}_{c}(z_{\rm
new})\left[1+o(\frac{1}{m})\right].$ (12)

The posterior predictive distribution is instrumental for Bayesian model
diagnosis, which is critical to ensure robust Bayesian FDR control. The usage
of $\hat{f}_{c}$ for Bayesian model checking will be illustrated in the next
section.

###  3.3 Robustness and model diagnosis

Model misspecification and inadequate fitting can lead to either conservative
or anti-conservative behaviors in FDR control. While the former (e.g.,
overestimating $\pi_{0}$) is typically accepted, the latter (e.g.,
underestimating $\pi_{0}$) inflates FDR and is considered dangerous. It is
therefore critical to ensure the robustness of the FDR control by either
theoretical arguments and/or practical model diagnosis procedures. In
particular, we expect an effective diagnosis approach could not only flag
improper modeling assumptions and/or inadequate fitting but also provide a
_directional_ indication of anti-conservative FDR control results.

The frequentist FDR control procedures mainly rely on the $p$-values derived
from the null hypothesis and seemingly avoid the detailed specifications of
alternative models. Nevertheless, a _precise_ frequentist FDR estimation
requires an accurate estimate of $\pi_{0}$, which inevitably needs the
knowledge of the alternative model. A common strategy to circumvent this
difficulty is to use an upper bound estimate of $\pi_{0}$ and acknowledge that
the conservativeness of the resulting procedure. For example, the B-H
procedure simply plugs in the most natural upper bound by setting $\pi_{0}=1$
in Equation (5). The $q$-value procedure takes advantages of the high-
dimensional setting of multiple hypothesis testing and finds a tighter upper
bound estimate for $\pi_{0}$ based on the sample quantile information. More
specifically, for any $\eta\in(0,1)$ and noting that $p$-values are uniformly
distributed under the null model, it follows from the Kolmogorov’s strong law
of large numbers that

$\frac{1}{m}\sum_{i=1}^{m}{\bf
1}\\{p_{i}\geq\eta\\}\xrightarrow{a.s.}\pi_{0}\,(1-\eta)+(1-\pi_{0})\,[1-F_{p^{*}}(\eta)],$
(13)

where $F_{p^{*}}(\cdot)$ denotes the CDF of the $p$-values under the (unknown)
alternative model. By ignoring the second term on the right-hand side, an
upper bound estimate of $\pi_{0}$ is obtained by

$\tilde{\pi}_{0}=\min\left(1,\frac{\sum_{i=1}^{m}{\bf
1}\\{p_{i}\geq\eta\\}}{m\,(1-\eta)}\right),$ (14)

and it satisfies the upper bound condition in the sense,

$\tilde{\pi}_{0}\geq\pi_{0}\mbox{ almost surely},\mbox{ as }m\to\infty.$ (15)

Theoretically, the upper bound property of $\tilde{\pi}_{0}$ holds for
arbitrary $\eta$ values, although setting $\eta$ close to 0.50 reduces the
variance of the estimate. In practice, some use $\tilde{\pi}_{0}$ extrapolated
at $\eta\to 1$ and argue it provides a tighter bound assuming that most of
large $p$-values are generated from the null hypothesis. However, this
assumption can be false if the alternative scenario assumed by the unimodal
alternative (UA) assumption is indeed true, where the data generated from the
alternative models can be very much likely to produce large $p$-values (see
discussions in Stephens (2016)). Nevertheless, in all cases, $\tilde{\pi}_{0}$
should provide a reliable upper bound estimate and ensures the robust
performance of FDR control. Finally, we note that estimator (14) is not
restricted to the usage of $p$-values and can be generalized to almost any
test statistics with available quantile information from the null
distribution. For example, if the $p$-values are derived from the two-sided
tests of $z$-statistics, the following estimator based on the $\eta$-quantile
of $z^{2}$, denoted by $\chi^{2}_{\eta}$, provides the identical estimation as
(14),

$\tilde{\pi}_{0}=\min\left(1,\frac{\sum_{i=1}^{m}{\bf
1}\\{z^{2}_{i}\leq\chi^{2}_{\eta}\\}}{m\,\eta}\right).$

The locfdr procedure (Efron et al., 2001, Efron, 2008) utilizes a non-
parametric empirical Bayes (NPEB) inference framework, where the density of
the mixture distribution, $f_{c}$, is estimated non-paramterically. To
estimate $\pi_{0}$, it makes the “zero assumption” (ZA), i.e., $f_{1}(z)\ll
f_{0}(z)$ as $z\to 0$ in case of simple $z$-statistics. Similar to the case of
estimator (14), we view the ZA as a robust assumption for FDR control: the
violation of the assumption results in over-estimated $\pi_{0}$. In general,
locfdr ensures an accurate fit of $f_{c}$, if the ZA holds, the resulting
estimate is accurate and the subsequent FDR control is optimal or near-
optimal; otherwise, it overestimates $\pi_{0}$ and results in conservative
control of FDR.

In comparison, the standard Bayesian inference with the parametric
specification of $f_{\theta}$ is more susceptible to model misspecification
and/or inadequate fitting. Here we emphasize the critical importance of
conducting post-fitting model diagnosis utilizing the posterior predictive
distribution. In the case that Bayesian inference is approximated by the PEB
procedure, the standard posterior predictive checking in Bayesian inference
converges to the common practice of examining the goodness-of-fit of
$\hat{f}_{c}$, based on the argument provided by equation (12). In particular,
we find that it is effective to examine the alignment between the sample
quantiles from the observed data and the corresponding theoretical quantiles
computed (or estimated) from $\hat{f}_{c}$, which is an idea naturally
extended from Cook et al. (2006). In some of the examples that we show in
Section 4, the patterns of significant misalignment of the two quantiles also
provide an important indication if a particular inadequate fit leads to a
conservative or anti-conservative FDR control.

To provide a theoretical guidance to Bayesian model diagnosis, we extend
Theorem 1 to a practical setting where the inference of $\pi_{0}$ and $f_{1}$
is required. Specifically, we consider the estimates $\hat{\pi}_{0}$ and
$\hat{f}_{c}$ (or equivalently, $f_{\hat{\theta}}$) are obtained through an
inference procedure (e.g., PEB), and the local fdr for test $i$ is evaluated
by
$\hat{u}_{i}:=\hat{u}(z_{i};\hat{\pi}_{0})=\frac{\hat{\pi}_{0}f_{0}(z_{i})}{\hat{f}_{c}(z_{i})}$.
Furthermore, we denote the induced order statistics from
$\\{\hat{u}_{1},...,\hat{u}_{m}\\}$ by
$(\hat{u}_{(1)},\hat{u}_{(2)},...,\hat{u}_{(m)})$. Our main result is
summarized in Theorem 2.

######  THEOREM 2.

Assume that

  1. 1.

there exists a monotone mapping from the marginal likelihood ratio statistics,
$f_{0}/\hat{f}_{1}$, to the corresponding $p$-values;

  2. 2.

${\rm E}\,||\hat{f}_{c}-f_{c}||^{2}\to 0$.

It follows that

${\rm BFDR}\left(t_{b,i}=\hat{u}_{(i)};\hat{\pi}_{0}\right)\xrightarrow{p}{\rm
FDR}\left(t_{f,i}=p_{(i)};\hat{\pi}_{0}\right)\mbox{ for all }i.$

Furthermore, if $\hat{\pi}_{0}\xrightarrow{p}\pi_{0}$,

${\rm BFDR}\left(t_{b,i}=\hat{u}_{(i)};\hat{\pi}_{0}\right)\xrightarrow{p}{\rm
BFDR}\left(t_{b,i}=u^{*}_{(i)}\right)$

###### Proof.

See Appendix B. ∎

Remark. Because the identifiability issue associated with mixture model
inference can be extremely complicated, the overall measure of adequate
fitting represented by the second assumption of Theorem 2 does not generally
guarantee that $\hat{\pi}_{0}$ is accurat in theory. Thus, our convergence
result between ${\rm BFDR}$ and ${\rm FDR}$ is conditional on the
$\hat{\pi}_{0}$ estimate, and the accuracy of the $\pi_{0}$ estimate needs to
be further ensured to establish that the FDR control results are approaching
the theoretical optimal.

Theorem 2 implies that if a Bayesian model is adequately fit by the parametric
empirical Bayes approach, we expect a strong agreement between the Bayesian
and the corresponding frequentist FDR control procedures, provided that the
same estimate of $\hat{\pi}_{0}$ is applied to both procedures. Conversely,
notable discordance between the Bayesian and frequentist approaches indicates
poor modeling assumption and/or inadequate fitting of the assumed parametric
model.

Note that Theorem 2 is not directly applicable to the NPEB approach
implemented in the locfdr method, despite that the overall adequate fitting of
the mixture distribution, $f_{c}$, is usually achieved. This is because the
Bayes factors/likelihood ratio statistics are implicit (i.e., no explicit form
of $f_{\theta}$ is estimated), and the theoretical validation of the first
assumption is not plausible. Nevertheless, in practice, the rank correlation
between local fdr’s from the NPEB approach and the corresponding p-values can
be numerically examined. As the rank correlation $\to 1$, it is expected,
based on Theorem 2, that the two approaches yield concordant FDR control
results. This point will be fully illustrated in Section 4.

One of the important practical implications of Theorem 2 is that the goodness-
of-fit of the mixture distribution (i.e., the second assumption) is a
necessary condition to ensure the frequentist property of the Bayesian FDR
control. This result underscores the importance of model checking as a
critical part of an overall Bayesian FDR control procedure. Beyond the scope
of Theorem 2, exploring the detailed patterns of misfitting by $\hat{f}_{c}$
may offer valuable insights in determining whether the inadequate fit leads to
conservative or anti-conservative FDR control.

###  3.4 Multiple hypothesis testing with non-exchangeable data

In many practical settings of multiple hypothesis testing, there exist
ancillary data that make the observed data non-exchangeable. For example, in
the analysis of microarray/RNA-seq experiments for identifying differentially
expressed genes, the genes are naturally grouped into gene sets by their
biological relevance/pathway information. It is often suspected, a priori,
that genes in certain sets are more (or less) likely to be differentially
expressed under the specific experimental condition.

Without loss of generality, we assume that testing data can be classified into
$K$ mutually exclusive groups based on some known categorical annotation
$\mbox{\boldmath$d$}=(d_{1},\dots,d_{m})$, where $d_{i}$ is the corresponding
group label for the data involved in the $i$-th test. Under this setting, the
exchangeability assumption in the original two-group model is no longer valid
due to the availability of the covariates, $d$. Nevertheless, it is
straightforward to extend the CIHM from the two-group model into the following
form to accommodate the non-exchangeable data:

$\displaystyle\gamma_{i}\mid d_{i}\sim{\rm Bernoulli}(1-\pi_{d_{i},0}),$ (16)
$\displaystyle z_{i}\mid\gamma_{i}=0,\,d_{i}\sim f_{d_{i},0},$ $\displaystyle
z_{i}\mid\gamma_{i}=1,\,d_{i}\sim f_{d_{i},1}.$

The Bayesian FDR control is coherent in dealing with the non-exchangeable data
in multiple hypothesis testing. Specifically, the same decision rule (3) is
applied and the only notable difference is that the false discovery
probability for the $i$-th test,
$\Pr(\gamma_{i}=1\mid\mbox{\boldmath$z$},d_{i})$ is now additionally
conditional on the covariate data, $d_{i}$. Under the oracle condition, i.e.,
as the quantities $\\{(\pi_{k,0},f_{k,0},f_{k,1}):k=1,...,K\\}$ are known, the
desired local fdr is computed by

$\Pr(\gamma_{i}=1\mid\mbox{\boldmath$z$},d_{i})=\frac{\pi_{d_{i},0}\,f_{d_{i},0}(z_{i})}{\pi_{d_{i},0}\,f_{d_{i},0}(z_{i})+(1-\pi_{d_{i},0})\,f_{d_{i},1}(z_{i})},$

which is monotonic to the weighted likelihood ratio (wlr) statistic,

$wlr^{*}(z_{i};d_{i})=\frac{1-\pi_{d_{i},0}}{\pi_{d_{i},0}}\,\frac{f_{d_{i},1}(z_{i})}{f_{d_{i},0}(z_{i})}.$

Importantly, under the oracle condition for multiple testing with non-
exhangeable data, Sun and Cai (2007) shows that the Bayesian FDR control is
optimal in the sense that it minimizes FNR while controlling the desired FDR
level. A natural question to ask is: is there asymptotically equivalent
frequentist procedure corresponding to the Bayesian procedure?

Let $F_{k,wlr^{*}}$ denote the CDF of the wlr statistic,
$\frac{1-\pi_{k,0}}{\pi_{k,0}}\,\frac{f_{k,1}}{f_{k,0}}$, for group $k$ under
the null. We find the following procedure seemingly provides a frequentist FDR
control equivalent.

######  ALGORITHM 1.

(Frequentist FDR control of multiple hypothesis testing with non-exchangeable
data)

  1. 1.

Obtain the CDF of the weighted likelihood ratio statistics under the null
hypothesis, $F_{k,wlr^{*}}$, for each group $k=1,...,K$.

  2. 2.

At a given threshold of the wlr statistics, namely t, reject the set of tests
$\\{i:wlr^{*}_{i}\geq t\\}$ across all groups.

  3. 3.

Evaluate the frequentist FDR at threshold $t$ by

${\rm
FDR}(t)=\frac{\sum_{i=1}^{m}\pi_{d_{i},0}[1-F_{d_{i},wlr^{*}}(t)]}{\left[\sum_{i=1}^{m}{\bf
1}_{\\{wlr^{*}_{i}\geq t\\}}\right]\vee 1}.$ (17)

In practice, to apply Algorithm 1, we rank the wlr statistics from all tests
(across groups) in descending order and evaluate ${\rm FDR}(t=wlr^{*}_{(i)})$
sequentially until the pre-defined FDR control level is achieved. We show, in
Proposition 1, that Algorithm 1 controls frequentist FDR and is equivalent to
its Bayesian counterpart as $m\to\infty$.

######  PROPOSITION 1.

Under the oracle condition, the frequentist procedure described in Algorithm 1
is asymptotically equivalent to the Bayesian FDR control procedure for testing
multiple hypotheses with non-exchangeable group structures, i.e.,

${\rm BFDR}\left(t_{b,i}=u^{*}_{(i)}\right)\xrightarrow{a.s.}{\rm
FDR}\left(t_{f,i}=wlr_{(i)}^{*}\right)\mbox{ for all }i.$

###### Proof.

See Appendix C. ∎

Because of the asymptotic equivalence, we conclude that Algorithm 1 also
shares the optimality of the Bayesian procedure under the setting. Although
the quantity $1-F_{d_{i},wlr^{*}}$ may be interpreted as a group-specific
$p$-value, it is not directly involved in the decision rule as in the case for
non-exchangeable data. Furthermore, note that the ranking of the wlr
statistics is generally different from the ranking of the group-specific
$p$-values across groups. Many authors (Genovese et al., 2006, Ignatiadis et
al., 2016) also have explicitly sought the decision rules of the following
form for testing with non-exchangeable data,

$\delta_{i}(t)={\bf 1}_{\\{\frac{p_{i}}{w_{d_{i}}}\leq t\\}},$ (18)

where $p_{i}$ represents the $p$-value derived from $z$ for the $i$-th test
and $w_{d_{i}}$ represents a group-specific weight. It should be noted that
weighting $p$-values generally results in a very different rejection path than
weighting the likelihood ratio statistics because of the irreconcilable
difference in ranking individual tests and is therefore unlikely, if not
impossible, to achieve the optimal performance.

In a practical setting when $\\{(\pi_{k,0},f_{k,0},f_{k,1}):k=1,...,K\\}$ are
unknown, the same approximate Bayesian inference (i.e., PEB) and model
diagnosis strategies apply for the Bayesian FDR control procedure. On the
other hand, extending the optimal oracle frequentist FDR control seems non-
trivial: unlike in the two-group model setting, it is required to make
explicit assumption and inference on the distributions from the alternative
scenarios for each group $k$. As a consequence, the operational procedure of
an ideal frequentist FDR control method in this setting may inevitably
resemble its Bayesian counterpart.

Finally, we note that the parametric Bayesian approach can be consistently
extended to account for continuous covariate data $d$. As a starting point, a
CIHM can be specified as

$\displaystyle\gamma_{i}\mid d_{i}\sim{\rm
Bernoulli}\,\left(\frac{\exp[\alpha_{0}+\alpha_{1}d_{i}]}{1+\exp[\alpha_{0}+\alpha_{1}d_{i}]}\right),$
(19) $\displaystyle z_{i}\mid\gamma_{i}=0\,\sim f_{0},$ $\displaystyle
z_{i}\mid\gamma_{i}=1\,\sim f_{1}.$

Note that the logistic prior model for $\gamma_{i}\mid d_{i}$ is applicable
for the categorical covariate data: it is simply an alternative
parametrization of the $\gamma_{i}$ prior in (16). Furthermore, it is also
possible to specify parametric distributions of $f_{0}$ and $f_{1}$ that are
dependent on the continuous covariate data. However, these modeling choices
are largely context-dependent and outside the scope of this paper.

##  4 Numerical illustrations

###  4.1 locfdr vs. $q$-value

The locfdr and the $q$-value procedures are two commonly applied procedures
for Bayesian and frequentist FDR controls, respectively. The implementation of
the locfdr method represents a frequentist inference alternative to estimate
the Bayesian FDR, as it circumvents prior specification of $\pi_{0}$ and the
explicit computation of parametric likelihood for alternatives. For the very
reason, it is difficult to evaluate the rejection path of the local fdr method
via Theorem 2. Here we conduct numerical experiments to compare the behaviors
of the locfdr and the $q$-value methods.

To ensure both the locfdr and the $q$-value methods utilize the equivalent
test statistics, we directly simulate $z^{2}$-statistics and obtain the
$p$-values for the $q$-value method according to the $\chi^{2}$-distribution
with 1 degree of freedom. We then compute $z_{i}={\rm
sgn}(z_{i})\sqrt{z_{i}^{2}}$, for each test $i$, where ${\rm sgn}(z_{i})$ is
independently drawn from the set $\\{+1,-1\\}$ with the probability 0.5 for
each choice. This transformation results in an overall symmetric $z$-score
distribution, and the corresponding two-sided $p$-values derived from the
$z$-statistics are identical to what is used in the $q$-value method. We
employ the two-groups model to simulate $z^{2}$ statistics with $\pi_{0}=0.55$
and $m=20,000$. The null data are simulated from the theoretical
$\chi^{2}_{1}$-distribution. The $z^{2}$ statistics of the alternative models
are generated from a family of $\Gamma(k,\theta)$ distributions, where $k$ and
$\theta$ denote the shape and scale parameters, respectively. We find that
this family of distributions is convenient to describe a wide spectrum of
alternative scenarios under both the UA and the ZA. In particular, if
$z\sim{\rm N}(0,\phi^{2})$, it follows that $z^{2}\sim\Gamma(0.5,2\phi^{2})$.
Adjusting the shape parameter $k$ in the range of $(0,0.5)$ results in heavy-
tailed distributions of $z$-scores that resemble $t$ and double-exponential
distributions. In comparison, setting $k\gg 0.5$ leads to the distribution of
the simulated $z$-scores with two modes that are symmetric to 0, and the
larger $k$ value results in greater separation of the two modes. Under this
setting, the ZA seems reasonable for $k>0.5$, otherwise, the generated data
are better described by the UA assumption. For inference, we apply the same
non-parametric estimator of $\pi_{0}$ (14) for both the locfdr and the
$q$-value methods by setting $\eta=0.5$ and assuming the true theoretical null
distribution of $z$-scores, i.e., ${\rm N}(0,1)$.

The simulation results show that when $k\geq 0.5$, the two methods always
yield highly concordant rejection paths, and the practical difference is
typically negligible. In comparison, as $k$ is decreased from 0.5 toward 0,
the rejection paths become noticeably different: although both approaches
behave conservatively, the locfdr method is significantly more conservative
than the $q$-value method. In light of Theorem 2, the observation indicates
that the ranking of the data becomes increasingly different by the two
approaches as $k$ decreases from 0.5 towards 0. We examine the empirical rank
correlations of $p$-values and the values of the inferred local fdr’s from the
simulated data with respect to a range of $k$ values. Figure 3 clearly
verifies the observed pattern. This can be explained by the fact that
likelihood ratio statistics $f_{1}(|z|)/f_{0}(|z|)$ may not be monotonically
decreasing with respect to $|z|$ for heavy-tailed distributions.

The most important finding from this exercise is that comparing the rejection
paths (or the rankings of test statistics) between the locfdr and the
$q$-value method may be helpful in diagnosing the zero assumption, as the
concordance of the rejections seemingly suggests the validity of the ZA. It is
worth noting that in most cases the divergence between the rejection paths
becomes noticeable only when the pre-defined control level is large. For
practically applied stringent FDR control levels (e.g., $\alpha=0.05$), the
two approaches generally agree well in all situations. (This is because the
ranking of data points at extreme tails are typically consistent in both
$p$-values and likelihood ratio statistics.) Nevertheless, the complete
rejection path, rather than the part under some stringent threshold, is
informative for validating the model assumptions.

|
---|---
Figure 3: Numerical comparison between the locfdr and the $q$-value methods.
The top-left panel shows the disconcordance of the rejection paths between the
two approaches as the ZA is severely violated (shape parameter set to 0.3 for
the alternative distribution). The top-right panel indicates the concorddance
of the two rejection paths when the ZA is reasonable (shape parameter $=1.0$
for the alternative distribution). The bottom panel shows the rank correlation
between $p$-values (used by the $q$-value method) and the local fdr’s computed
by the locfdr method. Clearly, as the ZA is severely violated, the rank
correlation significantly deviates from 1. Based on Theorem 2, this explaines
the disconcordance of the rejection paths between the locfdr and the $q$-value
methods

Finally, we apply the comparison to the Hedenfalk data distributed in the
$q$-value package. The data set includes 3,220 pre-computed $p$-values from a
differential gene expression study. We transform the provided $p$-values to
the corresponding $z$-scores according to the distribution function of the
standard normal and apply both the locfdr and the $q$-value methods using the
same $\hat{\pi}_{0}$ estimation method described above. The comparison of the
rejection paths indicates the results of the two approaches have an excellent
agreement (Figure 4).

Figure 4: comparison of rejection paths of the lfdr and the $q$-value methods
in Hedenfalk data. Although the number of tests is modest ($m=3170$), the two
rejection paths appear to be highly concordant.

###  4.2 Model Diagnosis for Bayesian FDR Control Procedures

We perform simulations to illustrate the importance of post-fitting model
diagnosis for Bayesian FDR control procedures, especially for the PEB-based
approaches. Adopting a similar simulation scheme described in the previous
section (i.e., simulating $z$ or $z^{2}$-statistics from various alternative
distributions including heavy-tailed, bi-modal distributions at the two
different extremes), we generate data from the two-group model with a spectrum
of alternative scenarios. Our aim is to examine the robustness of the PEB-
based Bayesian FDR control procedures under different modes of model
misspecification. We are also interested in investigating the ability of
standard model diagnosis techniques for detecting the inadequate model fitting
and assessing its consequence on FDR control.

For demonstration, we focus on a recently proposed Bayesian FDR control
method, ASH, a PEB implementation assuming the UA (Stephens, 2016). Briefly,
ASH, implemented in the R package ashr, has the ability to flexibly model uni-
modal alternative distributions from a mixture of a rich class of base
distributions. The weights for each basis distribution function are estimated
from data by an EM algorithm. To serve our purpose in this demonstration, we
intentionally limit the base distributions to the class of Normal
distributions with mean 0 and variable values of variance parameters, i.e.,
the alternative model assumes

$f_{1}=\sum_{k=1}^{K}\omega_{k}\,{\rm N}(0,\sigma_{k}^{2}),$

where $\omega_{k}$’s are the weight parameters to be estimated by the EM
algorithm and a grid of $\sigma_{k}$’s values are pre-selected in a data-
driven way (Stephens, 2016). It is also important to note that the Bayes
factors computed under the assumed alternative model are monotonic to the two-
sided $p$-values derived from the $z$-statistics. We anticipate that this
modeling strategy would be robust (although not always perfect) if the UA
assumption is indeed correct, but potentially problematic if the ZA actually
holds.

To perform model diagnosis, we compare the fitted mixture distribution,

$\hat{f}_{c}=\hat{\pi}_{0}f_{0}+(1-\hat{\pi}_{0})\sum_{k=1}^{K}\hat{\omega}_{k}\,{\rm
N}(0,\sigma_{k}^{2})$

with the observed mixture distribution from data. In particular, we compare a
full range of theoretical quantiles computed from $\hat{f}_{c}$ with the
observed sample quantiles. Because of the large sample size employed in the
simulation ($m=10,000$), we are able to assess an asymptotic $p$-value for
each pair of quantile values compared (Ferguson, 2017), which we use to
evaluate the severity of the mismatch between the fitted theoretical quantiles
and the observed ones. This strategy is indeed a standard approach to Bayesian
model diagnosis by noting that $\hat{f}_{c}$ is an approximation of posterior
predictive distribution as shown in Equation (12).

For each simulated data set, we also construct an oracle rejection path by
computing the likelihood ratios/Bayes factors assuming the true alternative
model and plugging in the true $\pi_{0}$. We also similarly construct an
expected rejection path by utilizing a two-sided $p$-value for each test data
point and the estimated value of $\hat{\pi}_{0}$ from the EM algorithm.
Comparing the rejection path from the Bayesian approach to the expected
rejection path should provide some assessment of goodness-of-fit according to
the theoretical result of Theorem 2.

First, we highlight the results from three distinct alternative distributions.

##### 1\. Scaled $t$-distribution with 10 degrees of freedom.

This alternative distribution follows the UA but modestly deviates from the
specific parametric assumption. Despite the imperfect model specification, the
model diagnosis indicates that the fit by ASH is mostly adequate (Table 1),
i.e., there are no obvious mismatches between the estimated and sample
quantiles. The comparison between the ASH rejection path and the oracle path
(Figure 5) indicates that the ASH result is slightly conservative, i.e., FDRs
are over-estimated by ASH, especially at more relaxed threshold values. The
$q$-value method behave much more conservatively, judging by the comparison of
its rejection path with the oracle path. Finally, we find that the ASH
rejection path is highly concordant to the expected rejection path, which
validates the result of model diagnosis and Theorem 2.

Table 1:  Model diagnosis for alternative scenario 1 ($t$-distribution)  | 10% | 20% | 30% | 40% | 50% | 60% | 70% | 80% | 90%
---|---|---|---|---|---|---|---|---|---
sample quantile | 0.051 | 0.214 | 0.495 | 0.962 | 1.710 | 2.955 | 5.116 | 8.811 | 16.81
fitted quantile | 0.052 | 0.214 | 0.509 | 0.983 | 1.728 | 2.928 | 4.974 | 8.723 | 16.76
$p$-value | 0.781 | 0.935 | 0.415 | 0.446 | 0.700 | 0.723 | 0.249 | 0.667 | 0.910

##### 2\. Scaled double-exponential distribution.

This alternative distribution follows the UA but has a much heavier tail than
the assumed alternative model could sufficiently capture. The model diagnosis
procedure can effectively detect the inadequate fitting in this case.
Importantly, we identify a pattern in all severe mismatches between the
estimated and sample quantiles: the theoretical quantiles of the fitted
distribution are consistently over-estimated (Table 2). Furthermore, there is
also noticeable discordance between the ASH rejection path and the expected
rejection path, indicating that the goodness-of-fit condition of Theorem 2 is
violated (Figure 5). Comparing to the oracle rejection path, ASH, local fdr
and $q$-value approaches all display significantly conservative behaviors,
which can be explained by the fact that all three methods over-estimate
$\pi_{0}$. However, ASH exhibits a less degree of conservativeness and much-
improved power over the existing the ZA-based approaches.

Table 2:  Model diagnosis for alternative scenario 2 (double exponential distribution) | 10% | 20% | 30% | 40% | 50% | 60% | 70% | 80% | 90%
---|---|---|---|---|---|---|---|---|---
sample quantile | 0.024 | 0.103 | 0.238 | 0.461 | 0.826 | 1.400 | 2.328 | 4.175 | 9.144
fitted quantile | 0.028 | 0.114 | 0.266 | 0.503 | 0.859 | 1.403 | 2.305 | 4.093 | 9.331
$p$-value | $0.048^{*}$ | $0.021^{*}$ | $0.001^{*}$ | $0.003^{*}$ | 0.129 | 0.938 | 0.672 | 0.447 | 0.523

##### 3\. Bimodal distribution induced by $\Gamma$ distribution with $k=0.7$

This alternative distribution departs from the UA and is more properly
described by the ZA. In this scenario, the model diagnosis procedure also
identifies the severely poor fit of the data. Importantly, we note that the
pattern of mismatch between the estimated and observed quantiles are very
different from the case of heavy tail alternative distributions: the
theoretical quantiles of the fitted distribution are consistently under-
estimated (Table 3). Furthermore, there is also noticeable discordance between
the ASH rejection path and the expected rejection path, indicating that the
goodness-of-fit condition in Theorem 2 is violated. In comparison to the
oracle rejection path, we note that ASH is severely under-estimate $\pi_{0}$,
which leads to dangerous anti-conservative behavior in FDR control, whereas
the $q$-value method behaves properly under this scenario (Figure 5).

Table 3:  Model diagnosis for alternative scenario 3 (bimodal distribution)  | 10% | 20% | 30% | 40% | 50% | 60% | 70% | 80% | 90%
---|---|---|---|---|---|---|---|---|---
sample quantile | 0.099 | 0.371 | 0.823 | 1.481 | 2.418 | 3.774 | 5.827 | 9.023 | 14.66
fitted quantile | 0.079 | 0.321 | 0.750 | 1.405 | 2.357 | 3.728 | 5.747 | 8.929 | 14.94
$p$-value | $<0.001^{*}$ | $<0.001^{*}$ | $0.003^{*}$ | $0.048^{*}$ | 0.281 | 0.570 | 0.480 | 0.570 | 0.279
|  |
---|---|---
|  |
|  |
Figure 5: evaluation and diagnosis of PEB-based Bayesian FDR control procedure
assuming UA. The inference is conducted by the ASH method. Each row represents
a unique alternative distribution: the $t$ and double exponential
distributions satisfy the UA (but not the particular parametric assumptions by
ASH), the bimodal alternative severely violates the UA assumption. The first
column shows the comparison of the observed rejection paths by ASH and the
oracle rejection paths for each alternative distribution; the second column
compares the observed rejection paths by ASH and the expected rejection paths;
and the third column provides the reference comparison between the rejection
paths of the $q$-value method and the oracle rejection paths.

Summarizing from these observations, we conclude that i) model diagnosis
method is essential and effective in validating Bayesian FDR control methods;
ii) even the data may be inadequately fitted by the assumed parametric model,
identification of specific patterns of misfitting can be helpful to determine
the conservativeness/anti-conservativeness of a method in FDR control. It is
worth emphasizing that specific patterns of misfitting typically depends on
the inference methods/parametric assumptions. In the case of ASH, the specific
mismatch patterns can be reasonably explained by the following observation:
with flexible base distributions, ASH typically provides very good fit at the
tails of the mixture distributions, where the data points have higher
influence on the overall (marginal) likelihood; this behavior, nevertheless,
causes inadequate fit for data of the non-tail areas and the specific patterns
can be detected by examining a full range of quantiles.

To verify this finding, we conduct additional simulations of $z$-scores using
the alternative models from a family of $\Gamma$ distributions with the shape
parameter varying from 0.1 to 0.9. For each shape parameter, we generate 100
data sets with $m=10,000$. We apply both ASH and $q$-value methods to compute
the rejection paths and perform the described model diagnosis procedures.
Figure 6 highlights the difference of $\pi_{0}$ estimation by the UA (ASH) and
the ZA ($q$-value) assumptions. Under the ZA, the $\hat{\pi}_{0}$ estimates
remain upper-bounds of the true value, but noticeably more conservative when
the UA assumption is true. On the other hand, we confirm that the UA
assumption can lead to severe under-estimation of $\pi_{0}$ when the ZA is
indeed appropriate. Finally, the simulation results indicate that the simple
model diagnosis method has good sensitivity and specificity in detecting the
anti-conservative FDR control behaviors by ASH. In particular, we flag a test
case if at least one theoretical quantile under-estimates the corresponding
sample quantile and the mismatch exceeds a pre-defined $p$-value threshold.
For a suggestive threshold of $p$-value $=0.05$, the model diagnosis procedure
flags 70.1% of the test cases where the UA assumption is violated (i.e., the
shape parameter $\geq 0.6$). In comparison, it only flags 4.2% of the test
cases when UA assumption is correctly specified (i.e., the shape parameters
$\leq 0.5$)

Figure 6:  Comparison of UA and ZA estimations of $\pi_{0}$ in different
settings. The true $\pi_{0}$, indicated by the dotted horizontal line, is
0.60. The ZA estimates are obtained by the nonparametric quantile estimator
(14), and the UA estimates are obtained by the EM algorithm assuming a mixture
normal alternative distribution. The ZA estimates are uniformly conservative
in all scenarios. The UA estimates, in comparison, are more accurate when the
assumption holds (i.e., shape parameter $\leq 0.5$); but can be significantly
under-estimated if the assumption is severely violated (i.e., shape parameter
$>0.5$.

##  5 Conclusion and discussion

Our results described in this paper provide a theoretical connection between
Bayesian and frequentist FDR control methods in the presence of a large number
of tests. The setting of large-scale multiple testing creates a near
asymptotic scenario where the frequentist property of the Bayesian method can
be practically examined. Our results also have profound implications for
applying parametric Bayesian FDR control methods in a practical setting, in
particular, the critical importance of post-fitting model diagnosis. From our
theoretical arguments and numerical experiments, we conclude that thorough and
careful model checking can indeed improve the robustness of the parametric
Bayesian FDR control procedures. More specifically, we have showcased the
numerical recipes in examining both the ZA (by examning the concordance of
rejection paths between the locfdr and the $q$-value methods) and the UA (by
Bayesian model diagnosis) assumptions. Beyond the settings of the two-groups
model, we argue that the Bayesian approaches show superior advantages in
extensibility and flexibility over the $p$-value based frequentist approaches.
In fact, we have shown that the ideal/optimal frequentist procedure under such
setting should very much resemble the corresponding Bayesian approach.

It has been previously shown in the literature (e.g., Sun and Cai (2007)) that
Bayesian FDR control procedures can be generalized to a broad class of binary
probabilistic classification problem for prioritizing the control of one
specific type of misclassification errors, i.e., the false positives. From
this viewpoint, hypothesis testing is simply a special case with null
hypotheses being framed as non-discoveries, and as a generalization, applying
Bayesian procedures help overcome some practical limitations associated with
formulating and solving hypothesis testing problems. Here we demonstrate this
point in two cases.

##### Case 1: testing complex null hypothesis

In standard hypothesis testing practice, a key step is to set up a null
hypothesis with a well-defined $f_{0}$ distribution. With the emerging complex
data structures in applications, this can be difficult as some ideal null
models may involve very complex scenarios. An example of this kind is to test
colocalization in genetic association studies where the scientific goal is to
statistically assess if a genetic variant is simultaneously associated two
separate complex traits. Let latent binary indicators $\gamma_{i}$ and
$\delta_{i}$ denote the association status of the $i$-th genetic variant with
trait 1 and trait 2, respectively. In the formulation of hypothesis testing
problem, the goal is to test

$\displaystyle H_{0}:\gamma_{i}\neq 1\mbox{ or }\delta_{i}\neq 1,$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}vs.$ $\displaystyle
H_{1}:\gamma_{i}=1\mbox{ and }\delta_{i}=1.$

Although the marginal associations with respect to $\gamma_{i}$ and
$\delta_{i}$ are relatively easy to test, the null hypothesis for
colocalization involves multiple compatible and mutually exclusive scenarios.
Consequently, the corresponding $f_{0}$ distribution is a mixture distribution
for which $p$-value is difficult to compute or calibrate. In comparison, the
Bayesian approach only needs to focus on evaluating
$\Pr(\gamma_{i}=1,\delta_{i}=1\mid\mbox{ data})$ for the alternative scenario,
and the quantity, $1-\Pr(\gamma_{i}=1,\delta_{i}=1\mid\mbox{ data})$, is
indeed the desired local fdr. (A practical Bayesian solution for the
colocalization problem can be found in Wen et al. (2017).) Another related
example is the estimation of the empirical null distribution in the local fdr
framework Efron (2012). Similarly, such idea is extremely difficult, if not
impossible, to implement in the frequentist FDR control framework.

##### Case 2: flexible definition of discovery

The Bayesian framework allows the flexible definition of discoveries in
probabilistic terms, which can lead to more interpretable results for better
scientific communication. A longstanding scientific problem that motivates the
development of FDR control technique is the analysis of microarray and RNA-seq
data analysis for identifying differentially expressed genes in two
experimental conditions. The common practice in this context is to test a
“sharp null” hypothesis which states that the difference of the expression
levels of a target gene from the two conditions is exactly 0. Nevertheless,
many have pointed out the caveat of the sharp null hypothesis (e.g., it is
scientifically plausible that all genes are differentially expressed under the
different conditions) and argued that the real purpose of the analysis is to
identify genes with the convincingly large non-zero difference in expression
levels, denoted by $\beta$. Recently, Stephens (2016) derives the local false
sign rate (lfsr) as the posterior error probability in estimating the sign of
$\beta$, which enables probabilistic classification of large versus small
effects and a straightforward application of Bayesian FDR control procedure.
In another example, Li et al. (2011) defines a discovery in high-throughput
genomic experiments if a noticeable observation is replicable in a follow-up
validation study. They derive a probabilistic measure, irreproducible
discovery rate (IDR), and apply the Bayesian FDR control procedure to identify
those convincingly reproducible signals. In both examples, there seem no
straightforward applications of the frequentist FDR control procedure.

Finally, we conclude that our work in this paper should motivate more
applications of Bayesian FDR control approaches involving large-scale and
complex data: the Bayesian approaches can be flexible, powerful, robust, and
most importantly, verifiable by examining their frequentist properties.

##  6 Resources

The code implementing the illustrations and simulations used in this paper can
be found in the github repo: https://github.com/xqwen/unified_fdr.

##  Appendix A Proof of Theorem 1

###### Proof.

By equation (4),

${\rm
BFDR}\left(t_{b,i}=u^{*}_{(i)}\right)=\frac{\left(\sum_{k=1}^{m}u_{k}^{*}\cdot{\bf
1}\\{u^{*}_{k}\leq u^{*}_{(i)}\\}\right)\big{/}m}{i/m}.$ (20)

Note that, individually, because $u_{k}^{*}(z_{k})={\rm E}({\bf
1}\\{\gamma_{k}=0\\}\mid\pi_{0},z_{k})$, it follows that

${\rm E}\left[u_{k}^{*}\cdot{\bf 1}{\\{u_{k}^{*}\leq
u_{(i)}^{*}\\}}\,\bigl{|}\,\pi_{0},u_{(i)}^{*}\right]={\rm E}\left[{\bf
1}{\\{u_{k}^{*}\leq u_{(i)}^{*}\mbox{ and
}\gamma_{k}=0\\}}\,\bigl{|}\,\pi_{0},u_{(i)}^{*}\right]=\pi_{0}\,F_{u^{*}}(u_{(i)}^{*}),\forall~{}k,$
(21)

where $F_{u^{*}}(\cdot)$ denotes the CDF of $u^{*}$’s when data are generated
from $f_{0}$. It then follows directly from the strong law of large number
(SLLN),

$\frac{\left(\sum_{k=1}^{m}u_{k}^{*}\cdot{\bf 1}\\{u^{*}_{k}\leq
u^{*}_{(i)}\\}\right)\big{/}m}{i/m}\xrightarrow{a.s.}\frac{m\pi_{0}F_{u^{*}}(u^{*}_{(i)})}{i}~{}\mbox{
for all $i$}.$

Note, Equation (8) indicates that $u^{*}(z)$ is a monotone function of ${\rm
BF}^{*}(z)$, and implies that there exists a monotone mapping, denoted by
$G(\cdot)$, between $u^{*}$ and the corresponding $p$-value by the assumption
2. As a consequence, $G(u^{*}_{(i)})=p_{(i)}$. It follows that

$F_{u^{*}}(u^{*}_{(i)})=\Pr\left[u^{*}\leq u^{*}_{(i)}\mid
H_{0}\right]=\Pr\left[G(u^{*})\leq G(u^{*}_{(i)})\mid H_{0}\right]=p_{(i)}.$

Finally, by equation (7) we conclude that

${\rm BFDR}\left(t_{b,i}=u^{*}_{(i)}\right)\xrightarrow{a.s.}{\rm
FDR}\left(t_{f,i}=p_{(i)}\right)~{}\mbox{for all $i$}.$

∎

##  Appendix B Proof of Theorem 2

###### Proof.

The proof is built upon the same technical arguments used in Sun and Cai
(2007) (SC in brief, henceforward).

We first define

$\tilde{u}(z):=\frac{\hat{\pi}_{0}f_{0}(z)}{f_{c}(z)}={\rm E}\left({\bf
1}\\{\gamma_{i}=0\\}\mid\hat{\pi}_{0},z\right).$

Given that ${\rm E}||\hat{f}_{c}-f_{c}||^{2}\to 0$, it can be shown, by Lemma
A.1 and Lemma A.2 of SC, that

$\hat{u}\xrightarrow{p}\tilde{u},\forall i,$ (22)

and consequently, for any given order statistics from observed sequences

$\hat{u}_{(k)}\xrightarrow{p}\tilde{u}_{(k)},\forall k,$ (23)

Equation (22) and the Helly-Bray theorem imply

$g(t)={\rm E}\left[\hat{u}_{k}\cdot{\bf 1}{\\{\hat{u}_{k}\leq
t\\}}\,\bigl{|}\,\hat{\pi}_{0},t\right]\to{\rm E}\left[\tilde{u}_{k}\cdot{\bf
1}{\\{\tilde{u}_{k}\leq
t\\}}\,\bigl{|}\,\hat{\pi}_{0},t\right]=h(t),~{}\forall t\in[0,1].$ (24)

Furthermore by Equation (23),
$g(\hat{u}_{(k)})\xrightarrow{p}h(\tilde{u}_{(k)})$, i.e.,

${\rm E}\left[\hat{u}_{k}\cdot{\bf
1}{\\{\hat{u}_{k}\leq\hat{u}_{(i)}\\}}\,\bigl{|}\,\hat{\pi}_{0},\hat{u}_{(i)}\right]\xrightarrow{p}{\rm
E}\left[\tilde{u}_{k}\cdot{\bf
1}{\\{\tilde{u}_{k}\leq\tilde{u}_{(i)}\\}}\,\bigl{|}\,\hat{\pi}_{0},\tilde{u}_{(i)}\right],~{}\forall~{}k.$
(25)

Next, we show that

${\rm E}\left[\hat{u}_{k}\cdot{\bf
1}{\\{\hat{u}_{k}\leq\hat{u}_{(i)}\\}}\,\bigl{|}\,\hat{\pi}_{0},\hat{u}_{(i)}\right]\xrightarrow{p}\hat{\pi}_{0}p_{(i)}.$
(26)

This is because

$\displaystyle{\rm E}\left[\tilde{u}_{k}\cdot{\bf
1}{\\{\tilde{u}_{k}\leq\tilde{u}_{(i)}\\}}\,\bigl{|}\,\hat{\pi}_{0},\tilde{u}_{(i)}\right]$
$\displaystyle={\rm E}\left({\bf 1}\\{\tilde{u}_{k}\leq\tilde{u}_{(i)}\mbox{
and }\gamma_{k}=0\\}\mid\hat{\pi}_{0},\tilde{u}_{(i)}\right)$ (27)
$\displaystyle=\hat{\pi}_{0}F_{\tilde{u}}\left(\tilde{u}_{(i)}\right),$

where $F_{\tilde{u}}\left(\cdot\right)$ denotes the cdf of $\tilde{u}(z)$, and
similarly we denote $F_{\hat{u}}\left(\cdot\right)$ as the cdf of $\hat{u}(z)$
under the null. Note that, by (22) and (23),

$F_{\tilde{u}}\left(\tilde{u}_{(i)}\right)={\rm E}\left({\bf
1}\\{\tilde{u}\leq\tilde{u}_{(i)}\\}\mid
H_{0},\tilde{u}_{(i)}\right)\xrightarrow{p}{\rm E}\left({\bf
1}\\{\hat{u}\leq\hat{u}_{(i)}\\}\mid
H_{0},\tilde{u}_{(i)}\right)=F_{\hat{u}}\left(\hat{u}_{(i)}\right)$ (28)

Under the assumption ${\rm BF}(z)=\frac{f_{0}(z)}{\hat{f}_{1}(z)}$ is
monotonic to the corresponding $p$-value, it follows that for any given
$\hat{\pi}_{0}$, there exists a monotone transformation, $\hat{G}(\cdot)$,
between $\hat{u}$ and its corresponding $p$-value, such that
$\hat{G}(\hat{u})=p$, and $\hat{G}(\hat{u}_{(i)})=p_{(i)}$ for any given order
statistics. It follows that,

$F_{\hat{u}}\left(\hat{u}_{(i)}\right)=\Pr\left[\hat{u}\leq\hat{u}_{(i)}\mid
H_{0}\right]=\Pr\left[\hat{G}(\hat{u})\leq\hat{G}(\hat{u}_{(i)})\mid
H_{0}\right]=p_{(i)}.$ (29)

Jointly by (25), (27), (28) and (29), we have established (26), i.e.,

${\rm E}\left[\hat{u}_{k}\cdot{\bf
1}{\\{\hat{u}_{k}\leq\hat{u}_{(i)}\\}}\,\bigl{|}\,\hat{\pi}_{0},\hat{u}_{(i)}\right]\xrightarrow{p}\hat{\pi}_{0}p_{(i)}.$

Finally, note that

${\rm
BFDR}\left(t_{b,i}=\hat{u}_{(i)};\hat{\pi}_{0}\right)\xrightarrow{a.s.}\frac{{\rm
E}\left[\hat{u}_{k}\cdot{\bf
1}{\\{\hat{u}_{k}\leq\hat{u}_{(i)}\\}}\,\bigl{|}\,\hat{\pi}_{0}\right]}{i/m},$
(30)

and

${\rm
FDR}\left(t_{f,i}=p_{(i)};\hat{\pi}_{0}\right)\xrightarrow{a.s.}\frac{\hat{\pi}_{0}p_{(i)}}{i/m},$
(31)

we therefore conclude that

${\rm BFDR}\left(t_{b,i}=\hat{u}_{(i)};\hat{\pi}_{0}\right)\xrightarrow{p}{\rm
FDR}\left(t_{f,i}=p_{(i)};\hat{\pi}_{0}\right),\forall i.$ (32)

In addition, if the condition $\hat{\pi}_{0}\xrightarrow{p}\pi_{0}$ is
satisfied, i.e., $\pi_{0}$ is accurately estimated by $\hat{\pi}_{0}$. It is
straightforward to show that

${\rm
BFDR}\left(t_{b,i}=\hat{u}_{(i)};\hat{\pi}_{0}\right)\xrightarrow{p}\frac{\pi_{0}p_{(i)}}{i/m}.$
(33)

By Theorem 1, it is sufficient to conclude

${\rm BFDR}\left(t_{b,i}=\hat{u}_{(i)};\hat{\pi}_{0}\right)\xrightarrow{p}{\rm
BFDR}(t_{b,i}=u^{*}_{(i)}),\forall i.$ (34)

∎

##  Appendix C Proof of Proposition 1

###### Proof.

Note that

${\rm
BFDR}\left(t_{b,i}=u^{*}_{(i)}\right)=\frac{\left(\sum_{k=1}^{m}u_{k}^{*}\cdot{\bf
1}\\{u^{*}_{k}\leq u^{*}_{(i)}\\}\right)\big{/}m}{i/m}.$

Individually,

${\rm E}\left[u_{k}^{*}\cdot{\bf 1}_{\\{u_{k}^{*}\leq
u_{(i)}^{*}\\}}\,\bigl{|}\,d_{i}\right]={\rm E}\left[{\bf
1}_{\\{wlr^{*}_{k}\geq wlr_{(i)}^{*}\mbox{ and
}\gamma_{i}=0\\}}\,\bigl{|}\,d_{i}\right]=\pi_{d_{i},0}\,[1-F_{d_{i},wlr^{*}}(u_{(i)}^{*})],$

for all $k$, and

${\rm Var}\left[{\bf 1}_{\\{wlr^{*}_{k}\geq wlr_{(i)}^{*}\mbox{ and
}\gamma_{i}=0\\}}\,\bigl{|}\,d_{i}\right]\mbox{ is obviously bounded.}$

It follows from the Kolmogorov’s SSLN that

${\rm BFDR}\left(t_{b,i}=u^{*}_{(i)}\right)\xrightarrow{a.s.}{\rm
FDR}\left(t_{f,i}=wlr_{(i)}^{*}\right)\mbox{ for all }i.$

∎

## References

  * Benjamini and Hochberg (1995) Y. Benjamini and Y. Hochberg.  Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , 57(1):289–300, 1995.
  * Cook et al. (2006) S. R. Cook, A. Gelman, and D. B. Rubin.  Validation of software for bayesian models using posterior quantiles.  _Journal of Computational and Graphical Statistics_ , 15(3):675–692, 2006.
  * Cox and Hinkley (1979) D. R. Cox and D. V. Hinkley.  _Theoretical statistics_.  CRC Press, 1979.
  * Efron (2008) B. Efron.  Microarrays, Empirical Bayes and the Two-Groups Model.  _Statistical Science_ , 23(1):1–22, Feb. 2008\.
  * Efron (2012) B. Efron.  _Large-scale inference: empirical Bayes methods for estimation, testing, and prediction_ , volume 1.  Cambridge University Press, 2012.
  * Efron et al. (2001) B. Efron, R. Tibshirani, J. D. Storey, and V. Tusher.  Empirical bayes analysis of a microarray experiment.  _Journal of the American Statistical Association_ , 96(456):1151–1160, 2001.
  * Ferguson (2017) T. S. Ferguson.  _A course in large sample theory_.  Routledge, 2017.
  * Genovese and Wasserman (2004) C. Genovese and L. Wasserman.  A stochastic process approach to false discovery control.  _Annals of Statistics_ , 32(3):1035–1061, 2004\.
  * Genovese et al. (2006) C. R. Genovese, K. Roeder, and L. Wasserman.  False discovery control with p-value weighting.  _Biometrika_ , 93(3):509–524, 2006.
  * Ghosh et al. (2006) D. Ghosh, W. Chen, and T. Raghunathan.  The false discovery rate: a variable selection perspective.  _Journal of Statistical Planning and Inference_ , 136(8):2668–2684, Aug. 2006.
  * Ignatiadis et al. (2016) N. Ignatiadis, B. Klaus, J. B. Zaugg, and W. Huber.  Data-driven hypothesis weighting increases detection power in genome-scale multiple testing.  _Nature methods_ , 13(7):577–580, 2016.
  * Kass and Steffey (1989) R. E. Kass and D. Steffey.  Approximate bayesian inference in conditionally independent hierarchical models (parametric empirical bayes models).  _Journal of the American Statistical Association_ , 84(407):717–726, 1989.
  * Li et al. (2011) Q. Li, J. B. Brown, H. Huang, P. J. Bickel, et al.  Measuring reproducibility of high-throughput experiments.  _The annals of applied statistics_ , 5(3):1752–1779, 2011.
  * Müller et al. (2004) P. Müller, G. Parmigiani, C. Robert, and J. Rousseau.  Optimal Sample Size for Multiple Testing: The Case of Gene Expression Microarrays.  _Journal of the American Statistical Association_ , 99(468):990–1001, Dec. 2004.
  * Müller et al. (2006) P. Müller, G. Parmigiani, and K. Rice.  FDR and Bayesian multiple comparisons rules.  In _Bayesian Statistics 8_ , volume 0, pages 349–370. Oxford University Press, 2006.
  * Muralidharan (2010) O. Muralidharan.  An empirical Bayes mixture method for effect size and false discovery rate estimation.  _The Annals of Applied Statistics_ , 4(1):422–438, Mar. 2010.
  * Newton et al. (2004) M. A. Newton, A. Noueiry, D. Sarkar, and P. Ahlquist.  Detecting differential gene expression with a semiparametric hierarchical mixture method.  _Biostatistics_ , 5(2):155–76, Apr. 2004.
  * Stephens (2016) M. Stephens.  False discovery rates: a new deal.  _Biostatistics_ , 18(2):275–294, 2016.
  * Storey (2002) J. D. Storey.  A direct approach to false discovery rates.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 64(3):479–498, Aug. 2002.
  * Storey (2003) J. D. Storey.  The positive false discovery rate: a Bayesian interpretation and the q-value.  _The Annals of Statistics_ , 31(6):2013–2035, Dec. 2003.
  * Storey (2007) J. D. Storey.  The optimal discovery procedure: a new approach to simultaneous significance testing.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 69(3):347–368, June 2007.
  * Storey et al. (2004) J. D. Storey, J. E. Taylor, and D. Siegmund.  Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 66(1):187–205, Feb. 2004.
  * Sun and Cai (2007) W. Sun and T. T. Cai.  Oracle and adaptive compound decision rules for false discovery rate control.  _Journal of the American Statistical Association_ , 102(479):901–912, 2007.
  * Wakefield (2009) J. Wakefield.  Bayes factors for genome-wide association studies: comparison with P-values.  _Genetic Epidemiology_ , 33(1):79–86, 2009.
  * Wen et al. (2017) X. Wen, R. Pique-Regi, and F. Luca.  Integrating molecular qtl data into genome-wide genetic association analysis: Probabilistic assessment of enrichment and colocalization.  _PLoS genetics_ , 13(3):e1006646, 2017.
  * Whittemore (2007) A. S. Whittemore.  A Bayesian False Discovery Rate for Multiple Testing.  _Journal of Applied Statistics_ , 34(1):1–9, Jan. 2007.