# Generalized estimators for multiple testing: proportion of true nulls and
false discovery rate

Xiongzhi Chen  Princeton University, Princeton, USA.  R.W. Doerge
††footnotetext: Address for correspondence: R.W. Doerge, Department of
Statistics, Purdue University, West Lafayette, IN 47907-1399, USA.
Email<EMAIL_ADDRESS>Purdue University, West Lafayette, USA.

###### Abstract

Two new estimators are proposed: the first for the proportion of true null
hypotheses and the second for the false discovery rate (FDR) of one-step
multiple testing procedures (MTPs). They generalize similar estimators
developed jointly by Storey, Taylor and Siegmund and can be applied also to
discrete $p$-values whose null distributions dominate the uniform
distribution. For the new estimator of the FDR, we establish its simultaneous
asymptotic conservativeness and justify formally the stopping time property of
its threshold for $p$-values not necessarily independent or continuous. Our
empirical studies show that, when applied to the aforementioned $p$-values,
both of our estimators usually outperform their competitors (considered in
this work) except that the first one may under-estimate the target when the
needed tuning parameters are inappropriately chosen. The methodology of our
work easily extends to other types of discrete $p$-values. We illustrate the
improvement in power our estimators can induce by applying them to next-
generation sequencing count data.

###### keywords:

Discrete $p$-values; False discovery rate; Generalized estimators; Multiple
testing; Stopping time property

##  1 Introduction

In a typical multiple testing (MT) scenario, there are $m$ null hypotheses
$H_{i}$ with associated $p$-values $p_{i}$, $i=1,...,m$, such that only
$m_{0}$ of them are true nulls but the rest $m_{1}$ false nulls. Further, the
proportion $\pi_{0}=m_{0}/m$ is unknown. To better balance the overall
statistical error and power in MT, Benjamini and Hochberg (1995) proposed the
concept of the false discovery rate (FDR) and a linear step-up MTP (BH
procedure) whose FDR is no larger than a prespecified level. Let $I_{0}$ be
the set of indices $i$ for the true nulls and $I_{1}$ that for the false
nulls. A one-step MTP that uses a threshold $t\in\left[0,1\right]$ to decide
the status of each null hypothesis such that $H_{i}$ is false if and only if
$p_{i}\leq t$ induces

$\left\\{\begin{array}[]{l}V\left(t\right)=\\#\left\\{i\in I_{0}:p_{i}\leq
t\right\\}\text{,}\\\ S\left(t\right)=\\#\left\\{i\in I_{1}:p_{i}\leq
t\right\\}\text{,}\\\
R\left(t\right)=V\left(t\right)+S\left(t\right)\text{,}\end{array}\right.$

and its FDR is defined as

$FDR\left(t\right)=E\left[\left.\dfrac{V\left(t\right)}{R\left(t\right)}\right|R\left(t\right)>0\right]\text{.}$

Theorem 5.1 of Benjamini and Yekutieli (2001) shows that the FDR of the BH
procedure, $FDR_{BH}$, satisfies

$FDR_{BH}\leq\pi_{0}\alpha$ (1)

for independent $p$-values, where $\alpha\in\left[0,1\right]$ is the
prespecified FDR level. The upper bound in (1) shows the role of the unknown
$\pi_{0}$ in quantifying $FDR_{BH}$. Specifically, using an estimator
$\hat{\pi}_{0}$ of $\pi_{0}$, such that $\pi_{0}\leq\hat{\pi}_{0}<1$ when
$\pi_{0}<1$, gives a tighter upper bound on $FDR_{BH}$ than using $\pi_{0}=1$,
and this in turn yields a rejection threshold no smaller than that induced by
the BH procedure at the same FDR level. Consequently, incorporation of a
proper $\hat{\pi}_{0}$ potentially increases the power of the BH procedure.
This observation has produced various estimators (denoted by $\hat{\pi}_{0}$)
of $\pi_{0}$ and resulted in many adaptive FDR control and/or estimation
procedures (Storey et al.; 2004; Benjamini et al.; 2006; Blanchard and
Roquain; 2009). However, almost all $\hat{\pi}_{0}$’s and adaptive FDR
procedures are developed for continuous $p$-values whose cumulative
distribution functions (c.d.f.’s) under true nulls (i.e., null $p$-value
distributions) are uniform on the unit interval.

Besides continuous $p$-values, discrete $p$-values from various other tests,
such as the Fisher’s exact test (FET), Binomial test (BIN) (Lehmann and
Romano; 2005) and the exact Negative Binomial test (ENT; Robinson and Smyth;
2008 and Di et al.; 2011), are widely used when (conditional) inference is
conducted with discrete test statistics in the analysis of data produced by
next-generation sequencing (NGS) technologies and other sources. The null
$p$-value distributions induced by these tests are discrete and heterogeneous.
Without proper adjustment, these characteristics are very likely to make more
conservative $\hat{\pi}_{0}$ and FDR procedures that were originally designed
for continuous $p$-values. Pounds and Cheng (2006) discusses the effects of
discrete, non-uniform $p$-value distributions on $\hat{\pi}_{0}$ and FDR
procedures. For permutation $p$-values, they showed that the BH procedure and
the procedure in Storey (2002) can be erroneous, and proposed two robust
estimators $\hat{\pi}_{0}^{P}$ and $\tilde{\pi}_{0}^{P}$ of $\pi_{0}$ whose
resultant FDR estimators outperform the previous two.

Despite these results and the ubiquity of discrete data, little progress has
been made in developing better $\hat{\pi}_{0}$’s for discrete $p$-values. In
Section 2, we show that Storey’s estimator $\hat{\pi}_{0}^{S}$ in Storey
(2002), $\hat{\pi}_{0}^{P}$ and $\tilde{\pi}_{0}^{P}$ all have excessive
upward biases when applied to certain types of discrete $p$-values. In Section
3, we develop a new estimator $\hat{\pi}_{0}^{G}$ of $\pi_{0}$ for discrete
$p$-values such that $\hat{\pi}_{0}^{G}$ inherits all excellent properties of
$\hat{\pi}_{0}^{S}$ for continuous $p$-values. Unlike existing estimators of
$\pi_{0}$, $\hat{\pi}_{0}^{G}$ explicitly accounts for the deviations of
discrete $p$-value distributions from the uniform distribution and thus has
less upward bias. We provide conditions to ensure that $\hat{\pi}_{0}^{G}$ is
less conservative than $\hat{\pi}_{0}^{S}$ and $\tilde{\pi}_{0}^{P}$ and
discuss the choice of tuning parameters for $\hat{\pi}_{0}^{G}$. In Section 4,
we propose a generalized FDR estimator for one-step MTPs, establish its
simultaneous asymptotic conservativeness, and show the stopping time property
of the threshold of the generalized FDR estimator. In Section 5, we compare
the generalized estimators with some of their popular counterparts using two-
sided $p$-values from FETs, BINs and ENTs using simulated data. Our simulation
studies demonstrate the improved performance of the new estimators in terms of
accuracy, power and stability. However, we caution that $\hat{\pi}_{0}^{G}$
can under-estimate $\pi_{0}$ when its tuning parameters are not adaptively
chosen. The generalized estimators are applied to NGS count data for
Arabidopsis thaliana and shown to result in improvement. We conclude the paper
with a discussion in Section 6 and relegate all proofs into Appendix A.

##  2 Excessive bias of three existing estimators

Let the probability space be $\left(\Omega,\mathcal{A},P\right)$ with $\Omega$
the sample space, $\mathcal{A}$ the sigma-algebra on $\Omega$, and $P$ the
probability measure on $\mathcal{A}$. Without loss of generality (WLOG), we
assume $0<\pi_{0}<1$. The conservativeness of adaptive FDR procedures relies
on a crucial property of an estimator $\hat{\pi}_{0}$ of $\pi_{0}$, i.e.,
$E\left[\hat{\pi}_{0}\right]\geq\pi_{0}$ or $\hat{\pi}_{0}\geq\pi_{0}$,
$P$-a.s., called “conservativeness” or “total conservativeness” of the
estimator, respectively. Without (total) conservativeness, an adaptive FDR
procedure may fail to control or conservatively estimate the FDR of an MTP. We
show that some existing estimators of $\pi_{0}$ can be too conservative (i.e.,
too upwardly biased) for discrete $p$-values.

To estimate $\pi_{0}$, Storey (2002) introduces

$\hat{\pi}_{0}^{S}\left(\lambda\right)=\left(1-\lambda\right)^{-1}m^{-1}\sum\limits_{i=1}^{m}\mathbf{1}_{\left\\{p_{i}>\lambda\right\\}}\text{,
}\lambda\in[0,1)\text{.}$ (2)

Let $U\left(\tilde{a},\tilde{b}\right)$ denote the uniform distribution on a
non-empty interval $\left[\tilde{a},\tilde{b}\right]$. When $p_{i}\sim
U\left(0,1\right)$ for all $i\in I_{0}$,

$E\left[\hat{\pi}_{0}^{S}\left(\lambda\right)\right]=\pi_{0}+\left(1-\lambda\right)^{-1}m^{-1}\sum_{i\in
I_{1}}\left\\{1-F_{i}\left(\lambda\right)\right\\}=\pi_{0}+b_{1}$

and $b_{1}\geq 0$ (usually positive) is the bias caused by the unknown
$F_{i}$, $i\in I_{1}$, where $F_{i}$, $i=1,..,m$, is the c.d.f. of $p_{i}$ and
any c.d.f. is taken to be right continuous with left limits (i.e., càdlàg).
However, when $F_{i}\left(t\right)\leq t$ for $t\in\left[0,1\right]$ Lebesgue
almost surely (denoted by “$F_{i}\preceq U\left(0,1\right)$” and referred to
as “$F_{i}$ dominates the uniform distribution”) for $i\in I_{0}$,

$E\left[\hat{\pi}_{0}^{S}\left(\lambda\right)\right]=\pi_{0}+\left(1-\lambda\right)^{-1}m^{-1}\sum\limits_{i\in
I_{0}}\left\\{\lambda-
F_{i}\left(\lambda\right)\right\\}+b_{1}=\pi_{0}+b_{0}+b_{1}$

and $b_{0}\geq 0$ (usually positive) is the extra bias.

Pounds and Cheng (2006) justified that, for randomized tests, FDR procedures
that do not account for the discreteness or non-uniformity of the $p$-values
are usually unreliable. They proposed two estimators of $\pi_{0}$ as

