#  A unified class of null proportion estimators with plug-in FDR control

Sebastian Döhler  Iqraa Meah

###### Abstract

Since the work of Storey et al., (2004), it is well-known that the performance
of the Benjamini-Hochberg (BH) procedure can be improved by incorporating
estimators of the number (or proportion) of null hypotheses, yielding an
adaptive BH procedure which still controls FDR. Several such plug-in
estimators have been proposed since then, for some of these, like Storey’s
estimator, plug-in FDR control has been established, while for some others,
e.g. the estimator of Pounds and Cheng, (2006), some gaps remain to be closed.
In this work we introduce a unified class of estimators, which encompasses
existing and new estimators and unifies proofs of plug-in FDR control using
simple convex ordering arguments. We also show that any convex combination of
such estimators once more yields estimators with guaranteed plug-in FDR
control. Additionally, the flexibility of the new class of estimators also
allows incorporating distributional informations on the $p$-values. We
illustrate this for the case of discrete tests, where the null distributions
of the $p$-values are typically known. In that setting, we describe two
generic approaches for adapting any estimator from the general class to the
discrete setting while guaranteeing plug-in FDR control. While the focus of
this paper is on presenting the generality and flexibility of the new class of
estimators, we also include some analyses on simulated and real data.

False discovery rate,

Adaptive BH procedure,

Convex ordering of random variables,

Discrete hypothesis testing,

Multiple hypothesis testing,

###### keywords:

###### Contents

  1. 1 Introduction
    1. 1.1 Background
    2. 1.2 Contributions
  2. 2 Framework
    1. 2.1 Distributional assumptions
    2. 2.2 FDR control for plug-in estimates
  3. 3 A unified class of plug-in estimators
  4. 4 Homogeneous estimators
    1. 4.1 Numerical results
    2. 4.2 More details on the Pounds and Cheng estimator
  5. 5 Adjusted estimators for discrete $p$-values
    1. 5.1 Transformations of discrete $p$-values
    2. 5.2 Adjusting the rescaling constants
    3. 5.3 A randomization approach
    4. 5.4 Simulation results
    5. 5.5 Real data analysis
  6. 6 Discussion
  7. A Auxiliary definitions and results
  8. B Complements to Section 4.1
  9. C Complements to Section 4.2
  10. D Additional Figures for simulated data of Section 4
  11. E Upper and lower bounds for the inverse moment of the uniform sum distribution

##  1 Introduction

###  1.1 Background

When many statistical tests are performed simultaneously, an ubiquituous way
to account for the false rejections is the false discovery rate (FDR), that is
the expected proportion of errors among the rejections. The seminal Benjamini
and Hochberg, (1995) procedure (abbreviated in the sequel as BH procedure)
works by rejecting $H_{(1)},\ldots,H_{(\hat{k})}$, where $\hat{k}$ is
determined in the following step-up manner

$\displaystyle\hat{k}$
$\displaystyle=\max\left\\{\ell\in\\{0,\ldots,m\\}:p_{(\ell)}\leq\frac{\ell}{m}\cdot\alpha\right\\},$
(1)

where $p_{(1)}\leq\ldots\leq p_{(m)}$ denote the ordered $p$-values, and
$H_{(1)},\ldots,H_{(m)}$ the corresponding null hypotheses and $p_{(0)}:=0$.
According to results in Benjamini and Hochberg, (1995) and Benjamini and
Yekutieli, (2001), this procedure guarantees that
$\textnormal{FDR}\leq\pi_{0}\alpha$ when the $p$-values are independent or
positively dependent, while for arbitrarily dependent $p$-values the Benjamini
and Yekutieli, (2001) procedure is available. The simplicity of the BH
procedure and its many useful theoretical properties have made it an
indispensable tool in modern high dimensional data analysis, see e.g.
Benjamini, (2010). Much work has gone into analyzing, extending and adapting
this procedure to various settings.

In this context, Storey et al., (2004) showed that the plug-in BH procedure

$\displaystyle\hat{k}$
$\displaystyle=\max\left\\{\ell\in\\{0,\ldots,m\\}:p_{(\ell)}\leq\frac{\ell}{\widehat{m}_{0}}\cdot\alpha\right\\},$
(2)

obtained by replacing $m$ in (1) by an estimate $\widehat{m_{0}}$ of $m_{0}$
still provides so called adaptive or plug-in FDR control at level $\alpha$
while allowing for more power. Classical examples of such estimates were
proposed by Storey et al., (2004)

$\displaystyle\widehat{m}_{0}^{\text{Storey}}$
$\displaystyle=\frac{1+\sum_{i=1}^{m}{\mathbf{1}\\{p_{i}>\lambda\\}}}{1-\lambda},$
(3) where $\lambda\in[0,1)$ is a tuning parameter, and by Pounds and Cheng,
(2006) $\displaystyle\widehat{m}_{0}^{\text{PC,2006}}$
$\displaystyle=1\wedge\left(2\sum_{i=1}^{m}p_{i}\right).$ (4)

While the focus of this work is on plug-in FDR control, estimates of $m_{0}$
(or equivalently $\pi_{0}=m_{0}/m$) can also be used for FDR estimation
purposes (see Storey, (2002)). Thus, there is a large body of literature on
this topic and numerous methods for establishing plug-in FDR control are
available, see e.g. Benjamini et al., (2006); Sarkar, (2008); Blanchard and
Roquain, (2009); Heesen and Janssen, (2016); Ditzhaus and Janssen, (2019) and
references therein. In this paper, we use a classical condition proposed by
Blanchard and Roquain, (2009) for establishing plug-in FDR control for a
unified class of $m_{0}$-estimators, see Section 2.2 for more details.

Previous work on plug-in FDR control has focused on continuous test
statistics, for which the null $p$-values are distributed according to the
uniform distribution. Considering the abundance of super-uniform $p$-values in
real life applications, the uniformity assumption is often violated which may
lead to undesirable conservatism of the $m_{0}$-estimators, see Section 5 for
more details. Super-uniform $p$-values can be observed when testing composite
null hypotheses or when dealing with discrete tests, the latter being the
setting of our interest in this work.

Discrete tests often originate when the tests are based on counts or
contingency tables: for example in clinical studies, the efficacy or safety of
drugs is determined by counting patients who survive a certain period, or
experience a certain type of adverse drug reaction after being treated, see
e.g. Chavant et al., (2011); and also in biology, where the genotype effect on
the phenotype can be analyzed by knocking out genes and counting the number of
individuals with a changed phenotype, see e.g. Muñoz-Fuentes et al., (2018).
In discrete testing, each $p$-value is super-uniform and (potentially) has its
own support, thus producing heterogeneous $p$-values. Pounds and Cheng, (2006)
recognized the need for developing methods tailored to discrete $p$-values and
introduced $\widehat{m}_{0}^{\text{PC,2006}}$ as a simple and robust
$m_{0}$-estimate in this setting. They did not, however, provide a proof of
plug-in FDR control, not even in the uniform setting. Further works adressing
the discreteness and heterogeneity include Chen et al., (2018) who introduced
a $m_{0}$-estimator for discrete $p$-values based on averaging Storey type
estimators for plug-in control. Biswas and Chattopadhyay, (2020) pointed out
an error in the proof of Chen et al., (2018) and provided a corrected version.
However it is unclear whether this estimator actually provides an improvement
over the classical (uniform) Storey estimator in practice. Thus, there are
gaps to be filled on $m_{0}$-estimation both for the uniform and discrete
case.

###  1.2 Contributions

In this paper we address some of the gaps and limitations mentioned above by
introducing a simple and flexible class of $m_{0}$ estimators which has the
following properties:

  * •

Plug-in FDR control is guaranteed for all estimators contained in this class
under independence of $p$-values. We give a unified proof using simple convex
ordering arguments (for the reader’s convenience we restate some definitions
and classical results on stochastic and convex ordering in Appendix A).

  * •

It provides a simple and flexible generic formulation which is useful for
designing new estimators. In particular, we obtain a simple modification of
$\widehat{m}_{0}^{\text{PC,2006}}$ with guaranteed plug-in FDR control.

  * •

Additional distributional information on the $p$-values like heterogeneity and
super-uniformity can be incorporated easily into estimators from the class. In
particular, the estimators of this class can be used in conjunction with
classical discrete $p$-value transformations like the mid-p transformation.

  * •

Combining several weighted estimators from the class preserves plug-in FDR
control.

The paper is organized as follows: Section 2 presents the statistical setting
and restates a classical sufficient criterion for plug-in FDR control. Section
3 introduces the new class of estimators, and presents the main mathematical
results on plug-in FDR control, followed by some numerical results in Section
4. In Section 5 we present approaches for adjusting estimators to
discreteness, and investigate their performance on simulated and real data.
The paper concludes with a discussion in Section 6. Technical details –
including classical results on stochastic and convex ordering– and further
analyses are deferred to the Appendices.

##  2 Framework

###  2.1 Distributional assumptions

We use a classical setting for multiple testing encompassing homogeneous and
heterogeneous nulls, see e.g. Döhler et al., (2018). We observe $X$, defined
on an abstract probabilistic space, valued in an observation space
$(\mathcal{X},\mathfrak{X})$ and generated by a distribution $P$ that belongs
to a set $\mathcal{P}$ of possible distributions. We consider $m$ null
hypotheses for $P$, denoted $H_{0,i}$, $1\leq i\leq m$, and we denote the
corresponding set of true null hypotheses by $\mathcal{H}_{0}(P)=\\{1\leq
i\leq m\>:\>\mbox{$H_{0,i}$ is satisfied by $P$}\\}$. We also denote by
$\mathcal{H}_{1}(P)$ the complement of $\mathcal{H}_{0}(P)$ in
$\\{1,\dots,m\\}$ and by $m_{0}(P)=m_{0}=|\mathcal{H}_{0}(P)|$ the number of
true nulls.
We assume that there exists a set of $p$-values that is a set of random
variables $\\{p_{i}(X),1\leq i\leq m\\}$, valued in $[0,1]$. We introduce the
following dependence assumptions between the $p$-values:

All the $p$-values $\\{p_{i}(X)$, $1\leq i\leq m\\}$ are mutually independent
in the model $\mathcal{P}$. (Indep)

The (maximum) null cumulative distribution function (c.d.f) of each $p$-value
is denoted

$F_{i}(t)=\sup_{P\in\mathcal{P}\>:\>i\in\mathcal{H}_{0}(P)}\\{\mathbf{P}_{X\sim
P}(p_{i}(X)\leq t)\\},\>\>t\in[0,1],\>\>1\leq i\leq m.$ (5)

We assume that the set of c.d.f $\mathcal{F}=\\{F_{i},1\leq i\leq m\\}$ is
known and we consider the following possible situations for the functions in
$\mathcal{F}$:

For all $i\in\\{1,\dots,m\\}$, $F_{i}$ is continuous on $[0,1]$ (Cont)
$\displaystyle\begin{array}[]{c}\mbox{For all $i\in\\{1,\dots,m\\}$, there
exists some finite set $\mathcal{A}_{i}\subset[0,1]$ such that}\\\
\mbox{$F_{i}$ is a step function, right continuous, that jumps only at some
points of $\mathcal{A}_{i}$.}\end{array}$ (Discrete)

The case (Discrete) typically arises when for all $P\in\mathcal{P}$ and
$i\in\\{1,\dots,m\\}$, $\mathbf{P}_{X\sim P}(p_{i}(X)\in\mathcal{A}_{i})=1$
for some given finite sets $\mathcal{A}_{i}\subset[0,1]$. Throughout the
paper, we will assume that we are either in the case (Cont) or (Discrete) and
we denote $\mathcal{A}=\cup_{i=1}^{m}\mathcal{A}_{i}$, with by convention
$\mathcal{A}_{i}=[0,1]$ when (Cont) holds. We will also make use of the
following classical assumption:

$\displaystyle\mbox{For all $i\in\\{1,\dots,m\\}$, $F_{i}(t)\leq t$ for all
$t\in[0,1]$}.$ (SuperUnif)

In this paper we will always assume that the $p$-values are mutually
independent ((Indep) holds) and super-uniform under the null ((SuperUnif)
holds).

###  2.2 FDR control for plug-in estimates

The following Theorem is a central result for plug-in FDR control, providing a
sufficient condition based on bounding the inverse moment of the estimator
$\widehat{m}_{0}$ by the inverse of $m_{0}$. Our presentation follows
Blanchard and Roquain, (2009), similar results can be found in Benjamini et
al., (2006); Sarkar, (2008); Zeisel et al., (2011).

######  Theorem 2.1.

Let $\widehat{m}_{0}=\widehat{m}_{0}(p_{1},\ldots,p_{m})$ be a coordinatewise
non-decreasing function of the $p$-values $(p_{1},\ldots,p_{m})$. Assume that
$(p_{1},\ldots,p_{m})$ are mutually independent (Indep) and (SuperUnif). For
$h\in\mathcal{H}_{0}$, denote by $p_{0,h}$ the set of $p$-values where $p_{h}$
has been replaced by $0$. If

$\displaystyle\mathbb{E}\left(\frac{1}{\widehat{m}_{0}(p_{0,h})}\right)$
$\displaystyle\leq\frac{1}{m_{0}}$ (IMC)

holds for all $h\in\mathcal{H}_{0}$, then the plug-in BH procedure given by
(2) controls FDR at level $\alpha$.

Throughout this paper, we will only consider coordinatewise non-decreasing
estimators and will always assume that the $p$-values are mutually
independent. Thus, when additionally (SuperUnif) holds true, the inverse
moment criterion (IMC) is sufficient for establishing plug-in FDR control in
our proofs. We mostly present results in term of the absolute number of null
hypotheses $m_{0}$, but clearly equivalent statements using the proportion of
null hypotheses $\pi_{0}=m_{0}/m$ hold, and in some cases we will present
results in terms of $\pi_{0}$ instead of $m_{0}$.

##  3 A unified class of plug-in estimators

In this section, we introduce a new class of estimators for $m_{0}$ (or
equivalently $\pi_{0}$) that mathematically guarantees plug-in FDR control. It
is based on sums of suitably transformed $p$-values, allowing us to recover
classical estimators, such as the Storey (3) and the PC (slightly modified)
(4) estimators, and also to define new estimators. We first present a general
result for single estimators and then show that plug-in FDR control is also
preserved for convex combinations of such estimators. To start, assume that
the $p$-values are transformed by certain functions $g\in\mathcal{G}$, with

$\displaystyle\mathcal{G}$ $\displaystyle=\\{g:[0,1]\rightarrow[0,1]:\text{
$g$ is non-decreasing and $\mathbb{E}g(U)>0$, where
$U\sim\mathcal{U}[0,1]$}\\}.$ (8)

Accordingly we define the class of estimators $\mathcal{F}_{0}$ as

$\displaystyle\mathcal{F}_{0}$
$\displaystyle=\left\\{\widehat{m}_{0}:[0,1]^{m}\rightarrow[0,\infty)|\>\widehat{m}_{0}(p_{1},\ldots,p_{m})=\frac{1}{{\nu(g)}}\left(1+\sum_{i=1}^{m}g(p_{i})\right),g\in\mathcal{G}\right\\},$
(9)

where $\nu(g)=\mathbb{E}g(U)$ for any $g\in\mathcal{G}$ with
$U\sim\mathcal{U}[0,1]$ (for brevity we sometimes omit the $g$ in $\nu$ when
there is no ambiguity concerning the function $g$). The class
$\mathcal{F}_{0}$ contains the classical estimator
$\widehat{m}_{0}^{\text{Storey}}$ (3) by taking
$g(u)={\mathbf{1}\\{u>\lambda\\}}$ and $\nu=1-\lambda$. It also contains a
slightly modified version $\widehat{m}_{0}^{\text{PC,new}}$ of the classical
estimator $\widehat{m}_{0}^{\text{PC,2006}}$ (4) obtained from taking $g(u)=u$
with $\nu=1/2$, i.e.

$\displaystyle\widehat{m}_{0}^{\text{PC,new}}$
$\displaystyle=2+2\sum_{i=1}^{m}p_{i}.$ (10)

In Section 4 we will introduce some additional estimators and discuss
$\widehat{m}_{0}^{\text{PC,new}}$ in more detail. The rationale behind the
definitions of the classes $\mathcal{G}$ and $\mathcal{F}_{0}$ is two-fold.
Requiring that $g$ is non-decreasing ensures that $\widehat{m}_{0}$ is
coordinatewise non-decreasing, allowing us to apply Theorem 2.1. The quantity
$g(p_{i})/\nu$ can be interpreted as the (local) contribution of $p_{i}$ to
the estimate of $m_{0}$. If we expect large $p$-values to provide evidence for
null hypotheses, then it seems reasonable to require $g$ to be non-decreasing.
Rescaling $g(p_{i})$ by $\nu=\mathbb{E}g(U)$ is a simple way of ensuring that
$\sum_{i=1}^{m}g(p_{i})/\nu$ is conservatively biased in the sense that
$\mathbb{E}(\sum_{i=1}^{m}g(p_{i})/\nu)\geq m_{0}$ in any constellation of
null and alternative hypotheses. This type of conservativeness may however not
be strong enough for plug-in control. As our main result – Proposition 3.2
below – shows, simply adding $1/\nu$ as a ’safety margin’ to the above
estimate is enough for ensuring plug-in FDR control.
In some situations $p$-values under the null may be heterogeneous, i.e. the
$p$-values may have different distributions under the null, so that using an
individual transformation for each $p$-value may be helpful. To this end, we
introduce the following richer and more flexible class of estimators.

$\begin{split}\mathcal{F}=\biggl{\\{}\widehat{m}_{0}:[0,1]^{m}&\rightarrow[0,\infty)\
\bigg{|}\
\widehat{m}_{0}(p_{1},\ldots,p_{m})=\frac{1}{\min(\nu_{1},\ldots,\nu_{m})}\\\
&+\sum_{i=1}^{m}\frac{g_{i}(p_{i})}{\nu_{i}},\quad\text{with
$g_{i}\in\mathcal{G}$ and $\nu_{i}=\mathbb{E}[g_{i}(U)]$,
$U\sim\mathcal{U}[0,1]$ for all $i$}\biggr{\\}}.\end{split}$ (11)

We state our main result on plug-in FDR control for this more general class
below. Clearly, $\mathcal{F}_{0}\subset\mathcal{F}$, so that the results
stated for $\mathcal{F}$ also hold for $\mathcal{F}_{0}$. We now present an
upper bound in the convex order for transformed uniform random variables in
terms of Bernoulli random variables, which is the main technical tool we use
for proving plug-in FDR control for the class $\mathcal{F}$.

######  Lemma 3.1.

For any $g\in\mathcal{G}$ we have
$g(U)\leqslant_{\text{cx}}\mathbf{Bin}(1,\nu)$ and $\nu=\mathbb{E}g(U)$, where
$U\sim\mathcal{U}[0,1]$, and the notation $\leqslant_{\text{cx}}$ denotes the
convex ordering (see Definition A.2).

###### Proof.

For $U\sim\mathcal{U}[0,1]$ define $X=g(U)$, so that $\mathbb{E}(X)=\nu$. Let
$l_{X}=\inf_{x\in[0,1]}g(x)$, $u_{X}=\sup_{x\in[0,1]}g(x)$ be the lower and
upper endpoints of the support of $X$, and define a two-point distribution $Y$
concentrated on $\\{l_{X},u_{X}\\}$ by $P(Y=l_{X})=(u_{X}-\nu)/(u_{X}-l_{X})$
and $P(Y=u_{X})=(\nu-l_{X})/(u_{X}-l_{X})$. By Lemma A.3 we have
$X\leqslant_{\text{cx}}Y$.
Now let $Z\sim\mathbf{Bin}\left(1,\nu\right)$ and denote the distribution
function of $Y$ and $Z$ by $F$ and $G$. Clearly,
$\mathbb{E}Y=\nu=\mathbb{E}Z$. Since $[l_{X},u_{X}]\subset[0,1]$ and both $Y$
and $Z$ are two-point distributions, the function $G-F$ posesses one crossing
point on $[0,1)$. Indeed, for $t\in[0,l_{X}),\quad F(t)=0$ while $G(t)=1-\nu$
so that $G-F$ is positive, and for $t\in[u_{X},1),F(t)=1$ while $G(t)=1-\nu$
so that $G-F$ is negative. For $t\in[l_{X},u_{X})$, $G-F$ can be positive or
negative depending on $\nu$. Overall, the sign sequence of $G-F$ is $+,-$ so
that Lemma A.4 implies that $Y\leqslant_{\text{cx}}Z$ and the claim follows. ∎

The following proposition is our main result on plug-in FDR control.

######  Proposition 3.2.

Assume that $p_{1},\ldots,p_{m}$ are mutually independent and (SuperUnif)
holds. Then (IMC) holds true for any estimator
$\widehat{m}_{0}\in\mathcal{F}$, where $\mathcal{F}$ is defined by (11). In
particular, the BH plug-in procedure (2) using $\widehat{m}_{0}$ controls FDR
at level $\alpha$.

###### Proof.

Since $\widehat{m}_{0}$ is coordinatewise non-decreasing, it is sufficient to
verify (IMC). For any $h\in\mathcal{H}_{0}$, monotonicity and super-uniformity
give us $\widehat{m}_{0}(p_{0,h})\geqslant_{\text{st}}1/\nu+S_{0}$, where
$\nu=min_{l\in\mathcal{H}_{0}\setminus\\{h\\}}\nu_{l}$, and
$S_{0}=\sum_{\ell\in\mathcal{H}_{0}\setminus\\{h\\}}g_{\ell}(U_{\ell})/\nu_{\ell}$
with $(U_{\ell})_{\ell\in\mathcal{H}_{0}}$ i.i.d random variables distributed
according to $\mathcal{U}[0,1]$. By Lemma 3.1 we have
$g_{\ell}(U_{\ell})\leqslant_{\text{cx}}\mathbf{Bin}(1,\nu_{\ell})$ and Lemma
A.6 gives
$\mathbf{Bin}(1,\nu_{i})/\nu_{i}\leqslant_{\text{cx}}\mathbf{Bin}(1,\nu)/\nu$.
Since the convex ordering is preserved under convolutions (see Lemma A.5) we
obtain $\nu S_{0}\leqslant_{\text{cx}}\mathbf{Bin}(m_{0}-1,\nu)$. Finally, the
mapping $x\mapsto\nu/(1+x)$ is convex on $[0,\infty)$ and therefore from the
Definition A.2 of $\leqslant_{\text{cx}}$ we obtain that

$\displaystyle\mathbb{E}\left(\frac{1}{\widehat{m}_{0}(p_{0,h})}\right)$
$\displaystyle\leq\mathbb{E}\left(\frac{1}{\frac{1}{\nu}+S_{0}}\right)=\mathbb{E}\left(\frac{\nu}{1+\nu
S_{0}}\right)\leq\mathbb{E}\left(\frac{\nu}{1+\mathbf{Bin}(m_{0}-1,\nu)}\right)\leq\frac{1}{m_{0}},$
(12)