$\tilde{\pi}_{0}^{P}=\left\\{\begin{array}[]{ccl}\min\left(1,2\bar{p}\right)&&\text{for
2-sided $p$-values,}\\\ \min\left(1,2\bar{a}\right)&&\text{for 1-sided,
continuous $p$-values,}\\\ \min\left(1,8\bar{a}\right)&&\text{for 1-sided,
discrete $p$-values,}\end{array}\right.$

where $\bar{p}=m^{-1}\sum\limits_{i=1}^{m}p_{i}$,
$a_{i}=2\min\left\\{p_{i},1-p_{i}\right\\}$,
$\bar{a}=m^{-1}\sum\limits_{i=1}^{m}a_{i}$, and

$\hat{\pi}_{0}^{P}=\min\left(1,\dfrac{1}{m}\sum_{i=1}^{m}\dfrac{p_{i}}{E\left[p_{i}|H_{i}=0\right]}\right)\text{,}$

where $H_{i}=0$ or $1$ is set depending on whether $H_{i}$ is a true or false
null. They claimed that $\tilde{\pi}_{0}^{P}$ is robust to various types of
$p$-value distributions and demonstrated its improvement over
$\hat{\pi}_{0}^{S}\left(\lambda\right)$. They further claimed
$\hat{\pi}_{0}^{P}$ to be better than $\tilde{\pi}_{0}^{P}$. Interestingly,
these two estimators still may have excessive upward bias when $F_{i}\preceq
U\left(0,1\right)$ for $i\in I_{0}$ since none of them adjusts for the
discreteness of $p$-values. WLOG, we consider two-sided discrete $p$-values
and assume neither of $\tilde{\pi}_{0}^{P}$ nor $\hat{\pi}_{0}^{P}$ is $1$.
For $\tilde{\pi}_{0}^{P}$, we see

$E\left[\tilde{\pi}_{0}^{P}\right]\geq 2m^{-1}\sum_{i\in
I_{0}}E\left[p_{i}\right]+2m^{-1}\sum_{i\in
I_{1}}E\left[p_{i}\right]\geq\pi_{0}+2m^{-1}\sum_{i\in
I_{1}}E\left[p_{i}\right]\text{.}$

So $\tilde{\pi}_{0}^{P}$ has extra (usually positive) bias
$b_{2}=2m^{-1}\sum_{i\in I_{0}}E\left[p_{i}\right]-\pi_{0}\geq 0$ besides the
bias $2m^{-1}\sum_{i\in I_{1}}E\left[p_{i}\right]>0$. Similarly,

$E\left[\hat{\pi}_{0}^{P}\right]=\pi_{0}+\dfrac{1}{m}\sum_{i\in
I_{1}}\dfrac{E\left[p_{i}|H_{i}=1\right]}{E\left[p_{i}|H_{i}=0\right]}=\pi_{0}+b_{3}\text{,}$

and if some ratios in the summation are large, so will be its positive bias
$b_{3}$.

##  3 A generalized estimator of $\pi_{0}$

Let $F_{i}^{\ast}$ denote the null $p$-value distribution for $p_{i}$,
$i=1,...,m$. Unless otherwise stated, it is assumed that $F_{i}^{\ast}$
dominates $U\left(0,1\right)$. To accommodate the discreteness of $p$-values
and reduce the upward bias of $\hat{\pi}_{0}$ as much as possible while
preserving its conservativeness, we define

$\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)=\max\left[0,\min\left\\{1,\hat{\pi}_{0}^{S}\left(\lambda\right)-\dfrac{\delta\left(\lambda,\mathbf{\epsilon}\right)}{\left(1-\lambda\right)m}\right\\}\right]\text{,}$
(3)

where

$\delta\left(\lambda,\mathbf{\epsilon}\right)=\sum_{i=1}^{m}\varepsilon_{i}\left\\{\lambda-
F_{i}^{\ast}\left(\lambda\right)\right\\}\text{,}$

and $\mathbf{\epsilon}=\left(\varepsilon_{1},...,\varepsilon_{m}\right)$ is a
vector of prespecified or adaptively estimated constants (usually in
$\left[0,1\right]$).

We denote $\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)$ by
$\hat{\pi}_{0}^{G}$. When $\varepsilon_{i}=\varepsilon$ for some
$\varepsilon\in\left[0,1\right]$ for all $1\leq i\leq m$, we write
$\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)$ as
$\hat{\pi}_{0}^{G}\left(\lambda,\varepsilon\right)$ and
$\hat{\pi}_{0}^{G}\left(\lambda,1\right)$ as
$\hat{\pi}_{0}^{G}\left(\lambda\right)$. Clearly
$\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)\leq\hat{\pi}_{0}^{S}\left(\lambda\right)$.
Since $\hat{\pi}_{0}^{S}\left(\lambda\right)\in\left[0,1\right]$ has to hold
even though (2) does not ensure this, we see either
$\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)=\hat{\pi}_{0}^{S}\left(\lambda\right)$
when each $F_{i}^{\ast}$, $i=1,...,m$ is $U\left(0,1\right)$ or
$\hat{\pi}_{0}^{G}\left(\lambda,0\right)=\hat{\pi}_{0}^{S}\left(\lambda\right)$.
Hence, $\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)$ generalizes
$\hat{\pi}_{0}^{S}\left(\lambda\right)$. Moreover, it distinguishes itself
from all existing estimators of $\pi_{0}$ by explicitly subtracting the
excessive upward bias caused by discrete $p$-values in its expectation, thus
giving more accurate and less conservative estimate of $\pi_{0}$.

###  3.1 Conservativeness of $\hat{\pi}_{0}^{G}$

We present the biases for $\hat{\pi}_{0}^{G}\left(\lambda,\varepsilon\right),$
$\hat{\pi}_{0}^{S}\left(\lambda\right)$ and $\tilde{\pi}_{0}^{P}$ respectively
and reveal their relationship.

###### Theorem 1

Regardless of the stochastic order of $F_{i}$ for $1\leq i\leq m$ relative to
$U\left(0,1\right)$, the following hold:

  1. 1.

The bias of $\hat{\pi}_{0}^{G}\left(\lambda,\varepsilon\right)$ is

$b^{G}\left(\varepsilon\right)=\dfrac{1-\varepsilon\lambda}{1-\lambda}-\dfrac{\sum\limits_{i=1}^{m}\left\\{F_{i}\left(\lambda\right)-\varepsilon
F_{i}^{\ast}\left(\lambda\right)\right\\}}{\left(1-\lambda\right)m}-\pi_{0}\text{.}$
(4)

  2. 2.

The bias of $\tilde{\pi}_{0}^{P}$ is

$b^{P}=\dfrac{2}{m}\sum_{i\in
I_{0}}\left(E\left[p_{i}\right]-0.5\right)+\dfrac{2}{m}\sum_{i\in
I_{1}}E\left[p_{i}\right]\text{.}$ (5)

  3. 3.

$\pi_{0}\leq
E\left[\hat{\pi}_{0}^{G}\left(\lambda\right)\right]\leq\min\left\\{E\left[\hat{\pi}_{0}^{S}\left(\lambda\right)\right],E\left[\tilde{\pi}_{0}^{P}\right]\right\\}$
if and only if

$0\leq b^{G}\left(1\right)\leq b^{P}\text{.}$ (6)

The proof of Theorem 1 follows from the definitions and is omitted. Though
Theorem 1 gives the necessary and sufficient condition to ensure that
$\hat{\pi}_{0}^{G}\left(\lambda\right)$ is the least conservative among the
three, both (4) and (6) are not available for $\lambda\neq 0$ due to the lack
of information on $F_{i}$, $i\in I_{1}$. Nevertheless, our simulation studies
in Section 5 show that $E\left[\hat{\pi}_{0}^{G}\left(0.5\right)\right]\leq
E\left[\tilde{\pi}_{0}^{P}\right]$ holds for various values of $\pi_{0}$.

###  3.2 Computing
$\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\varepsilon}\right)$ and choice of
tuning parameters

The definition of $\hat{\pi}_{0}^{G}$ implies the following algorithm for its
computation, which we state only for discrete $p$-values. Assume each
$F_{i}^{\ast}$ is discontinuous and let its support be

$S_{i}=\left\\{t\in\mathbb{R}:\int_{\left\\{t\right\\}}dF_{i}^{\ast}\left(t\right)>0\right\\}\text{.}$

Clearly, for each $\lambda$, there must be a unique $t_{i,\lambda}\in S_{i}$
such that

$t_{i,\lambda}=\sup\left\\{t\in S_{i}:\lambda\geq t\right\\}\text{.}$

Since $F_{i}^{\ast}\left(t\right)=t$ if $t\in S_{i}$ and that $F_{i}^{\ast}$
is càdlàg, it follows that

$\int_{(\lambda,1]}dF_{i}^{\ast}\left(t\right)=1-F_{i}^{\ast}\left(\lambda\right)=1-t_{i,\lambda}\geq
1-\lambda\text{.}$

Thus

$\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)=\dfrac{\sum\limits_{i=1}^{m}\left[\mathbf{1}_{\left\\{p_{i}>\lambda\right\\}}-\varepsilon_{i}\left(\lambda-
t_{i,\lambda}\right)\right]}{\left(1-\lambda\right)m}\text{,}$

where we have omitted the $\min$, $\max$ operators in (3) for notational
simplicity.

When $m$ is large but there is little or no information on $p$-value
distributions under the false null hypotheses, it is more feasible to set all
$\varepsilon_{i}=\varepsilon.$ So we focus on
$\hat{\pi}_{0}^{G}\left(\lambda,\varepsilon\right)$ and discuss choice of
tuning parameters $\lambda$ and $\varepsilon$. In the ideal case where
$b^{G}\left(\varepsilon\right)$ in (4) is non-negative for
$\left(\lambda,\varepsilon\right)$ in a set
$K_{0}\subseteq[0,1)\times\left[0,1\right]$ composed of finitely many
connected components, bootstrap method similar to that in Storey et al. (2004)
can be implemented to find the “best” tuning pair
$\left(\lambda_{\ast},\varepsilon_{\ast}\right)$. For completeness this method
is described below:

Step 1:

for a sufficiently small $d^{\ast}>0$, form a finite, compact set $K\subseteq
K_{0}$ such that its neighborhood of size $d^{\ast}$ covers $K_{0}$; compute
$\hat{\pi}_{0}^{G}\left(\lambda,\varepsilon\right)$ for each
$\left(\lambda,\varepsilon\right)\in K$.

Step 2:

for each $\left(\lambda,\varepsilon\right)\in K$, form $B$ bootstrap versions
$\hat{\pi}_{0,b}^{G}\left(\lambda,\varepsilon\right)$ of
$\hat{\pi}_{0}^{G}\left(\lambda,\varepsilon\right)$, $b=1,...,B$ by taking
bootstrap samples of the $p$-values.

Step 3:

for each $\left(\lambda,\varepsilon\right)\in K$, estimate its respective
mean-squared error (MSE) as

$\widehat{MSE}\left(\lambda,\varepsilon\right)=\dfrac{1}{B}\sum_{b=1}^{B}\left[\hat{\pi}_{0,b}^{G}\left(\lambda,\varepsilon\right)-\min_{\left(\tilde{\lambda},\tilde{\varepsilon}\right)\in
K}\left\\{\hat{\pi}_{0}^{G}\left(\tilde{\lambda},\tilde{\varepsilon}\right)\right\\}\right]^{2}\text{.}$

Step 4:

set
$\left(\lambda_{\ast},\varepsilon_{\ast}\right)=\arg\min_{\left(\lambda,\varepsilon\right)\in
K}\left\\{\widehat{MSE}\left(\lambda,\varepsilon\right)\right\\}$, and the
overall estimate of $\pi_{0}$ is
$\hat{\pi}_{0}^{G}\left(\lambda_{\ast},\varepsilon_{\ast}\right)$.

When $\varepsilon=0$ is set or the null $p$-values are distributed as
$U\left(0,1\right)$,
$\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)=\hat{\pi}_{0}^{S}\left(\lambda\right)$
and the above bootstrap method reduces to that in Storey et al. (2004) when
$K=\left\\{0,0.05,...,0.95\right\\}\times\left\\{0\right\\}$ is used as the
searching grid. However, knowing the impracticality of assuring
$b^{G}\left(\varepsilon\right)\geq 0$ in (4), we will not compute
$\hat{\pi}_{0}^{G}\left(\lambda_{\ast},\varepsilon_{\ast}\right)$. Instead, we
will empirically show in Section 5 that $\hat{\pi}_{0}^{G}\left(0.5\right)$
tends to be conservative and less upwardly biased for reasonable values of
$\pi_{0}$. We provide a heuristic explanation for the conservativeness of
$\hat{\pi}_{0}^{G}\left(0.5\right)$. Let $\pi=1-\pi_{0}$. Consider the choice
of $\lambda=1/2$ and $0\leq\varepsilon\leq 1$. Then
$E\left[\hat{\pi}_{0}^{G}\left(\lambda,\varepsilon\right)\right]=\pi_{0}+\beta_{0}+\beta_{1}$,
where $\beta_{0}=2m^{-1}\left(1-\varepsilon\right)\sum\limits_{i\in
I_{0}}\left\\{1/2-F_{i}\left(1/2\right)\right\\}\geq 0$ by definition and

$\displaystyle\beta_{1}$ $\displaystyle=$ $\displaystyle
2\left(1-\varepsilon/2\right)\pi-2m^{-1}\sum_{i\in
I_{1}}\left\\{F_{i}\left(1/2\right)-\varepsilon
F_{i}^{\ast}\left(1/2\right)\right\\}$ $\displaystyle>$ $\displaystyle
2\left(1-\varepsilon/2\right)\pi-2\pi=-\varepsilon\pi\text{.}$

When $\pi_{0}$ and $m$ are not too small and certain $F_{i}\left(1/2\right)$
are smaller than $1/2$, $\beta_{0}\geq\varepsilon\pi$ can hold for
$\varepsilon<1$; when $\varepsilon=1$, $\beta_{1}\geq 0$ can happen since
$F_{i}$, $i\in I_{1}$ do not have to all dominate $U\left(0,1\right)$.
Therefore, under such circumstance, choosing $\varepsilon$ close to or equal
to $1$ usually yields $\beta_{0}+\beta_{1}\geq 0$ and
$E\left[\hat{\pi}_{0}^{G}\left(\lambda,\varepsilon\right)\right]\geq\pi_{0}$.

##  4 FDR estimation using $\hat{\pi}_{0}^{G}$

With $\hat{\pi}_{0}^{G}$, we can develop better adaptive FDR procedures for
discrete $p$-values by extending existing ones. Storey’s FDR estimator in
Storey et al. (2004) is defined as

$\widehat{FDR}_{\lambda}\left(t\right)=\dfrac{\hat{\pi}_{0}^{S}\left(\lambda\right)t}{m^{-1}\left\\{R\left(t\right)\vee
1\right\\}}\text{,}$

whose variant is
$\widehat{FDR}_{\lambda}^{\ast}\left(t\right)=\dfrac{\tilde{\pi}_{0}^{S}\left(\lambda\right)t}{m^{-1}\left\\{R\left(t\right)\vee
1\right\\}}\mathbf{1}_{\left\\{t\leq\lambda\right\\}}+\mathbf{1}_{\left\\{t>\lambda\right\\}}$,
where
$\tilde{\pi}_{0}^{S}\left(\lambda\right)=\hat{\pi}_{0}^{S}\left(\lambda\right)+\left(1-\lambda\right)^{-1}m^{-1}$.
We remark that $\tilde{\pi}_{0}^{S}\left(\lambda\right)$ may have non-
negligible upward bias when $m$ is moderate because of its extra summand
$\left(1-\lambda\right)^{-1}m^{-1}$. Additionally, we point out that setting
$\widehat{FDR}_{\lambda}^{\ast}\left(t\right)=1$ for $t>\lambda$ in the
definition of $\widehat{FDR}_{\lambda}^{\ast}\left(t\right)$ is inappropriate
in certain situations. For example, we can easily construct $m>3$ $p$-values
such that $R\left(1/2\right)>2\left(m+1\right)/3$. For $1\geq t>1/2$ this
implies

$\dfrac{t\tilde{\pi}_{0}^{S}\left(1/2\right)}{m^{-1}\left\\{R\left(t\right)\vee
1\right\\}}=\dfrac{2t\left\\{m-R\left(1/2\right)+1\right\\}}{R\left(t\right)}<1\text{.}$

Consequently, setting $\widehat{FDR}_{1/2}^{\ast}\left(t\right)=1$ in this
case ignores the truth and may lead to erroneous conclusions.

In view of the above, we propose the generalized FDR estimator as

$\widetilde{FDR}_{\lambda}\left(t\right)=\min\left[1,\dfrac{\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)t}{m^{-1}\left\\{R\left(t\right)\vee
1\right\\}}\right]\text{.}$

Further, for a function $f$ with domain $\left[0,1\right]$, we define the
threshold of $f$ at level $\alpha$ as

$t_{\alpha}\left(f\right)=\sup\left\\{t\in\left[0,1\right]:f\left(t\right)\leq\alpha\right\\}\text{,}$

and let

$\mathcal{F}_{t}=\sigma\left(\mathbf{1}_{\left\\{p_{i}\leq s\right\\}},t\leq
s\leq 1,i=1,...,m\right)$

be the sigma-algebra generated by the indicator functions inside the
parentheses.

###  4.1 Stopping time property of the threshold of
$\widetilde{FDR}_{\lambda}\left(t\right)$

Recently, martingale methods have been applied to obtain the conservativeness
of FDR procedures using the stopping time property of relevant thresholds with
respect to certain filtrations. For example, the stopping time property of
$t_{\alpha}\left(\widehat{FDR}_{\lambda}^{\ast}\right)$ is claimed in Lemma 4
of Storey et al. (2004) without a formal proof, so is that of
$\alpha^{\ast}\left(q^{\ast}\right)$ in (7.2) of Pena et al. (2011). These
claims have not been proved in the setting of independent, uniform null
$p$-values, and it is not clear whether this property holds for the
generalized FDR estimator and for general $p$-values. We show that
$t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)$ is a stopping time with
respect to the backward filtration

$\mathcal{G}=\left\\{\mathcal{F}_{t\wedge\lambda}:1\geq t\geq 0\right\\}$

for general $p$-values, discrete or continuous, independent or dependent,
where $a_{0}\wedge a_{1}=\min\left\\{a_{0},a_{1}\right\\}$. Further, we show
that $\widetilde{FDR}_{\lambda}$ reaches level $\alpha$ at
$t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)$ under one simple condition.

To achieve this, we start with $L\left(t\right)=t\left\\{R\left(t\right)\vee
1\right\\}^{-1}$, $t\in\left[0,1\right]$, the scaled inverse rejection
process. Order the $p$-values into
$p_{\left(1\right)}<p_{\left(2\right)}<...<p_{\left(n\right)}$ distinctly,
where the multiplicity of $p_{\left(i\right)}$ is $n_{i}$ for $i=1,...,n$. Let
$p_{\left(n+1\right)}=\max\left\\{p_{\left(n\right)},1\right\\}$ and
$p_{\left(0\right)}=0.$ Define $T_{j}=\sum\limits_{l=1}^{j}n_{l}$ for
$j=1,...,n$.

###### Lemma 1

The process $\left\\{L\left(t\right),t\in[0,1]\right\\}$ is such that

$L\left(t\right)=\left\\{\begin{array}[]{lcl}t&\text{if}&t\in[0,p_{\left(1\right)})\text{,}\\\
tT_{j}^{-1}&\text{if}&t\in[p_{\left(j\right)},p_{\left(j+1\right)})\text{ for
}j=1,..,n-1\text{,}\\\
tm^{-1}&\text{if}&t\in[p_{\left(n\right)},p_{\left(n+1\right)}]\text{.}\end{array}\right.$
(7)

Moreover, it can only be discontinuous at $p_{\left(i\right)}$, $1\leq i\leq
n$, where it can only have a downward jump with size

$L\left(p_{\left(i\right)}-\right)-L\left(p_{\left(i\right)}\right)=\dfrac{p_{\left(i\right)}n_{i}}{R\left(p_{\left(i\right)}\right)\left\\{R\left(p_{\left(i\right)}\right)-n_{i}\right\\}}>0\text{.}$

The conclusion of Lemma 1 is right the contrary to the claim in the proof of
Theorem 2 in Storey et al. (2004) that “the process $mt/R\left(t\right)$ has
only upward jumps and has a final value of $1$”, since it says “the process
$mt/R\left(t\right)$ has only downward jumps”. We construct a counterexample
to their claim as follows. For a small increase $c$ in $t$ which results
increase $a_{c}$ in $R\left(t\right)$, we see that

$L\left(t+c\right)-L\left(t\right)=\dfrac{t+c}{R\left(t\right)+a_{c}}-\dfrac{t}{R\left(t\right)}=\dfrac{cR\left(t\right)-ta_{c}}{\left\\{R\left(t\right)+a_{c}\right\\}R\left(t\right)}<0$

if and only if $\dfrac{c}{a_{c}}<\dfrac{t}{R\left(t\right)}$. Construct $m$
$p$-values with $n\geq 4$ such that there exists some $1\leq j_{0}<n-2$ with
$n_{j_{0}+1}>T_{j_{0}}$ but $p_{\left(j_{0}+1\right)}<1$. Choose $c_{1}$ and
$c_{2}$ such that
$0<c_{1}<\left\\{p_{\left(j_{0}+1\right)}-p_{\left(j_{0}\right)}\right\\}/2$
and $0<c_{2}<p_{\left(j_{0}+1\right)}-2c_{1}$. Let
$t_{0}=p_{\left(j_{0}+1\right)}-c_{1}$ and $c=c_{1}+c_{2}$. Then
$p_{\left(j_{0}\right)}<t_{0}<p_{\left(j_{0}+1\right)}$,
$R\left(t_{0}\right)=T_{j_{0}}$ and
$R\left(t_{0}+c\right)=T_{j_{0}}+n_{j_{0}+1}$. Further, $0<c<t_{0}$ and
$a_{c}=n_{j_{0}+1}$. So $\dfrac{c}{n_{j_{0}+1}}<\dfrac{t_{0}}{T_{j_{0}}}$ and
$L\left(t_{0}+c\right)-L\left(t_{0}\right)<0$. Letting $c\rightarrow 0$ gives
$p_{\left(j_{0}+1\right)}$ as a point of downward jump for $L\left(t\right)$.

Equipped with the properties of $L\left(\cdot\right)$, we have

###### Theorem 2

$t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)$ is a stopping time for
$\mathcal{G}$. Further,

$\widetilde{FDR}_{\lambda}\left\\{t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\right\\}=\alpha$

whenever $\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)>\alpha$ and
$R\left(t\right)>0$.

The conclusion of Theorem 2 is illustrated by an application of
$\hat{\pi}_{0}^{G}$ and $\widetilde{FDR}_{\lambda}$ in Section 5. Define

$\hat{\pi}_{0}\left(\lambda,\varsigma\right)=\max\left[0,\min\left\\{1,\hat{\pi}_{0}^{S}\left(\lambda\right)-\dfrac{\varsigma}{\left(1-\lambda\right)m}\right\\}\right]$

and

$\widetilde{FDR}_{\lambda}^{\varsigma}\left(t\right)=\min\left[1,\dfrac{\hat{\pi}_{0}\left(\lambda,\varsigma\right)t}{m^{-1}\left\\{R\left(t\right)\vee
1\right\\}}\right]\text{,}$

where
$\varsigma\in\left[0,\left(1-\lambda\right)m\hat{\pi}_{0}^{S}\left(\lambda\right)\right]$
is deterministic and does not depend on $t$. We call
$\widetilde{FDR}_{\lambda}^{\varsigma}$ a “Storey-type FDR estimator”. Using
the same arguments in the proof of Theorem 2, we can show