where the last bound is a well-known result for the inverse moment of Binomial
distributions (see e.g. A.7 in Appendix) so that (IMC) is proved. The
statement on plug-in FDR control now follows from Theorem 2.1. ∎

By taking $g(u)={\mathbf{1}\\{u>\lambda\\}}$ we have $\nu
S_{0}\sim\mathbf{Bin}(m_{0}-1,\nu)$ and therefore the second inequality from
the right in (12) can be replaced by an equality. Thus, it may be tempting to
conclude that $\widehat{m}_{0}^{\text{Storey}}$ is optimal. In the case of a
Dirac-Uniform constellation of $p$-values (see Blanchard and Roquain, (2009))
this is indeed true, since $\nu\widehat{m}_{0}^{\text{Storey}}\sim
1+\mathbf{Bin}(m_{0}-1,\nu)$ and therefore the left inequality in (12) can
also be replaced by an equality. In more general settings however, other
choices of $g$ may be better, as the results in Section 4.1 show.

We highlight that introducing a general class of estimators as (11) with
Proposition 3.2 allows a unified proof of plug-in FDR control for known
estimators like $\widehat{m}_{0}^{\text{Storey}}$ and
$\widehat{m}_{0}^{\text{PC,new}}$ and also for new estimators that we will
define in Section 4. Additionally, the classes $\mathcal{F}_{0}$ and
$\mathcal{F}$ possess stability properties that make it possible to combine
various plug-in estimators while maintaining FDR control.

######  Proposition 3.3.

Let $\widehat{m}_{1},\widehat{m}_{2}\in\mathcal{F}$, where $\mathcal{F}$ is
defined by (11) and let $\lambda\in[0,1]$. Then the BH plug-in procedure (2)
using $\widehat{m}_{0}=\lambda\widehat{m}_{1}+(1-\lambda)\widehat{m}_{2}$
controls FDR at level $\alpha$.

###### Proof.

We show that $\widehat{m}_{0}$ satisfies (IMC). Let
$\widehat{m}_{1},\widehat{m}_{2}\in\mathcal{F}$ have the representation

$\displaystyle\widehat{m}_{1}$
$\displaystyle=\frac{1}{\nu}+\sum_{i=1}^{m}\frac{g_{i}(p_{i})}{\nu_{i}}\qquad\text{and}\qquad\widehat{m}_{2}=\frac{1}{\mu}+\sum_{i=1}^{m}\frac{h_{i}(p_{i})}{\mu_{i}},$
where $\nu=min(\nu_{1},\ldots,\nu_{m})$ and $\mu=min(\mu_{1},\ldots,\mu_{m})$
so that $\displaystyle\widehat{m}_{0}$
$\displaystyle=\frac{\lambda}{\nu}+\frac{1-\lambda}{\mu}+\sum_{i=1}^{m}\frac{\lambda
g_{i}(p_{i})}{\nu_{i}}+\frac{(1-\lambda)h_{i}(p_{i})}{\mu_{i}}$

and define weights
$\kappa_{i}=\frac{\lambda\mu_{i}}{\lambda\mu_{i}+(1-\lambda)\nu_{i}}$ and
transformations $f_{i}=\kappa_{i}g_{i}+(1-\kappa_{i})h_{i}$. Clearly,
$\kappa_{i}\in[0,1]$ and $f_{i}\in\mathcal{G}$ and we introduce
$\epsilon_{i}=\mathbb{E}(f_{i})=\kappa_{i}\nu_{i}+(1-\kappa_{i})\mu_{i}$. From
the above definitions we obtain with some straightforward algebra

$\displaystyle\lambda$
$\displaystyle=\frac{\kappa_{i}\nu_{i}}{\epsilon_{i}}\qquad\text{and}\qquad
1-\lambda=\frac{(1-\kappa_{i})\mu_{i}}{\epsilon_{i}}$ (13) which yields
$\displaystyle\frac{\lambda}{\nu_{i}}g_{i}$
$\displaystyle+\frac{(1-\lambda)}{\mu_{i}}h_{i}=\frac{\kappa_{i}}{\epsilon_{i}}g_{i}+\frac{(1-\kappa_{i})}{\epsilon_{i}}h_{i}=\frac{f_{i}}{\epsilon_{i}}.$
(14) From (13) we have $\displaystyle\frac{\lambda}{\nu}$
$\displaystyle=\max\left(\frac{\lambda}{\nu_{1}},\ldots,\frac{\lambda}{\nu_{m}}\right)=\max\left(\frac{\kappa_{1}}{\epsilon_{1}},\ldots,\frac{\kappa_{m}}{\epsilon_{m}}\right)\qquad\text{and}$
$\displaystyle\frac{1-\lambda}{\mu}$
$\displaystyle=\max\left(\frac{1-\lambda}{\mu_{1}},\ldots,\frac{1-\lambda}{\mu_{m}}\right)=\max\left(\frac{1-\kappa_{1}}{\epsilon_{1}},\ldots,\frac{1-\kappa_{m}}{\epsilon_{m}}\right)$
so that the sub-additivity of the $\max$ function now yields the bound
$\displaystyle\frac{\lambda}{\nu}$
$\displaystyle+\frac{1-\lambda}{\mu}\geq\frac{1}{\epsilon},$ (15)

where $\epsilon=\min(\epsilon_{1},\ldots,\epsilon_{m})$. Combining (14) and
(15) now gives us

$\displaystyle\widehat{m}_{0}$
$\displaystyle=\frac{\lambda}{\nu}+\frac{1-\lambda}{\mu}+\sum_{i=1}^{m}\frac{\lambda
g_{i}(p_{i})}{\nu_{i}}+\frac{(1-\lambda)h_{i}(p_{i})}{\mu_{i}}\geq\frac{1}{\epsilon}+\sum_{i=1}^{m}\frac{f_{i}(p_{i})}{\epsilon_{i}}:=\widetilde{m}_{0}$

with $\widetilde{m}_{0}\in\mathcal{F}$ and from Proposition 3.2 we know that
(IMC) holds true for $\widetilde{m}_{0}$ and therefore also for
$\widehat{m}_{0}$. ∎

The proof shows that $\mathcal{F}$ is “almost” convex in the sense that
whenever equality holds in (15) we have
$\widehat{m}_{0}=\widetilde{m}_{0}\in\mathcal{F}$. If
$\widehat{m}_{1},\widehat{m}_{2}\in\mathcal{F}_{0}$, i.e. each estimator uses
only a single transformation function then it is easy to see that equality
holds in (15) which leads to the following result:

######  Proposition 3.4.

The class of estimators $\mathcal{F}_{0}$ given by (9) is convex. In
particular this implies that for any
$\widehat{m}_{1},\widehat{m}_{2}\in\mathcal{F}_{0}$ and $\lambda\in[0,1]$ the
BH plug-in procedure (2) controls FDR at level $\alpha$ for the estimator
$\widehat{m}_{0}=\lambda\widehat{m}_{1}+(1-\lambda)\widehat{m}_{2}$.

##  4 Homogeneous estimators

In this section we focus on the class of homogeneous estimators
$\mathcal{F}_{0}$ given by (9), i.e. on estimators of the form

$\displaystyle\widehat{m}_{0}$
$\displaystyle=\widehat{m}_{0}(p_{1},\ldots,p_{m})=\frac{1}{\nu}\left(1+\sum_{i=1}^{m}g(p_{i})\right),$

where $g\in\mathcal{G}$ and $\nu=\nu(g)=\mathbb{E}g(U)>0$, with
$U\sim\mathcal{U}[0,1]$. As mentioned before, this class includes the
classical estimator $\widehat{m}_{0}^{\text{Storey}}$ (3) and the new
estimator $\widehat{m}_{0}^{\text{PC,new}}$ (10), and also gives the scientist
freedom to define new estimators with proven plug-in FDR control thanks to
Proposition 3.2 There are many conceivable ways in which this can be done. As
an ad hoc example, we define a polynomial estimator of degree $r\geq 0$ and
thresholding parameter $\lambda\in[0,1)$ by taking $\widehat{m}_{0}$ as above
in (9) with $g(u)=g_{r,\lambda}(u)=u^{r}\cdot{\mathbf{1}\\{u>\lambda\\}}$, so
that $\nu=\frac{1-\lambda^{r+1}}{r+1}$. This gives us

$\displaystyle\widehat{m}_{0}^{\text{Poly}}(r,\lambda)=\frac{r+1}{1-\lambda^{r+1}}+\frac{r+1}{1-\lambda^{r+1}}\sum_{i=1}^{m}p_{i}^{r}\cdot{\mathbf{1}\\{p_{i}>\lambda\\}}.$
(16)

It is easily seen that the classical estimators
$\widehat{m}_{0}^{\text{Storey}}$ and $\widehat{m}_{0}^{\text{PC,new}}$ are
particular instances of $\widehat{m}_{0}^{\text{Poly}}(r,\lambda)$ with $r=0$
for $\widehat{m}_{0}^{\text{Storey}}$ and $r=1$ and $\lambda=0$ for
$\widehat{m}_{0}^{\text{PC,new}}$. Taking $r=1$ and $\lambda>0$ yields a
hybrid estimator which combines $\widehat{m}_{0}^{\text{Storey}}$ and
$\widehat{m}_{0}^{\text{PC,new}}$ which has the potential to combine the
strengths of both methods. For all estimators
$\widehat{m}_{0}^{\text{Poly}}(r,\lambda)$ plug-in FDR control follows
immediately from Proposition 3.2. For illustrational purposes we effectively
only use $r$ as a parameter and set the thresholding parameter to the
classical value of $\lambda=1/2$ throughout the paper. These examples are
primarily meant to illustrate the freedom and flexibility Proposition 3.2
allows for the class $\mathcal{F}_{0}$ and should not be interpreted as
recommendations for optimal choices. These examples are investigated further
in the following section.

###  4.1 Numerical results

Here we compare the performance of several estimators from the class
$\mathcal{F}_{0}$ in a Gaussian one-sided testing setting. We assume that we
observe $X_{1},\dots,X_{m}$ independent random variables with
$X_{i}\sim\mathcal{N}(0,1)$ for $i\in\mathcal{H}_{0}$ and
$X_{i}\sim\mathcal{N}(\mu,1)$ for $i\in\mathcal{H}_{1}$, with $\mu>0$, and we
test $H_{0,i}:\mu=0$ vs. $H_{1,i}:\mu>0$. For a given signal strength $\mu>0$
under the alternative, closed-form expressions for the expectation and
variance of $\widehat{m}_{0}$ are available, see Appendix B for more details.
Thus we can numerically compare the mean squared error (MSE) of estimators
from $\mathcal{F}_{0}$.

For this analysis, we fix $m=10\>000$ and first compare the MSE w.r.t. to the
signal strength $\mu$ of the alternatives with a fixed proportion of true
nulls $\pi_{0}=0.6$. Then, we compare the MSE w.r.t. the proportion of true
nulls $\pi_{0}$ with a fixed signal strength $\mu=1.5$. The considered
estimators for this comparison are $\widehat{m}_{0}^{\text{Storey}}$ (3),
$\widehat{m}_{0}^{\text{PC,new}}$ (10),
$\widehat{m}_{0}^{\text{Poly}}(1,1/2)$, and
$\widehat{m}_{0}^{\text{Poly}}(2,1/2)$ (see (16) for both). The MSE is
evaluated in terms of $\pi_{0}$ for better readability and displayed in Figure
1. The qualitative comparison between the estimators remains consistent across
both panels of Figure 1: $\widehat{m}_{0}^{\text{PC,new}}$ has the poorest
performance, characterized by the largest MSE, followed by
$\widehat{m}_{0}^{\text{Storey}}$. While the polynomial approach shows some
benefits, the improvement is not particularly remarkable except for small to
moderate values of $\pi_{0}$. For larger values of $\pi_{0}$ or $\mu$, there
are no noticeable differences in performance.