###### Theorem 3

$t_{\alpha}\left(\widetilde{FDR}_{\lambda}^{\varsigma}\right)$ is a stopping
time for $\mathcal{G}$. Further,

$\widetilde{FDR}_{\lambda}^{\varsigma}\left\\{t_{\alpha}\left(\widetilde{FDR}_{\lambda}^{\varsigma}\right)\right\\}=\alpha$

whenever $\hat{\pi}_{0}\left(\lambda,\varsigma\right)>\alpha$ and
$R\left(t\right)>0$.

The proof of Theorem 3 is omitted. Theorem 3 implies that, regardless of
whether the $p$-values are independent or continuous, the stopping time
property of the thresholds of Storey-type FDR estimators is generic, and that
the FDR estimator reaches the prespecified FDR level at its threshold whenever
the estimated proportion of true nulls is larger than this FDR level and there
is at least one rejection. It provides the first formal and most general
support for the use of the stopping time property of the thresholds of a
certain family of adaptive FDR estimators.

###  4.2 Large sample properties of $\widetilde{FDR}_{\lambda}$

Even though Theorem 3 reveals that $\widetilde{FDR}_{\lambda}^{\varsigma}$
reaches the prespecified FDR level at its threshold, it does not tell whether
$\widetilde{FDR}_{\lambda}^{\varsigma}$ conservatively estimates the FDR of a
one-step MTP whose rejection threshold is this threshold. Surprisingly, under
appropriate conditions, a much stronger property of
$\widetilde{FDR}_{\lambda}^{\varsigma}$ can be obtained that extends similar
results in Storey et al. (2004). We illustrate this using
$\widetilde{FDR}_{\lambda}$ and assume that

A1)

$m_{0}^{-1}V\left(t\right)$ and $m_{1}^{-1}S\left(t\right)$ converge almost
surely to càdlàg functions $G_{0}\left(t\right)$ and $G_{1}\left(t\right)$
such that $0<G_{0}\left(t\right)\leq t$ for all $t\in(0,1]$ and
$\lim\limits_{m\rightarrow\infty}m_{0}/m=\pi_{0,\infty}\in\left(0,1\right)$
exists.

When $G_{i}\left(\cdot\right)$, $i=0,1$ in A1) are continuous, Storey et al.
(2004) proved that $\widehat{FDR}_{\lambda}$ is simultaneously asymptotically
conservative in the sense that, for any $0<\delta\leq 1$, $P$-a.s.

$\left\\{\begin{array}[]{c}\lim\limits_{m\rightarrow\infty}\inf_{t\geq\delta}\left\\{\widehat{FDR}_{\lambda}\left(t\right)-FDR\left(t\right)\right\\}\geq
0\text{,}\\\
\lim\limits_{m\rightarrow\infty}\inf_{t\geq\delta}\left\\{\widehat{FDR}_{\lambda}\left(t\right)-\dfrac{V\left(t\right)}{R\left(t\right)\vee
1}\right\\}\geq 0\text{.}\end{array}\right.$ (8)

We show that these properties are endowed to $\widetilde{FDR}_{\lambda}$.
Noticing that (8) hinges on the total conservativeness of
$\hat{\pi}_{0}^{S}\left(\lambda\right)$ in that
$\hat{\pi}_{0}^{S}\left(\lambda\right)\geq\pi_{0,\infty}$, $P$-a.s.
automatically once A1) holds, we assume

A2)

$\lim\limits_{m\rightarrow\infty}\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)\geq\pi_{0,\infty}$,
$P$-a.s..

Thus we have

###### Theorem 4

Suppose A1) and A2) hold. Then for each $\delta\in(0,1]$,

$\liminf_{m\rightarrow\infty}\inf_{1\geq
t\geq\delta}\left\\{\widetilde{FDR}_{\lambda}\left(t\right)-\dfrac{V\left(t\right)}{R\left(t\right)\vee
1}\right\\}\geq 0$ (9)

and

$\liminf_{m\rightarrow\infty}\inf_{1\geq
t\geq\delta}\left\\{\widetilde{FDR}_{\lambda}\left(t\right)-FDR\left(t\right)\right\\}\geq
0$ (10)

hold $P$-a.s..

Theorem 4 ensures the simultaneous conservativeness of
$\widetilde{FDR}_{\lambda}$ when $m$ is sufficiently large. On the other hand,
we have

$\widetilde{FDR}_{\lambda}\left(t\right)-\widehat{FDR}_{\lambda}\left(t\right)=\dfrac{\delta\left(\lambda,\mathbf{\epsilon}\right)t}{m^{-1}\left\\{R\left(t\right)\vee
1\right\\}}\leq 0$ (11)

and $\sup_{1\geq t\geq
0}\dfrac{\delta\left(\lambda,\mathbf{\epsilon}\right)t}{m^{-1}\left\\{R\left(t\right)\vee
1\right\\}}\geq\delta\left(\lambda,\mathbf{\epsilon}\right)$. Hence,
$\widetilde{FDR}_{\lambda}$ may fail to conservatively estimate the FDR of a
one-step MTP when the needed conditions are violated. However, since (11)
implies $t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\geq
t_{\alpha}\left(\widehat{FDR}_{\lambda}\right)$, it follows that, at the same
FDR level, a one-step MTP using the threshold of $\widetilde{FDR}_{\lambda}$
is usually more powerful than one that uses the threshold of
$\widehat{FDR}_{\lambda}$.

When more asymptotic uniformity of the $p$-value distributions is available,
such as

A3)

$\lim\limits_{m\rightarrow\infty}m^{-1}\sum\limits_{i=1}^{m}\varepsilon_{i}=\bar{\varepsilon}$
and $\lim\limits_{m\rightarrow\infty}m_{1}^{-1}\sum\limits_{i\in
I_{1}}F_{i}^{\ast}\left(\lambda\right)=h_{1}\left(\lambda\right)$,

then with A1), we have

$\lim_{m\rightarrow\infty}m^{-1}\sum_{i=1}^{m}\varepsilon_{i}\left\\{\lambda-
F_{i}^{\ast}\left(\lambda\right)\right\\}=h^{\ast}\left(\lambda\right)\text{,}$

where
$h^{\ast}\left(\lambda\right)=t\bar{\varepsilon}-\pi_{0,\infty}G_{0}\left(\lambda\right)-\pi_{1,\infty}h_{1}\left(\lambda\right)$
and $\pi_{1,\infty}=1-\pi_{0,\infty}$. Further, the limit of
$\widetilde{FDR}_{\lambda}\left(t\right)$ as $m\rightarrow\infty$ is

$\widetilde{FDR}_{\lambda}^{\infty}\left(t\right)=\dfrac{\left[\left\\{1-G_{0}\left(\lambda\right)\right\\}\pi_{0,\infty}+\left\\{1-G_{1}\left(\lambda\right)\right\\}\pi_{1,\infty}\right]t-h^{\ast}\left(\lambda\right)t}{\left(1-\lambda\right)\left\\{\pi_{0}G_{0}\left(t\right)+\pi_{1}G_{1}\left(t\right)\right\\}}\text{.}$

With these additional assumptions, we can show that a one-step MTP using
$t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)$ as the rejection threshold
maintains FDR control at level $\alpha$ asymptotically.

###### Theorem 5

Suppose A1) to A3) hold. If
$\widetilde{FDR}_{\lambda}^{\infty}\left(t\right)<\alpha$ for some
$t\in\left[0,1\right]$, then

$\limsup_{m\rightarrow\infty}FDR\left\\{t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\right\\}\leq\alpha\text{.}$

We remark that the conclusions in Theorem 4 and Theorem 5 hold for Storey-type
FDR estimators with $\hat{\pi}_{0}^{G}$ replaced by any totally conservative
estimator $\hat{\pi}_{0}$ of $\pi_{0}$ when the involved proportions and
empirical processes converge appropriately. When $\hat{\pi}_{0}^{G}$ is not
conservative or totally conservative, $\widetilde{FDR}_{\lambda}$ may not
conservatively estimate the FDR of a one-step MTP. Fortunately, we can
quantify the possible downward bias of $\widetilde{FDR}_{\lambda}$. Set
$T_{n+1}=T_{n}$. For $\delta\in(0,1]$, define

$L_{\inf}^{\infty}=\liminf_{m\rightarrow\infty}\min_{1\leq j\leq
n+1}\left\\{\dfrac{p_{\left(j\right)}}{m^{-1}T_{j}}:p_{\left(j\right)}\geq\delta\right\\}\text{,}$

and

$C_{\delta,\lambda}\left(t\right)=\liminf_{m\rightarrow\infty}\inf_{1\geq
t\geq\delta}\left\\{\widehat{FDR}_{\lambda}\left(t\right)-FDR\left(t\right)\right\\}\text{.}$

###### Theorem 6

Suppose A1) and A3) hold. Then for any $\delta\in(0,1]$,

$C_{\delta,\lambda}\left(t\right)-\dfrac{h^{\ast}\left(\lambda\right)}{1-\lambda}\leq\lim_{m\rightarrow\infty}\inf_{1\geq
t\geq\delta}\left\\{\widetilde{FDR}_{\lambda}\left(t\right)-FDR\left(t\right)\right\\}\leq
C_{\delta,\lambda}\left(t\right)-\dfrac{h^{\ast}\left(\lambda\right)L_{\inf}^{\infty}}{1-\lambda}$

with probability $1$.

With some algebra, the proof of Theorem 6 follows from Lemma 1, (10) and (11),
and is thus omitted.

##  5 Simulation study and an application

Now we assess the performance of $\hat{\pi}_{0}^{G}$ and
$\widetilde{FDR}_{\lambda}$ via simulation studies. Since $\hat{\pi}_{0}^{G}$
reduces to $\hat{\pi}_{0}^{S}$ for continuous $p$-values whose null
distribution is $U\left(0,1\right)$ and the performance of $\hat{\pi}_{0}^{S}$
has been well documented, we compare the competing estimators of $\pi_{0}$
based on discrete, two-sided $p$-values induced by BINs, FETs and ENTs. We do
not assess $\hat{\pi}_{0}^{G}$ or $\widetilde{FDR}_{\lambda}$ for dependent
discrete $p$-values, since the deterministic term
$\delta\left(\lambda,\mathbf{\epsilon}\right)$ as the difference between
$\hat{\pi}_{0}^{S}$ and $\hat{\pi}_{0}^{G}$ and the extensive study of
$\hat{\pi}_{0}^{S}$ under dependence (Storey et al.; 2004; Blanchard and
Roquain; 2009) have already characterized the performance of
$\hat{\pi}_{0}^{G}$ and $\widetilde{FDR}_{\lambda}$ in this situation.
Finally, even though we have justified the simultaneous asymptotic
conservativeness of $\widetilde{FDR}_{\lambda}$, we do not study
$\widetilde{FDR}_{\lambda}$ using dependent discrete $p$-values due to the
possible, practical divergence of the required empirical processes related to
heterogeneous, discrete $p$-value distributions.

Specifically, we will compare $\hat{\pi}_{0}^{G}\left(0.5\right)$ against
$\hat{\pi}_{0}^{S}\left(0.5\right)$, $\tilde{\pi}_{0}^{P}$ and

$\hat{\pi}_{0}^{B}=m^{-1}\left(m-\left[m/2\right]+1\right)\left(1-p_{\left(\left[m/2\right]\right)}\right)^{-1}$

of Benjamini et al. (2006), and compare $\widetilde{FDR}_{\lambda}$ against
$\widehat{FDR}_{\lambda}$ and the adaptive BH procedure $FDR_{BH}^{\ast}$ with
$\alpha$ in the BH procedure being replaced by $\alpha/\hat{\pi}_{0}^{B}$. We
choose $\hat{\pi}_{0}^{S}\left(0.5\right)$ since discrete $p$-values can take
value $1$ and the qvalue package of Storey (2002) reports error when
estimating $\pi_{0}$ from such $p$-values, whereas the choice of
$\hat{\pi}_{0}^{B}$ using the median is to keep the robustness
$\hat{\pi}_{0}^{B}$ across different settings of discrete $p$-values since
$\hat{\pi}_{0}^{B}$ does not adjust for the discreteness of such $p$-values.

The simulation is set up as follows:

  1. 1.

Set $m=1000$, $\pi_{0}\in\left\\{0.5,0.6,0.7,0.8,0.95\right\\}$ and
$\alpha\in\left\\{0.025,0.05,0.075,0.1\right\\}$.

  2. 2.

For each fixed pair $\left(\pi_{0},\alpha\right)$,

    1. (a)

Simulation I for Poisson data: generate $m$ $\theta_{i1}$’s independently from
a Pareto distribution Pareto$\left(l,s\right)$ with location $l=7$ and shape
$s=7$; generate $m_{1}$ $\rho_{i}$’s independently from $U\left(1.5,5\right)$;
set $\theta_{i2}=\theta_{i1}$ for $1\leq i\leq m_{0}$ but
$\theta_{i2}=\rho_{i}\theta_{i1}$ for $m_{0}+1\leq i\leq m$; generate $2m$
$\xi_{ig}$’s with $g=1,2$ independently from Poisson distributions
Poisson$\left(\theta_{ig}\right)$ with mean $\theta_{ig}$.

    2. (b)

Simulation II for Binomial data: generate $m$ $\tilde{r}_{i}$’s independently
from Negative Binomial distribution
$\mathsf{NegBinomial}\left(\sigma,\mu\right)$ with size $\sigma=3$ and mean
$\mu=8$ and set $r_{i}=\tilde{r}_{i}+2$; generate $m$ $\theta_{i1}$’s
independently from $U\left(0.08,0.65\right)$ and $m_{1}$ $\rho_{i}$’s
independently from $U\left(1.5,13\right)$; set set $\theta_{i2}=\theta_{i1}$
for $1\leq i\leq m_{0}$ but
$\theta_{i2}=\rho_{i}\dfrac{\theta_{i1}}{1-\theta_{i1}}$ for $m_{1}+1\leq
i\leq m$; generate $2m$ $\xi_{ig}$’s with $g=1,2$ independently from Binomial
distributions Binomial$\left(\theta_{ig},r_{i}\right)$ with probability of
success $\theta_{ig}$ and number of trials $r_{i}$.

    3. (c)

Simulation III for Negative Binomial data: from the data arab from the R
package NBPSeq developed by Di et al. (2011), take all rows each of whose row
sums is between $1$ and $30$; process these rows by R package edgeR developed
by Robinson et al. (2010) to estimate the common dispersion $\phi=\sigma^{-1}$
being $1.451$, the probabilities $\tau_{i1}$ and the means $\theta_{i1}$ for
the counts in the control group modelled as
$\mathsf{NegBinomial}\left(\sigma,\theta_{i1}\right)$ with
$\theta_{i1}=\dfrac{\tau_{i1}\sigma}{1-\tau_{i1}}$; generate independently
$m_{1}$ $\rho_{i}$’s from Pareto$\left(1.5,1.426\right)$; set
$\theta_{i2}=\theta_{i1}$ for $1\leq i\leq m_{0}$ but
$\theta_{i2}=\rho_{i}\theta_{i1}$ for $m_{0}+1\leq i\leq m$; for each $i$
between $1$ and $m$ and $g$ in $\left\\{1,2\right\\}$, generate $3$
$\xi_{ig}$’s independently from
$\mathsf{NegBinomial}\left(\sigma,\theta_{ig}\right)$.

  3. 3.

With $\xi_{ig}$, $g=1,2$ for each $i$, conduct the BIN for Poisson data, FET
for Binomial data, and ENT for Negative Binomial data to test

$H_{i0}:\theta_{i2}=\theta_{i1}\text{ versus
}H_{i1}:\theta_{i2}\neq\theta_{i1},$

where $H_{i0}$ denotes a true null and $H_{i1}$ its corresponding false null;
obtain the two-sided $p$-value $p_{i}$ of the test.

  4. 4.

Apply $\hat{\pi}_{0}^{S}\left(0.5\right)$,
$\hat{\pi}_{0}^{G}\left(0.5\right)$, $\tilde{\pi}_{0}^{P}$,
$\hat{\pi}_{0}^{B}$, and the competing MTPs to the $m$ $p$-values
$\left\\{p_{i}\right\\}_{i=1}^{m}$.

  5. 5.

Repeat Steps (b) to (d) $250$ times to obtain certain statistics for each
$\hat{\pi}_{0}$ and for the MTPs.

###  5.1 Summary from simulation studies

Figure 1 displays, for different values of $\pi_{0}$, the excess
$\hat{\pi}_{0}-\pi_{0}$ of a $\hat{\pi}_{0}$ when it is applied to the
$p$-values induced by three types of tests, namely, BIN, FET and ENT. It can
be seen that $\hat{\pi}_{0}^{G}$ is stable and the most accurate when applied
to the $p$-values induced by BINs and FETs. For these two cases, even though
$\hat{\pi}_{0}^{G}$ may under-estimate $\pi_{0}$, it is usually conservative
in that $E\left[\hat{\pi}_{0}^{G}\right]\geq\pi_{0}$. All estimators of
$\pi_{0}$, especially $\hat{\pi}_{0}^{G}$, display aggressiveness, i.e.,
under-estimation, when applied to the $p$-values of the ENTs when $\pi_{0}$ is
large. This is due to the erroneous distortion imposed by edgeR to the
distributions of the test statistics of the ENTs and hence to those of the
$p$-values, since edgeR severely under-estimates as $\hat{\phi}$ the
dispersion $\phi=\sigma^{-1}$ from the original simulated counts and
normalizes these counts using $\hat{\phi}$; see Figure 2 for plots of
$\hat{\phi}$. In view of this, we omitted the case of MT with $p$-values of
the ENTs when presenting in Figure 3, Figure 4 and Figure 5 the number of
rejections, rejection thresholds, and false discovery proportions (FDPs,
Genovese and Wasserman; 2002),
$FDP\left(t\right)=V\left(t\right)\left\\{R\left(t\right)\right\\}^{-1}\mathbf{1}_{\left\\{R\left(t\right)>0\right\\}}+\mathbf{1}_{\left\\{R\left(t\right)=0\right\\}}$,
of the competing MTPs. These figures show, for each nominal FDR level
$\alpha$, that $\widetilde{FDR}_{0.5}$ usually has a larger rejection
threshold and more rejections, that its FDR is no larger than $\alpha$, and
that $\widetilde{FDR}_{0.5}$ usually has competitive FDPs.

Our simulation studies demonstrate that, when the discrete $p$-values have
null distributions dominating the uniform distribution and that they
reasonably reflect the power of the discrete test statistics that induce them,
the generalized estimators are preferred and can achieve more discoveries
subject to an FDR level.

###  5.2 Application to Arabidopsis thaliana data

We now apply the generalized estimators to an NGS data set. Lister et al.
(2008) studied cytosine methylation of two unreplicated lines of Arabidopsis
thaliana, namely, wild-type (Col-0) and mutant defective (Met1-3). It is
biologically known that the proportion of methylated cytosines are relatively
large in such a study. To assess if a cytosine is differential methylated
(DM), Lister et al. (2008) used FET on the discrete counts indicating the
level of methylation for each cytosine in Col-0 and Met1-3. Among the $22265$
cytosines, we extracted those whose total counts for both lines are positive
but no larger than $25$. This yields $3945$ cytosines. We model the counts for
each cytosine by independent Binomial distributions and use FET to test for
each cytosine the null: the two Binomial distributions have the same success
probability (namely, a cytosine is non-DM). At an FDR level of $\alpha=0.05$,
$\widehat{FDR}_{\lambda}$, $\widetilde{FDR}_{\lambda}$, and the BH procedure
are applied to the two-sided $p$-values of the $3945$ tests. The results are
summarized in Table 1.

Compared to Storey’s procedure and the BH procedure, respectively $72$ and
$146$ more DM cytosines were found using the generalized estimators. The
generalized FDR estimator gives the largest rejection threshold
$t_{0.05}\left(\widetilde{FDR}_{1}\right)=0.009550824$ among the three
competing FDR procedures. This illustrates the improvement $\hat{\pi}_{0}^{G}$
and $\widetilde{FDR}_{\lambda}$ can induce when applied to discrete $p$-values
whose distributions may have a few large jumps and whose null distributions
dominate $U(0,1)$. We remark that $\widetilde{FDR}_{0.05}=0.049999678$ at its
threshold is the consequence of the step size used in numeric search and it
does not affect the validity of Theorem 2.

Table 1: Results from application of the generalized estimators to Arabidopsis
thaliana data. “New” means using
$\hat{\pi}_{0}^{G}\left(\lambda,\varepsilon\right)$ and
$\widetilde{FDR}_{\lambda}$, “Storey” means using
$\hat{\pi}_{0}^{S}\left(\lambda\right)=\hat{\pi}_{0}^{G}\left(\lambda,0\right)$
and $\widehat{FDR}_{\lambda}$, and “BH” is the BH procedure; $\hat{t}$ is the
threshold of a one-step MTP, and $\hat{\alpha}\left(\hat{t}\right)$ the
estimated FDR of such an MTP that uses rejection threshold $\hat{t}$. “NA”
means “not applicable”. It can be seen that $\widetilde{FDR}_{0.5}$ makes more
discoveries subject to FDR being $0.05$.

Method $(\lambda,\varepsilon)$ $\widehat{\pi}_{0}$ $\widehat{t}$
$\hat{\alpha}\left(\hat{t}\right)$ $R(\widehat{t})$ New (0.5,1) 0.5984886
0.009550824 0.049999678 451 Storey (0.5,0) 0.7609632 0.006312294 0.049998698
379 BH NA 1 0.003865653 NA 305

##  6 Discussion

We have proposed generalized estimators of the proportion of true nulls and of
the FDR of one-step MTPs. The new estimators usually outperform existing ones
for discrete $p$-values whose null distributions dominate the uniform, and
they reduce to Storey’s estimators for continuous $p$-values. We have
established the asymptotic conservativeness of the generalized FDR estimator
and provided the first formal and most general justification of the stopping
time property of the threshold of Storey-type FDR estimators for $p$-values
not necessarily continuous or independent. Our strategy to adjust estimators
of $\pi_{0}$ can be easily adapted to other types of discrete p-values.