Figure 1: MSE against signal strength $\mu\in[0,3]$ with $\pi_{0}=0.6$(left)
and MSE against $\pi_{0}\in[0.1,0.9]$ with $\mu=1.5$(right)

###  4.2 More details on the Pounds and Cheng estimator

While FDR control for $\widehat{m}_{0}^{\text{Storey}}$ is a classical result
following from Theorem 2.1 (Blanchard and Roquain, (2009); Benjamini et al.,
(2006)), much less is known about the validity of
$\widehat{m}_{0}^{\text{PC,2006}}$ as a plug-in estimator. Indeed, Pounds and
Cheng, (2006) introduced their estimator (4) primarily to obtain a robust
estimate of FDR. To the best of our knowlegde, the only related result on
plug-in FDR control was obtained by Zeisel et al., (2011), who defined the
following modified version of $\widehat{m}_{0}^{\text{PC,2006}}$:

$\displaystyle\widehat{m}_{0}^{\text{PC,ZZD}}$
$\displaystyle=C(m)\cdot\min\left[m,\max\left(s(m),2\cdot\sum_{i=1}^{m}p_{i}\right)\right],$
(17)

where the correction factors $C(m)$ and $s(m)$ are chosen in such a way so
that (IMC) holds. However, determining these factors is non-trivial and
requires extensive use of numerical integration and approximations methods
(see Supplement B in Zeisel et al., (2011) for further details) so that no
simple representation of $C(m)$ and $s(m)$ is available (for selected values
of $m$, Table S1 in Zeisel et al., (2011) lists values for the correction
factors).

By contrast, our new modification (10) is extremely simple and, as we show in
Section 5, can be adapted easily to e.g. discrete tests, thus confirming a
conjecture in Pounds and Cheng, (2006). Its validity for plug-in FDR control
follows directly from Proposition 3.2 and involves no sophisticated asymptotic
or numerical approximations. Supplementary material in Appendix C shows that
the two versions of the PC estimator behave more or less identically.
Nevertheless, we argue in favor of using $\widehat{m}_{0}^{\text{PC,new}}$
since it is both conceptually and computationally much simpler than
$\widehat{m}_{0}^{\text{PC,ZZD}}$.

##  5 Adjusted estimators for discrete $p$-values

In this section we assume – additionally to mutual independence of the
$p$-values – that the null distribution functions $F_{1},\ldots,F_{m}$ are
known. As a particular application we consider the setting of discrete
$p$-values (see Section 1.1 for more detailed references). The classical plug-
in estimators, like the Storey, (2002) estimator defined in (3), were
developed for uniformly distributed $p$-values under the nulls, and can thus
suffer of an inflated bias when computed under (SuperUnif) assumption.

To illustrate the problem, we compare the bias of an arbitrary estimator
$\widehat{m}_{0}\in\mathcal{F}_{0}$ under the uniform setting with the bias
under the super-uniform setting. In the classical uniform case, considering
marginally independent $p$-values $p_{i}\sim X_{0}$ for $i\in\mathcal{H}_{0}$,
and $p_{i}\sim X_{1}$ for $i\in\mathcal{H}_{1}$, for some variables
$X_{0},X_{1}$ defined on $[0,1]$, the bias is seen to be

$\displaystyle\textnormal{Bias}[\widehat{m}_{0}]$
$\displaystyle=\mathbb{E}[\widehat{m}_{0}]-m_{0}=\frac{1}{\nu}(1+m_{1}\mathbb{E}[g(X_{1})]).$
(18) In contrast, under the super-uniform setting, still considering
independent $p$-values under the null and the alternative, the bias is
$\displaystyle\textnormal{Bias}[\widehat{m}_{0}]$
$\displaystyle=\frac{1}{\nu}(1+m_{1}\mathbb{E}[g(X_{1})])+\frac{1}{\nu}m_{0}(\mathbb{E}[g(X_{0})]-\nu).$
(19)

Recall that $\nu=\mathbb{E}[g(U)]$ with $U\sim\mathcal{U}[0,1]$, thus under
super-uniformity $\mathbb{E}[g(X_{0})]\geq\nu$ (see the characterization of
the usual stochastic order in Appendix A), which shows that an additional
source of conservativeness is present in this case. In general, practitioners
use classical estimators without worrying about $p$-values distributions,
ingenuously expecting the estimator to perform according to the “uniform” bias
(18) when in fact it often performs according to the “super-uniform” bias
(19). This motivates the need for a correction in the estimator that will aim
at deflating (19).
Super-uniformity does not solely appear in the discrete setting, it also
occurs e.g. when testing composite nulls, however in the discrete setting
additional information on the $p$-values c.d.f (as defined by (5)) may be
available and leveraged to correct the over-conservativeness. In this section,
we present two such approaches that incorporate the available knowledge of
$F_{i}$ – the $p$-value c.d.f under the null – in the estimators. The standard
way of defining $p$-values for discrete tests leads to distribution functions
that satisfy (Discrete) and (SuperUnif). As we later introduce transformed
$p$-values with transformed distribution, for clearer distinguishability the
c.d.f associated to these standard discrete $p$-values are denoted by
$F^{\text{sd}}_{1},\ldots,F^{\text{sd}}_{m}$ (where the upper-script “sd”
denotes “standard discrete”).

###  5.1 Transformations of discrete $p$-values

In order to reduce the individual conservatism of $p$-values caused by super-
uniformity, various transformation of discrete $p$-values have been proposed,
see e.g. Habiger, (2015). Perhaps the most popular transformation is the so-
called mid-$p$-value. For the realization $x$ of the random variable $X$, let
$p(x)$ be the (realized) standard $p$-value. Now define the mid-$p$-value
(Rubin-Delanchy et al.,, 2019) $q(x)$ given the observation $x$ as

$\displaystyle q(x)$ $\displaystyle=p(x)-\frac{1}{2}P_{0}(p(X)=p(x)),$ (20)

where $P_{0}$ denotes the distribution of $X$ under the null (for simplicity,
we assume that such a unique distribution exists). Transforming the $p$-value
through (20) helps to mimic the behavior of a uniform random variable in
expectation. Indeed, we always have $\mathbb{E}[q(X)]=1/2$, see Berry and
Armitage, (1995) for more details. However, the distribution of the
mid-$p$-value is no longer super-uniform but shrunk toward 0 as displayed in
Figure 2. In what follows, we denote by
$F^{\text{mid}}_{1},\ldots,F^{\text{mid}}_{m}$ the distribution functions of
the mid-$p$-values associated with the distribution functions
$F^{\text{sd}}_{1},\ldots,F^{\text{sd}}_{m}$ of the standard $p$-values. In
Section 5.2 we show how the distribution functions of standard discrete or
mid-$p$-values can be used in $m_{0}$-estimators introduced in Section 3,
while preserving plug-in FDR control.

Another transformation to reduce the conservativeness of discrete $p$-values
uses so-called randomized $p$-values which are defined in our context by

$\displaystyle r(x,u)$ $\displaystyle=p(x)-u\cdot P_{0}(p(X)=p(x)),$ (21)

where $u$ is the realization of a uniform random variable
$U\sim\mathcal{U}[0,1]$, independent of $X$. Alternatively to the notation
$r(x,u)$, we will also use (with a slight abuse) $r(p,u)$, where $p=p(x)$ is
the standard $p$-value obtained from observation $x$. Randomized $p$-values
and mid-$p$ values are related via the conditional expectation on the
observations $q(x)=\mathbb{E}_{U}[r(x,U)|X=x]$. Randomization leads to an
(unconditional) uniform behavior, however at the cost of introducing an
additional source of randomness which makes its use controversial for
decisions on individual hypotheses, see e.g. Habiger and Pena, (2011) for a
discussion. We show in Section 5.3 that for estimation purposes however,
randomized $p$-values can be beneficial for obtaining an efficient non-
randomized estimator.

Figure 2: Distribution functions of a standard $p$-value and its associated
mid-$p$-value.The standard $p$-value has support $\\{0.3,0.55,0.7,0.9,1\\}$
(orange rug plot), the c.d.f (orange solid line) is super-uniform, i.e. below
the uniform c.d.f (gray line). The support of the associated mid-$p$-value
$\\{0.15,0.465,0.625,0.8,0.95\\}$ (green rug plot) is shifted to the left, but
the probabilities (given by the jumps in the green c.d.f) remain the same.

###  5.2 Adjusting the rescaling constants

The first approach for adjusting estimators to discrete $p$-values is tailored
to estimators from the class $\mathcal{F}_{0}$ and adjusts the rescaling
constant $\nu$ in $\widehat{m}_{0}$. In fact, this approach is not limited to
the discrete setting, and can also be applied for arbitrary $p$-value
distributions.

######  Proposition 5.1.

Assume that $p_{1},\ldots,p_{m}$ are mutually independent and the null
distribution functions $F_{1},\ldots,F_{m}$ are known. For any
$\widehat{m}_{0}\in\mathcal{F}_{0}$ (9) with

$\displaystyle\widehat{m}_{0}(p_{1},\ldots,p_{m})$
$\displaystyle=\frac{1}{\nu(g)}\left(1+\sum_{i=1}^{m}g(p_{i})\right)$ (22)
define the adjusted estimator
$\displaystyle\widehat{m}^{\text{adj}}_{0}(p_{1},\ldots,p_{m})$
$\displaystyle=\frac{1}{min(\nu^{\text{adj}}_{1},\ldots,\nu^{\text{adj}}_{m})}+\sum_{i=1}^{m}\frac{g(p_{i})}{\nu^{\text{adj}}_{i}}$
(23)

where $\nu^{\text{adj}}_{i}=\mathbb{E}_{p_{i}\sim F_{i}}[g(p_{i})]$, is the
expectation of the transformed p-value taken w.r.t. $F_{i}$. Then the BH plug-
in procedure (2) using $\widehat{m}^{\text{adj}}_{0}$ controls FDR at level
$\alpha$.

###### Proof.

Without loss of generality we assume that $\nu^{\text{adj}}_{i}>0$, otherwise
the sum and minimum in (23) is to be taken over the index set
$\\{i:\nu^{\text{adj}}_{i}>0\\}$.

For any $i\in\\{1,\dots,m\\}$ define $g_{i}:[0,1]\rightarrow[0,1]$ by
$g_{i}(y)=g\circ F^{-1}_{i}(y)$ for $y\in(0,1]$, where
$F^{-1}_{i}(y)=\inf\\{x\in\mathbb{R}:F_{i}(x)\geq y\\}$ is the generalized
inverse of $F_{i}$, and set $g_{i}(0)=g(0)$. Since $g\in\mathcal{G}$ and
$F^{-1}_{i}$ are both nondecreasing, so is $g_{i}$. For $i\in\mathcal{H}_{0}$,
with $U\sim\mathcal{U}[0,1]$, we have $p_{i}\sim F^{-1}_{i}(U)$ by Proposition
2 in Embrechts and Hofert, (2013), so that $g_{i}(U)\sim g(p_{i})$, which
implies that $\mathbb{E}[g_{i}(U)]=\mathbb{E}_{p_{i}\sim
F_{i}}[g(p_{i})]=\nu^{\text{adj}}_{i}$, so that (23) belongs to the class
$\mathcal{F}$.
Now let $z_{1},\dots z_{m}$ be independent random variables with
$z_{i}\sim\mathcal{U}[0,1]$ for $i\in\mathcal{H}_{0}$ and
$z_{i}\sim\delta_{0}$ for $i\in\mathcal{H}_{1}$ (Dirac-Uniform configuration).
Since $g_{i}(U)\sim g(p_{i})$ for $i\in\mathcal{H}_{0}$, we have

$\displaystyle\widetilde{m}^{\text{adj}}_{0}(z_{1},\ldots,z_{m})$
$\displaystyle=\frac{1}{min(\nu^{\text{adj}}_{1},\ldots,\nu^{\text{adj}}_{m})}+\sum_{i\in\mathcal{H}_{0}}\frac{g_{i}(z_{i})}{\nu^{\text{adj}}_{i}}$
$\displaystyle\sim\frac{1}{min(\nu^{\text{adj}}_{1},\ldots,\nu^{\text{adj}}_{m})}+\sum_{i\in\mathcal{H}_{0}}\frac{g(p_{i})}{\nu^{\text{adj}}_{i}}\leq\widehat{m}^{\text{adj}}_{0}(p_{1},\ldots,p_{m})\quad\mbox{(a.s.)}$

Since $\widetilde{m}^{\text{adj}}_{0}(z_{1},\ldots,z_{m})\in\mathcal{F}$, by
Proposition 3.2 (IMC) holds for
$\widetilde{m}^{\text{adj}}_{0}(z_{1},\ldots,z_{m})$, and since
$\widetilde{m}^{\text{adj}}_{0}(z_{1},\ldots,z_{m})\leq\widehat{m}^{\text{adj}}_{0}(p_{1},\ldots,p_{m})$
(a.s.), (IMC) also holds for
$\widehat{m}^{\text{adj}}_{0}(p_{1},\ldots,p_{m})$. ∎

Following Proposition 5.1, we define the discrete-uniform estimator using (23)
with standard discrete $p$-values $p_{1},\ldots,p_{m}$ and their distribution
functions $F^{\text{sd}}_{1},\ldots,F^{\text{sd}}_{m}$

$\displaystyle\widehat{m}_{0}^{\text{du}}(p_{1},\ldots,p_{m})$
$\displaystyle=\frac{1}{min(\nu^{\text{du}}_{1},\ldots,\nu^{\text{du}}_{m})}+\sum_{i=1}^{m}\frac{g(p_{i})}{\nu^{\text{du}}_{i}},$
(24)

where $\nu^{\text{du}}_{i}=\mathbb{E}_{p_{i}\sim
F^{\text{sd}}_{i}}[g(p_{i})]$.

######  Corollary 5.1.

Assume that $p_{1},\ldots,p_{m}$ are mutually independent and (SuperUnif)
holds with null distribution functions
$F^{\text{sd}}_{1},\ldots,F^{\text{sd}}_{m}$ that are known. Then the BH plug-
in procedure (2) using $\widehat{m}_{0}^{\text{du}}$ as in (24)controls FDR at
level $\alpha$. Moreover $\widehat{m}_{0}^{\text{du}}\leq\widehat{m}_{0}$
(a.s.), where $\widehat{m}_{0}$ is the base non-adjusted estimator (22).

The last statement of the corollary shows that for standard discrete
$p$-values the estimator $\widehat{m}_{0}^{\text{du}}$ is guaranteed to
perform better than $\widehat{m_{0}}$. This follows from the fact that
$\nu^{\text{du}}_{1},\ldots,\nu^{\text{du}}_{m}\geq\nu$ because $g$ is non-
decreasing and (SuperUnif) holds (see Appendix A).

For classical estimators, the adjusted rescaling constants can be computed
easily, using

  * •

$\nu^{\text{du-Storey}}_{i}=1-F^{\text{sd}}_{i}(\lambda)$;

  * •

$\nu^{\text{du-PC}}_{i}=\sum_{x\in\mathcal{A}_{i}}x\cdot P(p_{i}=x)$, where
$\mathcal{A}_{i}$ denotes the support of $F^{\text{sd}}_{i}$.

Similarly to (24), we define a mid $p$-value estimator using (23) with
mid-$p$-values $q_{1},\ldots,q_{m}$ and their distribution functions
$F^{\text{mid}}_{1},\ldots,F^{\text{mid}}_{m}$

$\displaystyle\widehat{m}_{0}^{\text{mid}}(q_{1},\ldots,q_{m})$
$\displaystyle=\frac{1}{min(\nu^{\text{mid}}_{1},\ldots,\nu^{\text{mid}}_{m})}+\sum_{i=1}^{m}\frac{g(q_{i})}{\nu^{\text{mid}}_{i}}$
(25)

where $\nu^{\text{mid}}_{i}=\mathbb{E}_{q_{i}\sim
F^{\text{mid}}_{i}}[g(q_{i})]$ is the expectation taken under the null using
the mid $p$-value distribution function $F^{\text{mid}}_{i}$. For
$\widehat{m}_{0}^{\text{Storey}}$ we have $\nu^{\text{mid-
Storey}}_{i}=1-F^{\text{mid}}_{i}(\lambda)\leq
1-F^{\text{sd}}_{i}(\lambda)=\nu^{\text{du-Storey}}_{i}$ and
$g(q_{i})={\mathbf{1}\\{q_{i}>\lambda\\}}\leq{\mathbf{1}\\{p_{i}>\lambda\\}}=g(p_{i})$
so that $\widehat{m}_{0}^{\text{mid}}$ can be smaller or larger than
$\widehat{m}_{0}^{\text{du}}$, depending on the specific constellation. In the
case of the PC estimator we have $g(x)=x$ and since $\mathbb{E}q_{i}=1/2$ for
any mid-$p$-value (see Berry and Armitage, (1995)) we have
$\nu^{\text{mid}}_{1}=\ldots=\nu^{\text{mid}}_{m}=1/2$ so that in this case
the mid-$p$ estimator has a particularly simple representation. Combining this
with the fact that $q_{i}\leq p_{i}$ $(a.s.)$ gives us the following result.

######  Corollary 5.2.

Assume that $p_{1},\ldots,p_{m}$ are mutually independent and super-uniform
under the null (i.e. (SuperUnif) holds), and let $q_{1},\ldots,q_{m}$ denote
the corresponding mid-$p$-values. Then the mid-$p$ estimator of
$\widehat{m}_{0}^{\text{PC,new}}$ is given by

$\displaystyle\widehat{m}_{0}^{\text{mid-PC}}(q_{1},\ldots,q_{m})$
$\displaystyle=2+2\cdot\sum_{i=1}^{m}q_{i}$

and the BH plug-in procedure (2) using $\widehat{m}_{0}^{\text{mid-
PC}}(q_{1},\ldots,q_{m})$ controls FDR at level $\alpha$. Moreover,
$\widehat{m}_{0}^{\text{mid-PC}}\leq\widehat{m_{0}}$ $(a.s.)$, where
$\widehat{m}_{0}$ is the base non-adjusted estimator (22).

This result implies that for the PC estimator with discrete data we can simply
use $2+2\cdot\sum_{i=1}^{m}q_{i}$ instead of the more conservative
$2+2\cdot\sum_{i=1}^{m}p_{i}$ estimator without losing plug-in FDR control. We
point out that the mid-$p$-values are used exclusively for estimating $m_{0}$
in the plug-in procedure defined by (2) while the (ordered) standard discrete
$p$-values $p_{(k)}$ are used in the final decision step.

###  5.3 A randomization approach

Here we describe an approach related to Dickhaus et al., (2012) who argue for
using randomization methods in estimating $m_{0}$ on discrete data. For any
estimator $\widehat{m}_{0}$, not necessarily belonging to $\mathcal{F}_{0}$
define the associated expected randomized estimator as

$\displaystyle\widehat{m}_{0}^{\text{rand}}(p_{1},\ldots,p_{m})$
$\displaystyle=\left[\mathbb{E}_{(U_{1},\ldots,U_{m})}\left(\frac{1}{\widehat{m}_{0}(r(p_{1},U_{1}),\ldots,r(p_{m},U_{m}))}\right)\right]^{-1}$
(26)

where $U_{1},\ldots,U_{m}\sim\mathcal{U}[0,1]$ denote i.i.d uniform random
variables independent of $(p_{1},\ldots,p_{m})$. Thus, for fixed
$(p_{1},\ldots,p_{m})$ this estimator is obtained by taking the expectation
over the randomized $p$-values associated with $(p_{1},\ldots,p_{m})$. In most
cases (26) is analytically intractable, we therefore use Monte-Carlo
approximation of $\widehat{m}_{0}^{\text{rand}}$ obtained by averaging a large
number of simulations of
$\widetilde{m}_{0}(r(p_{1},U_{1}),\ldots,r(p_{m},U_{m}))$ (the vector
$(U_{1},\ldots,U_{m})$ is simulated many times, while $(p_{1},\ldots,p_{m})$
is kept fixed). Again, this approach comes with guaranteed FDR plug-in
control.

######  Corollary 5.3.

Assume that $p_{1},\ldots,p_{m}$ are mutually independent and super-uniform
under the null (i.e. (SuperUnif) holds) and let $\widehat{m}_{0}$ satisfy the
conditions of Theorem 2.1. Then the BH plug-in procedure (2) using
$\widehat{m}_{0}^{\text{rand}}(p_{1},\ldots,p_{m})$ defined by (26) controls
FDR at level $\alpha$. Moreover
$\widehat{m}_{0}^{\text{rand}}\leq\widehat{m}_{0}$ (a.s.).

###### Proof.

The proof uses Theorem 2.1. First, we show that
$\widehat{m}_{0}^{\text{rand}}(p_{1},\ldots,p_{m})$ is coordinatewise non-
decreasing. For fixed $(u_{1},\ldots,u_{m})\in[0,1]^{m}$ each (realized)
randomized $p$-value $r_{i}=r(p_{i},u_{i})$ is non-decreasing in $p_{i}$.
Since $\widehat{m}_{0}\in\mathcal{F}$ is coordinatewise non-decreasing in
$(p_{1},\ldots,p_{m})$, the function
$1/\widehat{m}_{0}(r(\cdot,u_{1}),\ldots,r(\cdot,u_{m}))$ is coordinatewise
decreasing for all $(u_{1},\ldots,u_{m})\in[0,1]^{m}$ and so is its
expectation which implies that $\widehat{m}_{0}^{\text{rand}}$ is
coordinatewise non-decreasing. To establish (IMC), we denote for
$h\in\mathcal{H}_{0}$ by $r_{0,h}$ the set of randomized $p$-values
$(r_{1},\ldots,r_{m})$, where $r_{h}$ has been replaced by $0$. By the
definition of $\widehat{m}_{0}^{\text{rand}}$ we have

$\displaystyle\mathbb{E}_{(p_{1},\ldots,p_{m})}\left[\frac{1}{\widehat{m}_{0}^{\text{rand}}(p_{0,h})}\right]$
$\displaystyle=\mathbb{E}_{(p_{1},\ldots,p_{m})}\left[\mathbb{E}_{(U_{1},\ldots,U_{m})}\frac{1}{\widehat{m}_{0}(r_{0,h})}\right]=\mathbb{E}_{(r_{1},\ldots,r_{m})}\left[\frac{1}{\widehat{m}_{0}(r_{0,h})}\right].$