It should be noted that $\hat{\pi}_{0}^{G}$ can be non-conservative when its
tuning parameters are improperly chosen and that in this case
$\widetilde{FDR}_{\lambda}$ may fail to conservatively estimate the FDR of a
one-step MTP. However, due to the lack of information on the $p$-value
distribution under the false null (i.e., alternative) hypotheses, how to
choose these tuning parameters to ensure the conservativeness of
$\hat{\pi}_{0}^{G}$ seems to be a very tough issue. We leave the investigation
on this to future research but point out that the techniques in Liang and
Nettleton (2012) are not applicable to heterogenous null $p$-value
distributions treated here. We remark that when $\pi_{0}$ is close to $1$
and/or each null discrete $p$-value distribution is very close to the uniform
distribution, $\hat{\pi}_{0}^{G}$ and $\widetilde{FDR}_{\lambda}$ can not
result in moderate improvements. However, the advantage of the generalized
estimators will be well manifested when $\pi_{0}$ is not close to $1$ and the
c.d.f.’s of the involved discrete test statistics have small supports with a
few large jumps.

Driven by the accumulation of massive discrete data together with accompanying
statistical tests to be conducted, the need for better estimators of the
proportion of true nulls and of the FDR of an MTP for this type of data is
urgent. We expect our methodology to be useful in multiple testing with
discrete data. In this work, we have mainly tackled the issue of multiple
testing with discrete $p$-values from the perspective of obtaining a better
estimator of the proportion of true null hypotheses, we shall report in a
future article our novel strategy to deal directly with multiple testing under
heterogeneous null $p$-value distributions.

## Acknowledgements

This research was funded in part by a National Science Foundation Plant Genome
Research Program grant (No. IOS-1025976) to R.W. Doerge, and it produced as a
technical report Chen and Doerge (2012) a version of this work with different
simulation studies and application. X. Chen would like to thank Douglas
Baumann and Sanvesh Srivastava for helpful discussions on NGS data and he is
especially grateful to Douglas Baumann for formatting the counts data in
Lister et al. (2008) to which the generalized estimators were applied.

##  Appendix A Proofs

###  A.1 Proof of Lemma 1

Clearly

$R\left(t\right)=\left\\{\begin{array}[]{ccl}0&\text{if}&0\leq
t<p_{\left(1\right)}\text{,}\\\ T_{j}&\text{if}&p_{\left(j\right)}\leq
t<p_{\left(j+1\right)},j=1,..,n-1\text{,}\\\
m&\text{if}&p_{\left(n\right)}\leq t\leq
p_{\left(n+1\right)}\text{,}\end{array}\right.$

and (7) holds. Therefore, the points of discontinuities of
$L\left(\cdot\right)$ are the original distinct $p$-values. This justifies the
first part of the assertion.

Now we show that $L\left(\cdot\right)$ can only have downward jumps at points
of discontinuity. Define

$\varphi\left(t,\eta\right)=\dfrac{t+\eta}{R\left(t+\eta\right)}-\dfrac{t}{R\left(t\right)}=\dfrac{\eta
R\left(t\right)+t\left\\{R\left(t\right)-R\left(t+\eta\right)\right\\}}{R\left(t+\eta\right)R\left(t\right)}\text{.}$

From the fact

$R\left(p_{\left(j\right)}\right)-R\left(p_{\left(j\right)}-\right)=n_{j}>0$

but $R\left(p_{\left(j\right)}+\right)-R\left(p_{\left(j\right)}\right)=0$ for
each $1\leq j\leq n,$ it follows that
$\varphi\left(p_{\left(j\right)},0+\right)=\lim_{\eta\downarrow
0}\varphi\left(t,\eta\right)=0$ but

$\varphi\left(p_{\left(j\right)},0-\right)=\lim_{\eta\uparrow
0}\varphi\left(p_{\left(j\right)},\eta\right)=\dfrac{p_{\left(j\right)}n_{j}}{R\left(p_{\left(j\right)}\right)\left\\{R\left(p_{\left(j\right)}\right)-n_{j}\right\\}}>0\text{.}$

Thus
$L\left(p_{\left(j\right)}-\right)-L\left(p_{\left(j\right)}\right)=\varphi\left(p_{\left(j\right)},0-\right)>0$
and the proof is completed.

###  A.2 Proof of Theorem 2

Define

$\mathbf{X}_{t}^{\left(m\right)}\left(\omega\right)=\left(1_{\left\\{p_{1}\leq
t\right\\}}\left(\omega\right),...,1_{\left\\{p_{m}\leq
t\right\\}}\left(\omega\right)\right)\text{, }\omega\in\Omega\text{.}$

Then
$\mathcal{F}_{t}=\sigma\left(\mathbf{X}_{s}^{\left(m\right)}\left(\omega\right),1\geq
s\geq t\right)$, $t\in\left[0,1\right]$ and $\left\\{\mathcal{F}_{t}:0\leq
t\leq 1\right\\}$ is a decreasing sequence of sub-sigma-algebras of
$\mathcal{A}$. Write
$\widetilde{FDR}_{\lambda}\left(t\right)=\widetilde{FDR}_{\lambda}\left(t,\omega\right)$
and $t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)=\tilde{t}_{\alpha}$.
When $\tilde{t}_{\alpha}=1$ or $0$ corresponding to $\alpha=1$ or $\alpha=0$,
it is automatically a stopping time. So we only consider
$\alpha\in\left(0,1\right)$ and $0<\tilde{t}_{\alpha}<1$. By definition,

$\left\\{\omega\in\Omega:\tilde{t}_{\alpha}\leq
s\right\\}=\bigcap\limits_{\left\\{t:s<t\leq
1\right\\}}A_{t}=\tilde{A}_{s}\text{,}$

where
$A_{t}=\left\\{\omega\in\Omega:\widetilde{FDR}_{\lambda}\left(t,\omega\right)>\alpha\right\\}$,
$1\geq t>s$. We need to show $\tilde{A}_{s}\in\mathcal{F}_{s\wedge\lambda}$.

Since either
$\mathcal{F}_{s\wedge\lambda}=\mathcal{F}_{s}\supseteq\mathcal{F}_{\lambda}$
when $s\leq\lambda$ or
$\mathcal{F}_{s\wedge\lambda}=\mathcal{F}_{\lambda}\supseteq\mathcal{F}_{s}$
when $s\geq\lambda$, the stopping time property holds once we prove
$\tilde{A}_{s}\in\mathcal{F}_{s}$. With the decomposition

$A_{t}=\bigcup\nolimits_{r\in\mathbb{Q}}\left(A_{t,r}\cap B_{r}\right)$

and

$\tilde{A}_{s}=\bigcap\limits_{\left\\{t:s<t\leq
1\right\\}}\bigcup\limits_{r\in\mathbb{Q}}\left(A_{t,r}\cap
B_{r}\right)=\bigcup\limits_{r\in\mathbb{Q}}\bigcap\limits_{\left\\{t:s<t\leq
1\right\\}}\left(A_{t,r}\cap B_{r}\right)\text{,}$

where
$A_{t,r}=\left\\{\omega\in\Omega:\dfrac{t}{m^{-1}\left\\{R\left(t\right)\vee
1\right\\}}\geq r\right\\}$ and
$B_{r}=\left\\{\omega\in\Omega:\dfrac{\alpha}{\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)}<r\right\\}$,
it suffices to show

$\bigcap\limits_{\left\\{t:s<t\leq 1\right\\}}\left(A_{t,r}\cap
B_{r}\right)=\bigcap\limits_{\left\\{t:s<t\leq 1\right\\}}A_{t,r}\cap
B_{r}\in\mathcal{F}_{s}\text{,}$

or equivalently,

$A_{s,r}^{\ast}=\bigcap\limits_{\left\\{t:s<t\leq
1\right\\}}A_{t,r}\in\mathcal{F}_{s}\text{,}$

since $B_{r}\in\mathcal{F}_{\lambda}$ holds already.

Define $I_{i}=[p_{\left(i\right)},p_{\left(i+1\right)})$, $i=1,...,n-1$. We
will add $I_{0}=[p_{\left(0\right)},p_{\left(1\right)})$ and
$I_{n}=[p_{\left(n\right)},p_{\left(n+1\right)}]$ when $p_{\left(n\right)}<1$.
When $p_{\left(n\right)}=1$, we take $I_{n-1}$ to be
$\left[p_{\left(n-1\right)},p_{\left(n\right)}\right]$. Obviously there must
be a unique $j^{\ast}$ with $0\leq j^{\ast}\leq n$ such that $s\in
I_{j^{\ast}}$. Given $R\left(1\right)=m$ and
$\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)\in\left[0,1\right]$,
the properties of $L\left(\cdot\right)$ in Lemma 1 imply

$A_{s,r}^{\ast}=A_{s,r}\bigcap\left(\bigcap\nolimits_{j=j^{\ast}+1}^{n+1}A_{p_{\left(j\right)},r}\right)\text{.}$

Consequently $A_{s,r}^{\ast}\in\mathcal{F}_{s}$ and
$\tilde{A}_{s}\in\mathcal{F}_{s}$, which validates the first part of the
assertion.

To show $\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)=\alpha$, it
suffices to consider the case $\alpha>0$. Since $\alpha<1$ and
$R\left(t\right)>0$ by assumption, we can write
$\widetilde{FDR}_{\lambda}\left(t\right)$ as

$\widetilde{FDR}_{\lambda}\left(t\right)=\dfrac{\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)}{\rho_{m}^{\ast}\left(t\right)}\text{,}$

where
$\hat{F}_{m}\left(t\right)=m^{-1}\sum\limits_{i=1}^{m}\mathbf{1}_{\left\\{p_{i}\leq
t\right\\}}$ for $t\in\left[0,1\right]$ and

$\rho_{m}^{\ast}\left(t\right)=t^{-1}\left\\{\hat{F}_{m}\left(t\right)-\hat{F}_{m}\left(0\right)\right\\}\text{.}$

Therefore,

$\tilde{t}_{\alpha}=\sup\left\\{t\in\left[0,1\right]:\rho_{m}^{\ast}\left(t\right)\geq\alpha^{-1}\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)\right\\}\text{,}$

i.e., $\tilde{t}_{\alpha}$ is the last time the “slope”
$\rho_{m}^{\ast}\left(t\right)$ is no less than
$\alpha^{-1}\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)$.

Obviously, $\rho_{m}^{\ast}\left(\tilde{t}_{\alpha}\right)\geq\alpha$ and
there must be a unique $0\leq j^{\prime}\leq n$ such that
$\tilde{t}_{\alpha}\in I_{j^{\prime}}$. Since
$\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)>\alpha$,
$\tilde{t}_{\alpha}<1$ holds. Thus there must be some $d^{\prime}>0$ such that
$I^{\ast}=[\tilde{t}_{\alpha},\tilde{t}_{\alpha}+d^{\prime})\subseteq$
$I_{j^{\prime}}$. Noting that
$\rho_{m}^{\ast}\left(t\right)=m^{-1}L^{-1}\left(t\right)$ is continuous and
decreasing on $I_{j^{\prime}}$, we see that if
$\rho_{m}^{\ast}\left(\tilde{t}_{\alpha}\right)>\alpha$, then we can choose
$\hat{t}_{\alpha}\in I^{\ast}$ such that $\hat{t}_{\alpha}>\tilde{t}_{\alpha}$
and $\rho_{m}^{\ast}\left(\hat{t}_{\alpha}\right)>\alpha$. This contradicts
the definition of $\tilde{t}_{\alpha}.$ Hence
$\rho_{m}^{\ast}\left(\tilde{t}_{\alpha}\right)=\alpha$ must hold, and this is
equivalent to
$\widetilde{FDR}_{\lambda}\left(\tilde{t}_{\alpha}\right)=\alpha$, which
completes the proof.

###  A.3 Proof of Theorem 4