where the second equality follows from the fact that for super-uniform
$p$-value $p_{h}=0$, the associated randomized $p$-value $r(p_{h},u)=0$ (a.s)
by Definition (21). Since the $(r_{1},\ldots,r_{m})$ are mutually independent
and uniform under the null and $\widehat{m}_{0}\in\mathcal{F}$, the bound
(IMC) for $\widehat{m}_{0}(r_{0,h})$ now follows since $\widehat{m}_{0}$
satisfies the conditions of Theorem 2.1. Therefore, the r.h.s. of the last
equation can be bounded by $1/m_{0}$ and plug-in FDR control for
$\widehat{m}_{0}^{\text{rand}}$ now follows from Theorem 2.1. To see that the
last statement of the corollary holds true, observe that since
$\widehat{m}_{0}$ is coordinatewise non-decreasing and $r(p_{i},U_{i})\leq
r(p_{i},0)=p_{i}$ we have
$\widehat{m}_{0}(r(p_{1},U_{1}),\ldots,r(p_{m},U_{m}))\leq\widehat{m}_{0}(p_{1},\ldots,p_{m})$
and therefore the r.h.s. of (26) is bounded by
$\widehat{m}_{0}(p_{1},\ldots,p_{m})$ $(a.s.)$. ∎

Dickhaus et al., (2012) argue for using randomized $p$-values in (essentially)
Storey’s estimator, i.e. applying $\widehat{m}_{0}^{\text{Storey}}$ to
$(r_{1},\ldots,r_{m})$ instead of $(p_{1},\ldots,p_{m})$ which yields a random
estimate that should provide a better estimate for $m_{0}$. They show that
plugging this estimator into the Bonferroni procedure yields asymptotic
control of the Familywise Error Rate (FWER) under certain assumptions. They
also point out that if fully reproducible results are desired it may be more
appropriate to work with the conditional expectation w.r.t. randomization,
i.e. using
$\mathbb{E}_{(U_{1},\ldots,U_{m})}\left(\widehat{m}_{0}^{\text{Storey}}(r_{1},\ldots,r_{m})\right)$.
Corollary 5.3 shows that we can obtain similar guarantees w.r.t. to plug-in
FDR control in a finite-sample setting for any estimator $\widehat{m}_{0}$
satisfying the conditions of Theorem 2.1 and in particular for
$\widehat{m}_{0}\in\mathcal{F}_{0}$ by using conditional expectation w.r.t.
randomization. The slightly complicated form of (26) is a natural consequence
of Theorem 2.1, but if the variance of
$\widehat{m}_{0}(r(p_{1},U_{1}),\ldots,r(p_{m},U_{m}))$ w.r.t.
$U_{1},\ldots,U_{m}$ is small we have the approximation
$\widehat{m}_{0}^{\text{rand}}(p_{1},\ldots,p_{m})\approx\mathbb{E}_{(U_{1},\ldots,U_{m})}\widehat{m}_{0}(r(p_{1},U_{1}),\ldots,r(p_{m},U_{m}))$.

###  5.4 Simulation results

In this section, we analyze how the discrete adjustments can improve the base
estimators on simulated data. More specifically, we follow Döhler et al.,
(2018) by simulating a two-sample problem in which a vector of $m=500$
independent binary responses is observed for $N=25$ subjects in both groups.
The goal is to test the $m$ null hypotheses $H_{0,i}$: ’$p_{1i}=p_{2i}$’,
$i=1,...,m$ where $p_{1i}$ and $p_{2i}$ are the success probabilities for the
$i^{th}$ binary response in group A and B respectively. Thus, for each
hypothesis $i$, the data can be summarized by a $2\times 2$ contingency table,
and we use (two-sided) Fisher’s exact tests (FETs) for testing $H_{0i}$. The
$m$ hypotheses are split in three groups of size $m_{1}$, $m_{2}$, and $m_{3}$
such that $m=m_{1}+m_{2}+m_{3}$. Then, the binary responses are generated as
i.i.d Bernoulli of probability 0.01 ($\mathbf{Bin}(1,0.01)$) at $m_{1}$
positions for both groups, i.i.d $\mathbf{Bin}(1,0.10$) at $m_{2}$ positions
for both groups, and i.i.d $\mathbf{Bin}(1,0.10)$ at $m_{3}$ positions for one
group and i.i.d $\mathbf{Bin}(1,p_{3})$ at $m_{3}$ positions for the other
group. Thus, null hypotheses are true for $m_{1}+m_{2}$ positions, while they
are false for $m_{3}$ positions (set $\mathcal{H}_{1}$). We interpret $p_{3}$
as the strength of the signal and set it to 0.4, while
$\pi_{1}=\frac{m_{3}}{m}$, corresponds to the proportion of signal. Also,
$m_{1}$ and $m_{2}$ are both taken equal to $\frac{m-m_{3}}{2}$.

We first compare the base estimators $\widehat{m}_{0}^{\text{Storey}}$ (3),
$\widehat{m}_{0}^{\text{PC,new}}$ (10), and
$\widehat{m}_{0}^{\text{Poly}}(2,1/2)$ (16) with their standard discrete
rescaled versions. Figure 3 displays the estimation results for a grid of true
$\pi_{0}\in\\{0.1,\dots,0.9\\}$. We can see that incorporating discreteness
leads to considerable improvements for all estimators over the entire range of
$\pi_{0}$ values. This is particularly relevant for large values of $\pi_{0}$
where the base estimators may lead to a strong deterioration in the power of
the plug-in BH procedure. On another note, among the base estimators we can
see that $\widehat{m}_{0}^{\text{Poly}}(2,1/2)$ performs poorly compared to
the results of Section 4.1. This seems plausible since a large portion of
$p$-values are equal to 1 in the discrete setting in contrast to the Gaussian
setting. For these $p$-values, the contribution in the
$\widehat{m}_{0}^{\text{Poly}}(2,1/2)$ estimator is equal to the constant
$\nu=\frac{3}{1-1/2^{3}}=\frac{24}{7}$ which is much larger than the
corresponding contribution of $\nu=2$ in the Storey estimator.

Figure 3: Comparison between base estimators
$\widehat{m}_{0}^{\text{Storey}}$, $\widehat{m}_{0}^{\text{PC,new}}$ and
$\widehat{m}_{0}^{\text{Poly}}(2,1/2)$ and their standard discrete rescaled
versions on simulated data.

In a second step, we compare the different discrete adjustments
$\widehat{m}_{0}^{\text{du}}$ as in (24), $\widehat{m}_{0}^{\text{mid}}$ as in
(25), and $\widehat{m}_{0}^{\text{rand}}$ as in (26) in Figure 4, where we
display the estimation results for three values of true
$\pi_{0}\in\\{0.2,0.5,0.7\\}$. We can see that there are no relevant
differences between the different adjustments. Therefore, there is no strong
reason to advocate a specific type of adjustment since they yield similar
outcomes.

Figure 4: Comparison between discrete adjustments
$\widehat{m}_{0}^{\text{du}}$, $\widehat{m}_{0}^{\text{mid}}$ and
$\widehat{m}_{0}^{\text{rand}}$ on simulated data.

###  5.5 Real data analysis

Finally, we compare the performance of base and discrete estimators on three
different datasets. The first dataset consists of data provided by the
International Mice Phenotyping Consortium (IMPC) (Karp et al.,, 2017), which
coordinates studies on the genotype influence on mouse phenotype. This dataset
includes, for each of the $m=266952$ studied genes, the counts of normal and
abnormal phenotypes thus providing multiple two by two contingency tables,
which can be analysed using FETs. Then we analyze the methylation dataset for
cytosines of Arabidopsis in Lister et al. (2008) which is part of the
R-package fdrDiscreteNull of Chen and Doerge, (2015). This dataset contains
$m=3525$ counts for a biological entity under two different biological
conditions or treatments also analyzed using FETs. Finally, the third dataset,
provided by the Regulatory Agency in the United Kingdom, includes adverse drug
reactions due to medicines and healthcare products. It contains the number of
reported cases of amnesia as well as the total number of adverse events
reported for each of the $m=2446$ drugs in the database. For more details we
refer to Heller and Gur, (2011) and to the accompanying R-package discreteMTP
of Heller et al., (2012), which also contains the data. Heller and Gur, (2011)
investigate the association between reports of amnesia and suspected drugs by
performing for each drug (one-sided) FETs.

From the results in Table 1 we can see that taking discreteness into account
is always beneficial, regardless of the adjustment used. Depending on the type
of discreteness and the amount of signal contained in the data, adjusting for
discreteness can provide a great improvement in some cases. Indeed, as the
example of the IMPC data shows, base estimators may not be able to recognize
the presence of any alternatives. However, the discrete estimators clearly
suggest that a considerable amount of alternatives is present.

Table 1: $\pi_{0}$-estimates for base estimators and adjusted discrete estimators on three different datasets containing discrete data. |  | Dataset
---|---|---
Adjustment | Estimator | IMPC | Arabidopsis | Pharmacovigilance
standard (none) | Storey | 1.26 | 0.67 | 1.79
| PC | 1.26 | 0.73 | 1.79
| Poly$(2,1/2)$ | 2.16 | 0.75 | 2.97
rescaled (du) | Storey | 0.63 | 0.59 | 1.05
| PC | 0.63 | 0.64 | 1.04
| Poly$(2,1/2)$ | 0.63 | 0.57 | 1.10
rescaled (mid) | Storey | 0.63 | 0.63 | 1.03
| PC | 0.63 | 0.64 | 1.05
| Poly$(2,1/2)$ | 0.63 | 0.58 | 1.11
randomized | Storey | 0.63 | 0.58 | 1.08
| PC | 0.63 | 0.64 | 1.06
| Poly$(2,1/2)$ | 0.63 | 0.56 | 1.14

##  6 Discussion

In this paper we introduce a unified class of $m_{0}$-estimators with
mathematical guarantees on plug-in FDR control in a classical setting. We also
describe some general approaches for adjusting $m_{0}$-estimators constructed
for continuous $p$-values to discrete $p$-values. While we show that these
results are useful both from a methodological viewpoint and for practical
purposes, there a numerous possibilities for further investigations, some of
which we describe now.

While we focus on FDR control in this paper, it is clear, that our new
estimators can also be used for FDR estimation, see Storey, (2002). For the
discrete estimators from Section 5 that are uniformly better than their
classical counterparts, this implies that the corresponding FDR estimators are
uniformly better es well.

In Section 4, we describe the performance of various polynomial estimators in
a one-sided Gaussian testing framework with the aim of illustrating the
flexibility of our Proposition 3.2 on plug-in FDR control. The numerical
results in Section 4.1 show that $\widehat{m}_{0}^{\text{Poly}}(2,1/2)$
performs uniformly better in terms of MSE than the other estimators in this
specific framework. In a second step, it might be interesting to study whether
optimal estimators can be derived, either within the whole class
$\mathcal{F}_{0}$ or perhaps within some sub-class like polynomial estimators.
This way, it may eventually be possible to obtain more efficient estimators in
practice, or at least give the user some guidance for choosing estimators from
the class $\mathcal{F}_{0}$.

We would also like to point out some connections to the work of Heesen and
Janssen, (2016), who split the unit interval into an estimation region on
which an estimator of $m_{0}$ is constructed and a rejection region on which
the BH procedure is run. Thus these estimators do not use all available
$p$-values, in contrast to our approach. Heesen and Janssen, (2016) derive a
general sufficient criterion, similar to Theorem 2.1, for finite sample plug-
in FDR control which they apply to Storey-type estimators and histogram-type
estimators (see MacDonald et al., (2019)). In contrast to our approach, the
transformation function $g$ applied to $p$-values need not be monotone,
however it is unclear whether e.g. smooth functions fit into this framework.
Their approach also accommodates “dynamization”, which allows data-dependent
tuning of parameters, see MacDonald et al., (2019). We may wonder if this
approach can be used for the estimators of our class, but as this question
exceeds the scope of this paper, we leave it for future research.

While we use polynomial estimators as simple examples that include both the
classical Storey and Pounds and Cheng estimators, other choices are
conceivable. Taking, for instance, certain kernel-type transformation
functions would lead to estimators that are advantageous from an asymptotic
viewpoint, see Neuvial, (2013). We leave this topic for future research.

In Section 5 we illustrate how information on the null distribution functions
of discrete $p$-values can be used to obtain more efficient
$m_{0}$-estimators. This information can be seen as a special case of
auxiliary covariates, for which it is well-known that their incorporation into
multiple testing procedures e.g. by weighting, can be highly beneficial (see
Ignatiadis et al., (2016); Durand, (2019)). We would like to mention that our
methods for $m_{0}$-estimation are not limited to the special case of discrete
$p$-values but should be able to accomodate these types of heterogeneity as
well. Thus, they might also be useful in the settings described above for
obtaining more efficient $m_{0}$-estimators.

While we assume independence of the $p$-values throughout this paper, it is
well-known that dependence may adversely affect the performance of
$m_{0}$-estimators or may require re-adjustments of tuning parameters like
$\lambda$ in $\widehat{m}_{0}^{\text{Storey}}$ (see e.g. Blanchard and
Roquain, (2009)). Thus, it might be interesting to investigate the behaviour
of our new estimators under various types of dependency.

Constructing multiple testing procedures for discrete data that provide finite
sample plug-in FDR control is challenging. In this paper we make some progress
by obtaining improved discrete estimators for $m_{0}$. While using these
discrete estimators in the plug-in BH procedure provide more power than using
classical estimators (based on uniformly distributed $p$-values under the
null), it is still not ideal because the discreteness is ignored in the
rejection stage of the procedure. Döhler et al., (2018) propose discrete
variants of the standard (i.e. non plug-in) BH procedure. They also sketch a
possible plug-in method based on combining this procedure with estimators of
$m_{0}$, but caution that it comes without mathematical guarantees. Thus, as
MacDonald et al., (2019) pointed out, it still remains an open problem to
develop procedures that integrate discreteness of the data in both the
estimation of $m_{0}$ and the rejection of $p$-values.

## Acknowledgements

This work is part of project DO 2463/1-1, funded by the Deutsche
Forschungsgemeinschaft. The authors thank Etienne Roquain for insightful
discussions and helpful suggestions which significantly improved the
manuscript.

##  Appendix A Auxiliary definitions and results

In this appendix we recall some definitions and results of stochastic ordering
following the presentation in [SS], to which we also refer the reader for
further details. We also recall a well-known bound on the inverse moment of
the Binomial distribution.

######  Definition A.1 (Stochastic order).

Let $X$ and $Y$ be two random variables such that

$\displaystyle\mathbf{P}(X>x)\leq\mathbf{P}(Y>x)\quad\text{for all
}x\in(-\infty,\infty),$

Then $X$ is said to be smaller than $Y$ in the usual stochastic order denoted
by $X\leqslant_{\text{st}}Y$.

An equivalent characterization of the stochastic order is that
$X\leqslant_{\text{st}}Y$ $\Leftrightarrow$
$\mathbb{E}[g(X)]\leq\mathbb{E}[g(Y)]$, for all non-decreasing functions
$g:\mathbb{R}\rightarrow\mathbb{R}$ for which the expectations exist (see
(1.A.7) in [SS])).

######  Definition A.2 (Convex order).

Let $X$ and $Y$ be two random variables such that

$\displaystyle\mathbb{E}(\phi(X))\leq\mathbb{E}(\phi(Y))\quad\text{for all
convex functions }\phi:\mathbb{R}\rightarrow\mathbb{R},$

provided the expectations exist. Then $X$ is said to be smaller than $Y$ in
the convex order denoted as $X\leqslant_{\text{cx}}Y$.

The next results follows from the definition of convex ordering, see Chapter 3
of [SS].

######  Lemma A.3 (Theorem 3.A.24 in [SS]).

Let $X$ be a random variable with mean $\mathbb{E}X$. Denote the left (right)
endpoint of the support of $X$ by $l_{X}\left[u_{X}\right]$. Let $Z$ be a
random variable such that
$\mathbf{P}\left\\{Z=l_{X}\right\\}=\left(u_{X}-\mathbb{E}X\right)/\left(u_{X}-l_{X}\right)$
and $\mathbf{P}\left\\{Z=u_{X}\right\\}=(\mathbb{E}X-$
$\left.l_{X}\right)/\left(u_{X}-l_{X}\right)$. Then

$\mathbb{E}X\leqslant_{\text{cx}}X\leqslant_{\text{cx}}Z$

where $\mathbb{E}X$ denotes a random variable that takes on the value
$\mathbb{E}X$ with probability 1 (the left handside just restates Jensen’s
inequality).

######  Lemma A.4 (Theorem 3.A.44 in [SS]).

Let $X$ and $Y$ be two random variables with equal means, density functions
$f$ and $g$, distribution functions $F$ and $G$, and survival functions
$\bar{F}$ and $\bar{G}$, respectively. Denote by $S^{-}(a)$ the number of sign
changes for function $a$. Then $X\leq_{\mathrm{cx}}Y$ if any of the following
conditions hold:

$\displaystyle S^{-}(g-f)=2\mbox{ and the sign sequence is }+,-,+;$
$\displaystyle S^{-}(\bar{F}-\bar{G})=1\mbox{ and the sign sequence is }+,-;$
$\displaystyle S^{-}(G-F)=1\mbox{ and the sign sequence is }+,-.$

######  Proposition A.5 (Theorem 3.A.12 d) in [SS]).

Let $X_{1},X_{2},\ldots,X_{m}$ be a set of independent random variables and
let $Y_{1},Y_{2},\ldots,Y_{m}$ be another set of independent random variables.
If $X_{i}\leqslant_{\text{cx}}Y_{i}$ for $i=1,2,\ldots,m$, then

$\sum_{j=1}^{m}X_{j}\leqslant_{\text{cx}}\sum_{j=1}^{m}Y_{j}.$

That is, the convex order is closed under convolutions.

######  Lemma A.6 (Example 3.A.48 in [SS]).

Let $X$ and $Y$ be Bernoulli random variables with parameters $p$ and $q$,
respectively, with $0<p\leq q\leq 1$. Then

$\frac{X}{p}\geq_{\mbox{cx}}\frac{Y}{q}.$

######  Lemma A.7 (Inverse moment for the Binomial distribution).

Let $B_{1},\ldots,B_{k}\sim\mathbf{Bin}(1,q)$. Then
$\mathbb{E}[1/(1+\sum_{i=1}^{k}B_{i})]\leq 1/((k+1)q)$.

###### Proof.

See e.g. Benjamini et al., (2006). ∎

##  Appendix B Complements to Section 4.1

In the context of Gaussian one-sided testing described in Section 4.1, let
$\widehat{m}_{0}(p_{1},\ldots,p_{m})=\frac{1}{\nu}\left(1+\sum_{i=1}^{m}g(p_{i})\right)\in\mathcal{F}_{0}$.
Define $X_{0}\sim g(p_{i})$ for $i\in\mathcal{H}_{0}$ and $X_{1}\sim g(p_{i})$
where $i\in\mathcal{H}_{1}$. Then we have

$\displaystyle\textnormal{Bias}(\widehat{m}_{0})$
$\displaystyle=\frac{1}{\nu}\mathbb{E}\left(1+\sum_{i\in\mathcal{H}_{1}}g(p_{i})\right)=(1+(m-m_{0})\cdot\mathbb{E}X_{1})/\nu,$
$\displaystyle\textnormal{Var}(\widehat{m}_{0})$
$\displaystyle=\frac{1}{\nu^{2}}\textnormal{Var}\left(\sum_{i\in\mathcal{H}_{0}}g(p_{i})+\sum_{i\in\mathcal{H}_{1}}g(p_{i})\right)=(m_{0}\cdot\textnormal{Var}(X_{0})+(m-m_{0})\cdot\textnormal{Var}(X_{1}))/\nu^{2},\qquad\text{with}$
$\displaystyle\textnormal{Var}(X_{0})$
$\displaystyle=\int_{0}^{1}g(u)^{2}du-\left[\int_{0}^{1}g(u)du\right]^{2},$
$\displaystyle\textnormal{Var}(X_{1})$
$\displaystyle=\int_{0}^{1}g(u)^{2}f_{1}(u)du-\left[\int_{0}^{1}g(u)f_{1}(u)du\right]^{2}.$

where $f_{1}(t)=\exp\left(-\mu\cdot\Phi^{-1}(t)-\mu^{2}/2\right)$ denotes the
density of the $p$-values under the alternative.

##  Appendix C Complements to Section 4.2

Here we present some numerical results, comparing the performance of
$\widehat{m}_{0}^{\text{PC,new}}$ (see (10)) and
$\widehat{m}_{0}^{\text{PC,ZZD}}$ (see (17)) for $m=500$, where the correction
factors $C(500)=1.011709$ and $s(500)=98$ are taken from Table S1 in Zeisel et
al., (2011).

We first analyze the two estimators on simulated data in a one-sided Gaussian
testing setting where we observe realizations of independent rv’s
$X_{1},\dots,X_{m_{0}}\sim N(0,1)$ and $X_{m_{0}+1},\dots,X_{500}\sim
N(1.5,1)$ for $1000$ Monte-Carlo simulation runs and a varying range of
$m_{0}=50,100,\dots,450$. We obtain $500$ $p$-values by testing the null
hypotheses $H_{0,i}:\mu=0$ vs. the alternatives $H_{1,i}:\mu>0$ simultaneously
for all $i\in\\{1,\dots,500\\}$ and calculate
$\widehat{m}_{0}^{\text{PC,new}}$ and $\widehat{m}_{0}^{\text{PC,ZZD}}$ as
well as the number of rejections obtained from the plug-in BH procedure in (2)
with $\alpha=0.05$.

|
---|---
Figure 5: Box plots for the difference between
$\widehat{m}_{0}^{\text{PC,new}}$ and $\widehat{m}_{0}^{\text{PC,ZZD}}$ for
point estimation (left) and rejection numbers for the plug-in BH procedures
(right) against a range of true $m_{0}=50,100,\dots,450$.

Figure 5 shows that over a wide range of true $m_{0}$ values,
$\widehat{m}_{0}^{\text{PC,new}}$ and $\widehat{m}_{0}^{\text{PC,ZZD}}$ yield
comparable results both w.r.t. the point estimates and for the number of
rejections. In fact, $\widehat{m}_{0}^{\text{PC,new}}$ appears to be slightly
more efficient than $\widehat{m}_{0}^{\text{PC,ZZD}}$.

Another comparison can be obtained when we assume that the signal under the
alternative is strong and that most hypotheses are nulls. In this case we have
$2\sum_{i=1}^{m}p_{i}\approx 2\sum_{i\in\mathcal{H}_{0}}p_{i}=:S$ so that we
can use the Central Limit Theorem to quantify the probability that
$\widehat{m}_{0}^{\text{PC,new}}$ is more conservative than
$\widehat{m}_{0}^{\text{PC,ZZD}}$