Define
$K_{\alpha}=\left\\{t\in\left[0,1\right]:\widetilde{FDR}_{\lambda}\left(t\right)<1\right\\}$
and $K_{\alpha,\delta}=K_{\alpha}\cap\left[\delta,1\right]$. Let
$FDP\left(t\right)=\dfrac{V\left(t\right)}{R\left(t\right)\vee 1}$. We verify
the claims for $t\in K_{\alpha,\delta}$ since they automatically hold when the
supremum is taken over $t\in\left[0,1\right]\backslash K_{\alpha}$. By the
Glivenko-Cantelli theorem,

$\lim_{m\rightarrow\infty}\sup_{0\leq t\leq
1}\left|m^{-1}V\left(t\right)-\pi_{0,\infty}G_{0}\left(t\right)\right|=0$ (12)

and

$\ \lim_{m\rightarrow\infty}\sup_{0\leq t\leq
1}\left|m^{-1}\left\\{R\left(t\right)\vee
1\right\\}-\left\\{\pi_{0,\infty}G_{0}\left(t\right)+\pi_{1,\infty}G_{1}\left(t\right)\right\\}\right|=0$
(13)

hold $P$-a.s.. Since
$\pi_{0,\infty}G_{0}\left(\delta\right)+\pi_{1,\infty}G_{1}\left(\delta\right)>0$
and $G_{i}\left(\cdot\right)$, $i=0,1$ are both non-decreasing functions,
using (13) it can be easily shown that

$\lim_{m\infty}\sup_{1\geq t\geq\delta}\left|\dfrac{m}{R\left(t\right)\vee
1}-\dfrac{1}{\pi_{0,\infty}G_{0}\left(t\right)+\pi_{1,\infty}G_{1}\left(t\right)}\right|=0\text{,
}P\text{-a.s..}$ (14)

Observing

$\displaystyle\widetilde{FDR}_{\lambda}\left(t\right)-FDP\left(t\right)$
$\displaystyle=$
$\displaystyle\dfrac{\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)t-\pi_{0,\infty}G_{0}\left(t\right)}{\left\\{R\left(t\right)\vee
1\right\\}/m}-\dfrac{m^{-1}V\left(t\right)-\pi_{0,\infty}G_{0}\left(t\right)}{\left\\{R\left(t\right)\vee
1\right\\}/m}$ $\displaystyle=$ $\displaystyle
I_{1}\left(t\right)-I_{2}\left(t\right)\text{,}$

we see, from (12), (13) and properties of $G_{i}\left(\cdot\right)$, $i=0,1$,
that

$\lim_{m\rightarrow\infty}\sup_{t\in
K_{\alpha,\delta}}\left|I_{2}\left(t\right)\right|\leq\lim_{m\rightarrow\infty}\sup_{t\in
K_{\alpha,\delta}}\dfrac{m}{R\left(t\right)\vee
1}\lim_{m\rightarrow\infty}\sup_{t\in
K_{\alpha,\delta}}\left|\dfrac{V\left(t\right)}{m}-\pi_{0,\infty}G_{0}\left(t\right)\right|=0\text{.}$

Since
$\lim\limits_{m\rightarrow\infty}\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)\geq\pi_{0}$,
$P$-a.s. and $G_{0}\left(t\right)\leq t$ for all $t\in(0,1]$, then

$\liminf_{m\rightarrow\infty}\inf_{1\geq
t\geq\delta}\left\\{\hat{\pi}_{0}^{G}\left(\lambda,\mathbf{\epsilon}\right)t-\pi_{0}G_{0}\left(t\right)\right\\}\geq
0$

and $\liminf\limits_{m\rightarrow\infty}\inf_{t\in
K_{\alpha,\delta}}I_{1}\left(t\right)\geq 0.$ This justifies (9).

To show (10), in view of (9) and

$\widetilde{FDR}_{\lambda}\left(t\right)-FDR\left(t\right)=\left\\{\widetilde{FDR}_{\lambda}\left(t\right)-FDP\left(t\right)\right\\}+\left\\{FDP\left(t\right)-FDR\left(t\right)\right\\}\text{,}$

it suffices to show

$\lim_{m\rightarrow\infty}\sup_{1\geq
t\geq\delta}\left|FDP\left(t\right)-FDR\left(t\right)\right|=0\text{,
}P\text{-a.s..}$ (15)

With (14), we have, $P$-a.s.,

$\displaystyle\lim_{m\rightarrow\infty}\sup_{1\geq
t\geq\delta}\left|FDP\left(t\right)-\dfrac{\pi_{0,\infty}G_{0}\left(t\right)}{\pi_{0,\infty}G_{0}\left(t\right)+\pi_{1,\infty}G_{1}\left(t\right)}\right|$
$\displaystyle\leq$ $\displaystyle\lim_{m\rightarrow\infty}\sup_{1\geq
t\geq\delta}\left|\dfrac{V\left(t\right)/m-\pi_{0,\infty}G_{0}\left(t\right)}{\left\\{R\left(t\right)\vee
1\right\\}/m}\right|$ $\displaystyle+\lim_{m\rightarrow\infty}\sup_{1\geq
t\geq\delta}\left|\left[\dfrac{1}{\left\\{R\left(t\right)\vee
1\right\\}/m}-\dfrac{1}{\pi_{0,\infty}G_{0}\left(t\right)+\pi_{1,\infty}G_{1}\left(t\right)}\right]\pi_{0,\infty}G_{0}\left(t\right)\right|$
$\displaystyle=$ $\displaystyle 0\text{.}$

By (A.3), the dominated convergence theorem and Fatou’s lemma, it follows that

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\lim_{m\rightarrow\infty}\sup_{1\geq
t\geq\delta}E\left[\left|FDP\left(t\right)-\dfrac{\pi_{0,\infty}G_{0}\left(t\right)}{\pi_{0,\infty}G_{0}\left(t\right)+\pi_{1,\infty}G_{1}\left(t\right)}\right|\right]$
$\displaystyle=$ $\displaystyle E\left[\lim_{m\rightarrow\infty}\sup_{1\geq
t\geq\delta}\left|FDP\left(t\right)-\dfrac{\pi_{0,\infty}G_{0}\left(t\right)}{\pi_{0,\infty}G_{0}\left(t\right)+\pi_{1,\infty}G_{1}\left(t\right)}\right|\right]$
$\displaystyle\geq$ $\displaystyle\lim_{m\rightarrow\infty}\sup_{1\geq
t\geq\delta}\left|FDR\left(t\right)-\dfrac{\pi_{0,\infty}G_{0}\left(t\right)}{\pi_{0,\infty}G_{0}\left(t\right)+\pi_{1,\infty}G_{1}\left(t\right)}\right|\text{.}$

Thus (15) holds, and so does (10), which completes the proof.

###  A.4 Proof of Theorem 5

Let $t^{\prime}>0$ be such that
$\alpha-\widetilde{FDR}_{\lambda}^{\infty}\left(t^{\prime}\right)=\varepsilon^{\prime}$
for some $\varepsilon^{\prime}>0$. We can take $m$ sufficiently large such
that
$\left|\widetilde{FDR}_{\lambda}^{\infty}\left(t^{\prime}\right)-\widetilde{FDR}_{\lambda}\left(t^{\prime}\right)\right|<\varepsilon^{\prime}/2$,
which implies $\widetilde{FDR}_{\lambda}\left(t^{\prime}\right)<\alpha$ and
$t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\geq t^{\prime}.$ Therefore
$\liminf\limits_{m\rightarrow\infty}t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\geq
t^{\prime}$, $P$-a.s.. By (9), it follows that, $P$-a.s.,

$\displaystyle\liminf_{m\rightarrow\infty}\left[\widetilde{FDR}_{\lambda}\left\\{t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\right\\}-FDP\left\\{t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\right\\}\right]$
$\displaystyle\geq$ $\displaystyle\liminf_{m\rightarrow\infty}\inf_{1\geq
t\geq
t^{\prime}/2}\left\\{\widetilde{FDR}_{\lambda}\left(t\right)-FDP\left(t\right)\right\\}\geq
0\text{.}$

Since
$\widetilde{FDR}_{\lambda}\left\\{t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\right\\}\leq\alpha,$
we have
$\limsup\limits_{m\rightarrow\infty}FDP\left\\{t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\right\\}\leq\alpha$,
$P$-a.s.. By Fatou’s lemma,

$\limsup_{m\rightarrow\infty}E\left[FDP\left\\{t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\right\\}\right]\leq
E\left[\limsup_{m\rightarrow\infty}FDP\left\\{t_{\alpha}\left(\widetilde{FDR}_{\lambda}\right)\right\\}\right]\leq\alpha\text{.}$

This completes the proof.

## References

  * (1)
  * Benjamini and Hochberg (1995) Benjamini, Y. and Hochberg, Y. (1995).  Controlling the false discovery rate: a practical and powerful approach to multiple testing, J. R. Statist. Soc. Ser. B 57(1): 289–300.
  * Benjamini et al. (2006) Benjamini, Y., Krieger, A. M. and Yekutieli, D. (2006).  Adaptive linear step-up procedures that control the false discovery rate, Biometrika 93(3): 491–507.
  * Benjamini and Yekutieli (2001) Benjamini, Y. and Yekutieli, D. (2001).  The control of the false discovery rate in mutliple testing under dependency, Ann. Statist. 29(4): 1165–1188.
  * Blanchard and Roquain (2009) Blanchard, G. and Roquain, E. (2009).  Adaptive false discovery rate control under independence and dependence, JMLR 10: 2837–2871.
  * Chen and Doerge (2012) Chen, X. and Doerge, R. (2012).  Generalized estimators formultiple testing: proportion of true nulls and false discovery rate, Technical report. Department of statistics, Purdue University, West Lafayette .
  * Di et al. (2011) Di, Y., Schafer, D. W., Cumbie, J. S. and Chang, J. H. (2011).  The NBP negative binomial model for assessing differential gene expression from RNA-Seq, Statistical Applications in Genetics and Molecular Biology 10(1): 24 pages.
  * Genovese and Wasserman (2002) Genovese, C. and Wasserman, L. (2002).  Operating characteristics and extensions of the false discovery rate procedure, J. R. Statist. Soc. Ser. B 64(3): 499–517.
  * Lehmann and Romano (2005) Lehmann, E. and Romano, J. P. (2005).  Testing Statistical Hypotheses, Springer, New York.
  * Liang and Nettleton (2012) Liang, K. and Nettleton, D. (2012).  Adaptive and dynamic adaptive procedures for false discovery rate control and estimation, J. R. Statist. Soc. Ser. B 74(1): 163–182.
  * Lister et al. (2008) Lister, R., O’Malley, R., Tonti-Filippini, J., Gregory, B. D., Berry, Charles C. Millar, A. H. and Ecker, J. R. (2008).  Highly integrated single-base resolution maps of the epigenome in arabidopsis, Cell 133(3): 523–536.
  * Pena et al. (2011) Pena, E. A., Habiger, J. D. and Wu, W. (2011).  Power-enhanced multiple decision functions controlling family-wise error and false discovery rates, Ann. Statist. 39(1): 556–583.
  * Pounds and Cheng (2006) Pounds, S. and Cheng, C. (2006).  Robust estimation of the false discovery rate, Bioinformatics 22(16): 1979–1987.
  * Robinson et al. (2010) Robinson, M. D., J., M. D. and Smyth, G. K. (2010).  edgeR: a Bioconductor package for differential expression analysis of digital gene expression data, Bioinformatics 26(1): 139–140.
  * Robinson and Smyth (2008) Robinson, M. D. and Smyth, G. K. (2008).  Small-sample estimation of negative binomial dispersion, with applications to SAGE data, Biostatistics 9(2): 321–332.
  * Storey (2002) Storey, J. D. (2002).  A direct approach to false discovery rates, J. R. Statist. Soc. Ser. B 64(3): 479–498.
  * Storey et al. (2004) Storey, J. D., Taylor, J. E. and Siegmund, D. (2004).  Strong control, conservative point estimation in simultaneous conservative consistency of false discover rates: a unified approach, J. R. Statist. Soc. Ser. B 66(1): 187–205.

Figure 1: The excess $\delta\left(\hat{\pi}_{0}\right)=\hat{\pi}_{0}-\pi_{0}$
of an estimator $\hat{\pi}_{0}$ of $\pi_{0}$, where the dashed line marks
$\hat{\pi}_{0}=\pi_{0}$. “BIN”, “ENT” or “FET” indicates the type of test that
produces the $p$-values to which a $\hat{\pi}_{0}$ is applied; in the vertical
strip names “pi0” refers $\pi_{0}$. “New”, “Storey”, “Benjamini” and “Pounds”
respectively refer to $\hat{\pi}_{0}^{G}\left(0.5\right)$,
$\hat{\pi}_{0}^{S}\left(0.5\right)$, $\hat{\pi}_{0}^{B}$ and
$\tilde{\pi}_{0}^{P}$. The legend are the sample standard deviations (Std Dev)
of the $\hat{\pi}_{0}$’s displayed as gradients of color intensities. The
under-estimation in $\hat{\pi}_{0}$’s when applied to $p$-values of the ENTs
is a consequence of the extremely under-estimate $\hat{\phi}$ of the
dispersion $\phi=\sigma^{-1}$ produced by edgeR; see Figure 2. Note that a
cyan boxplot indicates a $\hat{\pi}_{0}\equiv 1$. It can be seen that for
$p$-values of BINs and FETs, $\hat{\pi}_{0}^{G}$ is stable and the most
accurate, and that even though $\hat{\pi}_{0}^{G}<\pi_{0}$ may happen in these
two cases, the crucial property $E\left[\hat{\pi}_{0}\right]\geq\pi_{0}$ is
retained.

Figure 2: Estimated common dispersion $\hat{\phi}$ by edgeR, where the dashed
line marks $10^{-1}\phi$ and the true dispersion $\phi=1.451$. It can be seen
that edgeR severely under-estimates $\phi$ and consequently distorts the true
distributions of the ENTs and those of the induced $p$-values. It is clear
from our simulation study (see also Figure 1) that a very inaccurately
estimate of $\phi$ may affect drastically an estimator $\hat{\pi}_{0}$ of
$\pi_{0}$ and the adaptive MTPs that incorporate such an $\hat{\pi}_{0}$.

Figure 3: Number of rejections made by the competing MTPs. The horizontal
strip names indicate nominal FDR levels; in the vertical strip names “pi0”
refers to $\pi_{0}$. “New”, “Storey”, “Benjamini” respectively refer to the
generalized FDRestimator $\widetilde{FDR}_{0.5}$ with
$\hat{\pi}_{0}^{G}\left(0.5\right)$, Storey’s FDR estimator
$\widehat{FDR}_{0.5}$ with $\hat{\pi}_{0}^{S}\left(0.5\right)$, and the
adaptive BH procedure $FDR_{BH}^{\ast}$ with $\hat{\pi}_{0}^{B}$. The legend
indicates the type of test that produces the $p$-values. It can be seen that
$\widetilde{FDR}_{0.5}$ usually has more rejections at a nominal FDR level
$\alpha$ and that FET is not very powerful when $\alpha$ is very small.

Figure 4: Rejection thresholds of the competing MTPs. The horizontal strip
names indicate nominal FDR levels; in the vertical strip names “pi0” refers to
$\pi_{0}$. “New”, “Storey”, “Benjamini” respectively refer to the generalized
FDRestimator $\widetilde{FDR}_{0.5}$ with $\hat{\pi}_{0}^{G}\left(0.5\right)$,
Storey’s FDR estimator $\widehat{FDR}_{0.5}$ with
$\hat{\pi}_{0}^{S}\left(0.5\right)$, and the adaptive BH procedure
$FDR_{BH}^{\ast}$ with $\hat{\pi}_{0}^{B}$. The legend indicates the type of
test that produces the $p$-values. It can be seen that $\widetilde{FDR}_{0.5}$
usually has a larger rejection threshold at a nominal FDR level.

Figure 5: False discovery proportions(FDPs) of the competing MTPs. The
horizontal strip names indicate nominal FDR levels; in the vertical strip
names “pi0” refers to $\pi_{0}$. “New”, “Storey”, “Benjamini” respectively
refer to the generalized FDRestimator $\widetilde{FDR}_{0.5}$ with
$\hat{\pi}_{0}^{G}\left(0.5\right)$, Storey’s FDR estimator
$\widehat{FDR}_{0.5}$ with $\hat{\pi}_{0}^{S}\left(0.5\right)$, and the
adaptive BH procedure $FDR_{BH}^{\ast}$ with $\hat{\pi}_{0}^{B}$. The legend
indicates the type of test that produces the $p$-values. It can be seen that
the FDRs of $\widetilde{FDR}_{0.5}$ are usually well maintained at or under
the nominal levels and that $\widetilde{FDR}_{0.5}$ usually has competitive
FDPs.


The Foresight benchmark has two parameters: the initial weight given to the insertion costs in the current route $B_{t-1}$, $\theta_0^H$, and by how much this weight is decreased after every new customer insertion, $\Delta\theta^H$. We found that setting $\theta_0^H=1.0$ and $\Delta\theta^H=\frac{\theta_0^H}{\mathbb{E}[\mathcal{D}]}$ yields the best results, which is the same setting as in Yang et al., 2016.

For PPO, we study both (i) learning the second moment of the Gaussian distribution, \(\sigma\), and (ii) setting \(\sigma\) to a constant value. In Table <ref>, a \(\dag\) indicates that \(\sigma\) was learned by the actor. We maintain a learning rate for the actor that is ten times smaller than that of the critic, to guarantee that the critic's provided values are up-to-date. Both the actor and critic use a fully connected neural network with three hidden layers and ReLU activation.

Apart from the DSPO hyperparameters in Table <ref>, we use the following default settings for the CNN. The first convolutional layer has 32 filters for the synthetic case and 64 filters for the Amazon case. The second convolutional layer has two times more filters compared to the first layer. The kernel has dimension $(3,3)$ for the convolution layers and $(2,2)$ for the pooling layer. The stride of the kernel is $(1,1)$, padding is $(1,1)$, and we do not use dilation.

Hyperparameter tuning, FB = Foresight benchmark, LB = Linear benchmark.

2cSelected values

Hyperparameters     Set of values     Synthetic case     Amazon case

2* FB      $\theta^H_0$ (initial insertion weight)     $\{1.0,0.75,0.25\}$     $1.0$     $1.0$

$\Delta\theta^H$ (insertion weight update)      $\{\theta^H_0 \setminus \mathbb{E}[\mathcal{D}],0.05,0.1\}$     $\theta^H_0\setminus\mathbb{E}[\mathcal{D}]$     $\theta^H_0\setminus\mathbb{E}[\mathcal{D}]$

2* LB       Huber loss $\delta$     $\{0.5,0.75,1.0,1.35,1.5\}$     1.0     1.0

$\alpha^{LB}$ (learning rate)     $\{10^{-2},10^{-3},10^{-4}\}$     $10^{-2}$     $10^{-3}$

11* PPO       $\gamma$ (discount factor)     $\{0.9,0.99\}$     0.99     0.99

$\alpha^{PPO}_{cr}$ (learning rate critic)     $\{10^{-2},10^{-3},10^{-4},10^{-5}\}$     $10^{-4}$     $10^{-4}$

$\alpha^{PPO}_{ac}$ (learning rate actor)     $\{10^{-2},10^{-3},10^{-4},10^{-5}\}$     $10^{-5}$     $10^{-5}$

$\sigma$     $\{\dag,0.25,0.5,1\}$     0.25     0.5

Huber loss $\delta$     $\{0.5,0.75,1.0,1.35,1.5\}$     1.0     1.0

# actor NN nodes/layer     $\{8,16,32\}$     8     16

# critic NN nodes/layer     $\{8,16,32,64\}$     16     32

Batch size     $\{32,64,128\}$     128     128

$f$ (Fourier order)     $\{2,3,4\}$     3     3

Clipping factor     $\{0.1,0.2,0.3\}$     0.2     0.2

GAE $\lambda$     $\{0.9,0.95,0.99,1.0\}$     0.95     0.95

7* DSPO
M (# grids state spatial dimension)     {25,100,900,1600,3600}     100      1600

D (# layers state temporal dimension)     {1,2,3,4,5,6,7,9,10}     2     8

Batch size     $\{32,64,128\}$     64     128

$\alpha^{DSPO}$ (learning rate)      $\{10^{-2},10^{-3},10^{-4},10^{-5}\}$     $10^{-3}$     $10^{-3}$

Huber loss $\delta$     $\{0.5,0.75,1.0,1.35,1.5\}$     1.0     1.0

# NN nodes/FC layer     $\{32,64,128,256\}$     128     256

Dropout rate FC layer     $\{0,0.01,0.05,0.1\}$     0.0     0.05

§ COMPLEMENTARY FIGURES

In this section, we provide complementary figures. Figure <ref> provides a heatmap of the number of delivery stops per cell over the complete service area of Seattle on an average day. Not that the figure is tilted by $90$ degrees such that the north is oriented to the right.

Seattle heatmap of delivery stops on an average day.

Figure <ref> shows a boxplot of the online step time in milliseconds for all dynamic pricing methods. The online step time is the time required to obtain a single decision. This is relevant in light of website load times. Although the route insertion calculation is in essence relatively cheap in terms of computational effort, it does require the calculation of new travel times, i.e., in a practical situation it requires an API request. The linear benchmark, PPO, and DSPO do not require this calculation and, hence, provide faster results. We note that in absolute terms the differences between the policies are relatively small.

boxplot/draw direction=y,
/pgfplots/boxplot/box extend=0.4,
xticklabels=Hindsight, Foresight, Linear, PPO, DSPO,
ylabel = Step time (ms),
ylabel near ticks,
cycle list=black!85,
every axis plot/.append style=fill,fill opacity=0.75,
x tick label style=
text width=2.5cm,

+ [
area legend,boxplot prepared=
upper quartile=2.8332,
lower quartile=1.57577,
upper whisker=4.24387,
lower whisker=0.2386,
draw position=1
] coordinates ;

+ [
area legend,boxplot prepared=
upper quartile=3.5304,
lower quartile=1.97537,
upper whisker=4.64327,
lower whisker=0.30386,
draw position=2
] coordinates ;

+ [
area legend,boxplot prepared=
upper quartile=1.216435,
lower quartile=0.250655,
upper whisker=2.18336,
lower whisker=0.016432,
draw position=3
] coordinates ;

+ [
area legend,boxplot prepared=
upper quartile=0.927435,
lower quartile=0.2283,
upper whisker=2.18321,
lower whisker=0.00642,
draw position=4
] coordinates ;

+ [
area legend,boxplot prepared=
upper quartile=1.417635,
lower quartile=0.3709185,
upper whisker=2.48316,
lower whisker=0.02642,
draw position=5
] coordinates ;

Step time in milliseconds of the studied dynamic pricing methods (Amazon case).