$\displaystyle\mathbf{P}(\widehat{m}_{0}^{\text{PC,new}}>\widehat{m}_{0}^{\text{PC,ZZD}})$
$\displaystyle=\mathbf{P}(S>m\cdot
C(m)-2)\approx\overline{\Phi}\left(\sqrt{\frac{3}{m_{0}}}\cdot(m\cdot
C(m)-(m_{0}+2))\right).$

Figure 6 shows that this probability, for various values of the true $m_{0}$,
is quite small and even under the complete null ($m_{0}=500$) it is bounded by
$1/3$.

Figure 6: Approximate probabilities
$\mathbf{P}(\widehat{m}_{0}^{\text{PC,new}}>\widehat{m}_{0}^{\text{PC,ZZD}})$
for various values of true $m_{0}$, with fixed $m=500$.

##  Appendix D Additional Figures for simulated data of Section 4

We provide additional results on simulated data in the Gaussian one-sided
testing setting described in Section 4.1, with $m=10000$ and $\mu=1.5$. Figure
7 displays estimation results for $\pi_{0}$ over $1000$ Monte-Carlo
replications. They are in line with the analytical comparisons of the MSE
provided in Figure 1. Alongside, we also provide results on power, defined as
the ratio of the number of true discoveries to the number of alternatives, for
the corresponding plug-in BH (abbreviated in ABH for adaptive BH) procedures
using each of the estimators, the raw BH and oracle plug-in BH (using the true
$m_{0}$). The procedures are run for a fixed level $\alpha=0.05$. The power
enhancement among the different plug-in estimators’ is not striking except
perhaps for very small values of $\pi_{0}$ were where we recover the same
performance ranking as in Figure 1. For larger values of $\pi_{0}$, the
differences in power is not perceptible anymore, every procedure behaves
poorly as there is less and less signal.

|
---|---
Figure 7: Estimation results (left panel) for
$\widehat{m}_{0}^{\text{Storey}}$, $\widehat{m}_{0}^{\text{PC,new}}$,
$\widehat{m}_{0}^{\text{Poly}}(1,1/2)$, and
$\widehat{m}_{0}^{\text{Poly}}(2,1/2)$, and power results (right panel) for
the associated plug-in BH procedures on simulated data.

##  Appendix E Upper and lower bounds for the inverse moment of the uniform
sum distribution

The Pounds and Cheng estimator is closely related to the sum of independent
uniform random variables. This distribution plays a role in various contexts
and is also known as the Irwin-Hall distribution (for more details, see
Johnson et al., (1970)). As an auxiliary result, we give lower and upper
bounds for the inverse moment of this distribution.

######  Lemma E.1 (Inverse moments for Erlang distributions).

Let $E_{1},\ldots,E_{k}\sim\mathcal{E}(1)$ be independent exponentially
distributed random variables. Then $\mathbb{E}[1/\sum_{i=1}^{k}E_{i}]\leq
1/(k-1).$

###### Proof.

Since $X=\sum_{i=1}^{k}E_{i}$ is Gamma-distributed with shape $\alpha=k$ and
inverse scale parameter $\beta=1$ then $1/X$ is Inverse-gamma distributed with
mean $\beta/(\alpha-1)$, see Gelman et al., (2013). ∎

######  Proposition E.2.

(Inverse moment for sums of uniforms) For $k\geq 2$ let
$U_{1},U_{2},\ldots,U_{k}\sim\mathcal{U}[0,1]$ iid. Then we have

$\displaystyle\frac{2}{k}$
$\displaystyle\leq\mathbb{E}\left(\frac{1}{\sum_{i=1}^{k}U_{i}}\right)\leq\frac{2}{k-1}$
(27)

###### Proof.

Let $E_{1},E_{2},\ldots,E_{k}\sim\mathcal{E}(1)$ iid. From Theorems 3.A.24 and
3.A.46 in [SS] we have for $i=1,\ldots,,k$

$\displaystyle 1$
$\displaystyle\leqslant_{\text{cx}}2U_{i}\leqslant_{\text{cx}}E_{i}$ and since
the convex ordering is preserved under convolutions (see [SS] Theorem 3.A.12.)
we obtain $\displaystyle k$
$\displaystyle\leqslant_{\text{cx}}\sum_{i=1}^{k}2U_{i}\leqslant_{\text{cx}}\sum_{i=1}^{k}E_{i}.$
Together with the convexity of the mapping $x\mapsto 1/x$ on $(0,1)$ this
yields $\displaystyle\frac{1}{k}$
$\displaystyle\leq\mathbb{E}\left(\frac{1}{\sum_{i=1}^{k}2U_{i}}\right)\leq\mathbb{E}\left(\frac{1}{\sum_{i=1}^{k}E_{i}}\right)\leq\frac{1}{k-1},$

where the last inequality follows from Lemma E.1. ∎

## References

  * Benjamini, (2010) Benjamini, Y. (2010).  Simultaneous and selective inference: Current successes and future challenges.  Biometrical Journal, 52(6):708–721.

  * Benjamini and Hochberg, (1995) Benjamini, Y. and Hochberg, Y. (1995).  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  Journal of the Royal statistical society: series B (Methodological), 57(1):289–300.
  * Benjamini 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 multiple testing under dependency.  Annals of statistics, pages 1165–1188.
  * Berry and Armitage, (1995) Berry, G. and Armitage, P. (1995).  Mid-p confidence intervals: a brief review.  Journal of the Royal Statistical Society Series D: The Statistician, 44(4):417–423.
  * Biswas and Chattopadhyay, (2020) Biswas, A. and Chattopadhyay, G. (2020).  Modified estimator for the proportion of true null hypotheses under discrete setup with proven fdr control by the adaptive benjamini-hochberg procedure.  arXiv preprint arXiv:2009.03803.
  * Blanchard and Roquain, (2009) Blanchard, G. and Roquain, É. (2009).  Adaptive false discovery rate control under independence and dependence.  Journal of Machine Learning Research, 10(12).
  * Chavant et al., (2011) Chavant, F., Favrelière, S., Lafay-Chebassier, C., Plazanet, C., and Pérault-Pochat, M.-C. (2011).  Memory disorders associated with consumption of drugs: updating through a case/noncase study in the french pharmacovigilance database.  British journal of clinical pharmacology, 72(6):898–904.
  * Chen and Doerge, (2015) Chen, X. and Doerge, R. (2015).  fdrdiscretenull: False discovery rate procedure under discrete null distributions.  R package version, 1.
  * Chen et al., (2018) Chen, X., Doerge, R. W., and Heyse, J. F. (2018).  Multiple testing with discrete data: proportion of true null hypotheses and two adaptive fdr procedures.  Biometrical Journal, 60(4):761–779.
  * Dickhaus et al., (2012) Dickhaus, T., Straßburger, K., Schunk, D., Morcillo-Suarez, C., Illig, T., and Navarro, A. (2012).  How to analyze many contingency tables simultaneously in genetic association studies.  Statistical applications in genetics and molecular biology, 11(4).
  * Ditzhaus and Janssen, (2019) Ditzhaus, M. and Janssen, A. (2019).  Variability and stability of the false discovery proportion.  Electronic Journal of Statistics, 13(1):882–910.
  * Döhler et al., (2018) Döhler, S., Durand, G., and Roquain, E. (2018).  New fdr bounds for discrete and heterogeneous tests.  Electronic Journal of Statistics, 12(1):1867–1900.
  * Durand, (2019) Durand, G. (2019).  Adaptive p-value weighting with power optimality.
  * Embrechts and Hofert, (2013) Embrechts, P. and Hofert, M. (2013).  A note on generalized inverses.  Mathematical Methods of Operations Research, 77(3):423–432.
  * Gelman et al., (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013).  Bayesian data analysis.  Chapman and Hall/CRC Texts in Statistical Science Series. CRC, Boca Raton, Florida, third edition.
  * Habiger, (2015) Habiger, J. D. (2015).  Multiple test functions and adjusted p-values for test statistics with discrete distributions.  Journal of Statistical Planning and Inference, 167:1–13.
  * Habiger and Pena, (2011) Habiger, J. D. and Pena, E. A. (2011).  Randomised p-values and nonparametric procedures in multiple testing.  Journal of nonparametric statistics, 23(3):583–604.
  * Heesen and Janssen, (2016) Heesen, P. and Janssen, A. (2016).  Dynamic adaptive multiple tests with finite sample fdr control.  Journal of Statistical Planning and Inference, 168:38–51.
  * Heller and Gur, (2011) Heller, R. and Gur, H. (2011).  False discovery rate controlling procedures for discrete tests.  arXiv preprint arXiv:1112.4627, pages 00092–17.
  * Heller et al., (2012) Heller, R., Gur, H., and Yaacoby, S. (2012).  discretemtp: Multiple testing procedures for discrete test statistics.  R package version 0.1-2.
  * Ignatiadis et al., (2016) Ignatiadis, N., Klaus, B., Zaugg, J. B., and Huber, W. (2016).  Data-driven hypothesis weighting increases detection power in genome-scale multiple testing.  Nature methods, 13(7):577–580.
  * Johnson et al., (1970) Johnson, N. L., Kotz, S., and Johnson, N. L. (1970).  Continuous univariate distributions.
  * Karp et al., (2017) Karp, N. A., Mason, J., Beaudet, A. L., Benjamini, Y., Bower, L., Braun, R. E., Brown, S. D., Chesler, E. J., Dickinson, M. E., Flenniken, A. M., et al. (2017).  Prevalence of sexual dimorphism in mammalian phenotypic traits.  Nature communications, 8(1):1–12.
  * MacDonald et al., (2019) MacDonald, P. W., Liang, K., and Janssen, A. (2019).  Dynamic adaptive procedures that control the false discovery rate.
  * Muñoz-Fuentes et al., (2018) Muñoz-Fuentes, V., Cacheiro, P., Meehan, T. F., Aguilar-Pimentel, J. A., Brown, S. D., Flenniken, A. M., Flicek, P., Galli, A., Mashhadi, H. H., Hrabě de Angelis, M., et al. (2018).  The international mouse phenotyping consortium (impc): a functional catalogue of the mammalian genome that informs conservation.  Conservation genetics, 19:995–1005.
  * Neuvial, (2013) Neuvial, P. (2013).  Asymptotic results on adaptive false discovery rate controlling procedures based on kernel estimators.  Journal of Machine Learning Research, 14(44):1423–1459.
  * Pounds and Cheng, (2006) Pounds, S. and Cheng, C. (2006).  Robust estimation of the false discovery rate.  Bioinformatics, 22(16):1979–1987.
  * Rubin-Delanchy et al., (2019) Rubin-Delanchy, P., Heard, N. A., and Lawson, D. J. (2019).  Meta-analysis of mid-p-values: some new results based on the convex order.  Journal of the American Statistical Association, 114(527):1105–1112.
  * Sarkar, (2008) Sarkar, S. K. (2008).  On methods controlling the false discovery rate.  Sankhyā: The Indian Journal of Statistics, Series A (2008-), pages 135–168.
  * Shaked and Shantikumar, (2007) Shaked, M. and Shantikumar, J. G. (2007).  Stochastic orders.  Springer.
  * Storey, (2002) Storey, J. D. (2002).  A direct approach to false discovery rates.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 64(3):479–498.
  * Storey et al., (2004) Storey, J. D., Taylor, J. E., and Siegmund, D. (2004).  Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66(1):187–205.
  * Zeisel et al., (2011) Zeisel, A., Zuk, O., and Domany, E. (2011).  Fdr control with adaptive procedures and fdr monotonicity.  The Annals of applied statistics, pages 943–968.