I. INTRODUCTION

Consider the problem of estimating the proportion ${\pi }_0$ of true null hypotheses in a collection of m tests, given the observed p-values $p_1,\cdots ,p_m$ for the collection of m tests. This problem has attracted a lot of attentions in statistical literatures, attributing to its important role in dealing with multiple testing procedures, since the nominal paper Storey (2002). It has naturally arisen in assessing or controlling an overall false rejection rate, i.e., the false discovery rate (FDR) proposed by Benjamini and Hochberg (1995) in multiple null hypothesis testing problems and multiple comparisons. Benjamini and Hochberg proved that if $p_1,\cdots ,p_m$ are independent with continuous distributions, the popular Sime's multiple testing procedure (Sime 1986), where a null hypothesis is rejected whenever the observed p-value is less than $\alpha$, results in an FDR controlled by ${\pi }_0\alpha$. Therefore, a reliable estimate is essential to control the FDR in the Sime's multiple testing procedure and improve its testing power. On the other hand, the proportion ${\pi }_0$ is a quantity of interest in its own right. For example, researchers may wish to estimate the proportion of genes that are not differentially expressed in DNA microarray experiments (Parker et. al. 1988).

Store (2002) proposed to approach the estimating problem from a Bayesian point of view by treating $p_1,\cdots ,p_m$ as a random sample of size $m$ from a mixture distribution of a uniform distribution and a non-uniform distribution with mixing proportion ${\pi }_0$.

By Storey's approach, the observed $p$-values are considered to be a random sample of size $m$ from $f(x|{\pi }_0,h)$. As a consequence, the mixing proportion ${\pi }_0$ in the mixture model (1) represents the proportion of true-null $p$-value subpopulation and an estimate for ${\pi }_0$ based on a random sample from (1) thus yields an estimate for the proportion of true null hypotheses among the $m$ null hypotheses in the multiple testing problem.

The first well-studied estimator for ${\pi }_0$ was proposed by Storey (2002) as follows:

where $\lambda$ is an appropriately chosen value close to 1 and $W(\lambda )$ is the total number of the p-values greater than $\lambda$. This estimator is motivated by the fact that false-null p-values are typically small so that with $\lambda$ close to 1, $1-\lambda$ is close to the expected proportion of true-null p-values falling into $(\lambda ,1]$. Since the proportion of the p-values falling into $(\lambda ,1]$ is a mixture of $1-\lambda$ and ${\int }_{\lambda }^{1}h(x)dx$ with mixing proportions ${\pi }_0$, ${\hat{\pi }}_0^s$ is obtained mathematically by simply setting the proportion of false-null p-values falling into $(\lambda ,1]$ to be zero, i.e., ${\int }_{\lambda }^{1}h(x)dx=0$. As a consequence, ${\hat{\pi }}_0^s$ tends to overestimate ${\pi }_0$, as ${\int }_{\lambda }^{1}h(x)dx$ is typically positive. The biasedness can be too great to be acceptable when ${\pi }_0$ belongs to the lower part of [0, 1]. Storey (2002) proposed a bootstrap procedure to choose $\lambda$ that minimizes an upper bound of the mean square error of the resulting positive FDR estimator $\hat{Q}(\lambda )=\alpha {\hat{\pi }}_0^s(\lambda )/R(\alpha )$ where $R(\alpha )=\mathrm{\#}\{p_i\le \alpha \}$. See more discussion and some other work on how to select $\lambda$, see Storey, Taylor and Siegmund (2004), and Nettleton, Hwang, Caldo and Wise (2006).

Since the publication of Storey (2002), the mixture model approach described above has been adopted widely and many other estimators have been proposed. They include: Langaas, Lindqvist and Ferkingstad (2005)'s by a histogram approach, Wu, Guan and Zhao (2006)'s polynomial-type estimator, Jiang and Doerge (2008)'s estimator averaging over different $\lambda$ -values with Storey (2002)'s estimator, Zhao, Wu, Zhang and Chen (2012)'s estimator in exponential mixture model, Cheng, Gao and Tong (2015)'s estimator with reduction of estimating bias and standard deviation. For more references and discussion, see Cheng, Gao and Tong (2015), and Tong, Feng, Hilton and Zhao (2013) that deals with possibly dependent p-values. Nevertheless, when these estimators are intended to improve some aspects of Storey's estimator, the biasedness remains significant even with m as large as 2,000, especially when ${\pi }_0$ is not close to 1. Motivated by the histogram approach (see Mosig et al. 2001 and Nettleton et al. 2006), a new estimator is proposed in this paper via a likelihood approach with h being approximated by a modified histogram pdf. Akaike information criterion is used to select the number of categories in histogram construction. Simulation study demonstrates that the new estimator significantly improves popularly cited existing methods such as Storey (2002), Langaas, Lindqvist and Ferkingstad Langaas (2005), Jiang and Doerge (2008) and Nettleton et al.(2006).

The paper is organized as follows. The new method along with discussion of computing issues is presented in Section 2. Simulation results are reported in Section 3. A real-life example is given in Section 4.

II. METHOD

Let $p_1,\dots ,p_m$ be a random sample of size $m$ from the pdf $f(x|{\pi }_0,h)$ defined in (1). We propose ${\pi }_0$ to be estimated by the MLE when $h$ is subject to a histogram-type approximation. We shall proceed with a study on the identification problem in the mixture model.

a) Identification

Let $\mathrm{\Theta }$ be the parameter space consisting of all $({\pi }_0,h)$ 's with ${\pi }_0\in (0,1)$ and $h$ being a probability density function on $[0,1]$ and continuous in a neighborhood of 1 with $h(1-)=0$, where

For convenience and confirmation, when a pdf is used as a parameter throughout the paper, we should refer to the probability distribution, say the CDF specified by the pdf.

Lemma 1 The distribution $f(x|{\pi }_0,h)$ in the model (1) is identifiable in $({\pi }_0,h)\in \mathrm{\Theta }$.

Proof. Let two parameter sets $({\pi }_1,h_1)$ and $({\pi }_2,h_2)$ in $\mathrm{\Theta }$ define the same distribution, i.e., for any $x\in [0,1]$,

where $H_i$ is the CDF of $h_i,i=1,2$. Since $h_1$ and $h_2$ are continuous in a common neighborhood of 1, (2) implies that for any x in a neighborhood of 1,

So ${\pi }_1={\pi }_2$ because $h_1(1-)=h_2(1-)=0$. This in turn implies $H_1=H_2$ in (2) because and . The lemma is proved.

Remarks.

(a) First of all, note the simple fact that $h(1-)=0$ implies that $h$ is non-uniform. Of course, identification of the parameter ${\pi }_0$ of interest requires that $h$ is non-uniform.

(b) When ${\pi }_0=1$, $f$ is un-identifiable in $h$.

(c) The assumption that $h$ is continuous in a neighborhood of 1 with $h(1-)=0$ is critical. Without it, the lemma can fail. For a counterexample, let ${\pi }_1=1/3$ and $h_1(x)=1/4+(3/2)x$, and ${\pi }_2=1/2$ and $h_2(x)=2x$. We have

and for any $x$ in $[0,1]$

and

That is, ${\pi }_1+(1-{\pi }_1)h_1(x)={\pi }_2+(1-{\pi }_2)h_2(x)$ for all $x$ in $[0,1]$, but ${\pi }_1\ne {\pi }_2$ and $h_1\ne h_2$.

b) A

Motivated by the histogram approach (see Mosig et al. (2001), a histogram approximation to the alternative pdf $h(x)$ is proposed as follows. Let k > 2 be an integer. Define

Notes

where $0\le q_1\le 1,\cdots ,0\le q_{k-1}\le 1$ satisfying $q_1+\dots +q_{k-1}+q_{k-1}/2=1$. The linear modification over the right most subinterval $(1-1/k,1]$ warrants $\hat{h}\in \mathrm{\Theta }$ so that $f(x|{\pi }_0,\tilde{h})$ is identifiable in ${\pi }_0$ and $q=(q_1,\cdots ,q_{k-2})$. As Storey(2002) and Tong et al. (2013) remarked, the false-null p-values that follow the distribution h are typically small. In other words it is often the case in application that the alternative pdf h is highly skewed to the right. Therefore, with the linear modification on the rightmost subinterval $(1-1/k,1]$, the effect on the accuracy of the estimate for Storey (2002), Jiang and Doerge (2008) and some other methods assume that h is zero in a neighborhood of 1.

c) The new estimator with fixed k

Suppose that k > 2 is specific for now (a selection procedure will be described later). When h is approximated by or restricted to $\tilde{h}$, with a random sample $p_1,\cdots ,p_m$, the log-likelihood of the parameter ${\pi }_0$ of interest and the new nuisance parameter q becomes

where ${\xi }_{ij}$ is the indicator whether $p_i$ falls into the j-th category of the histogram with k bins or not, i.e.,

for $i=1,\cdots ,m,j=1,\cdots ,k.$

Let $({\hat{\pi }}_0(k),\hat{q})$ be the MLE of $({\pi }_0,q)$ with the log-likelihood function (3), subject to and , $q_k=q_{k-1}/2$ such that $\sum _{j=1}^k{q}_j=1$. This defines the new estimator ${\hat{\pi }}_0(k)$ for ${\pi }_0$ with specified k. A direct application of the standard theory of MLE yields the following theorem and hence a theory of consistency of the new estimator.

Theorem 1 For any $({\pi }_0,h)\in \mathrm{\Theta }$ and , let

Notes for $j=1,\cdots ,k-2$. Put $q=(q_1,\cdots ,q_{k-2})$ and $\theta =({\pi }_0,q)$. Then $\sqrt{m}({\hat{\pi }}_0(k)-{\pi }_0,\hat{q}-q)$ converges to $N(0,\mathrm{\Sigma })$ in distribution, as $m\to \mathrm{\infty }$, where

d) Computation of ${\hat{\pi }}_0(k)$

Maximizing the nonlinear log-likelihood function (3) can be complicating. However, the EM algorithm can be used to obtain an approximation to the MLE ${\hat{\pi }}_0(k)$ easily. To do so, introduce a latent Bernoulli variable $w$ indicating the component-ship of the $p$-value in the finite mixture distribution. That is, let $w$ be a binary random variable with $P(w=1)={\pi }_0$ and $P(w=0)=1-{\pi }_0$, and let a random variable $p$ follow the distribution as follows. Given $w=1$, $p$ follows the uniform distribution on $[0,1]$ and given $w=0$, $p\sim \tilde{h}$. It is clear that $p$ follows the mixture distribution $f(x|{\pi }_0,\tilde{h})$. Thus $\mathcal{P}=\{p_1,\cdots ,p_m\}$ can be viewed as the incomplete data of a random sample $(p_1,w_1),\cdots ,(p_m,w_m)$ from $(p,w)$ with missing values $w_1,\cdots ,w_m$.

Note that with the complete data $(p_1,w_1),\cdots ,(p_m,w_m)$, the likelihood function of $({\pi }_0,q)$ is

where ${\xi }_{ij}$ 's are defined in (4), and the log-likelihood function is thus

where ${\xi }_{ij}^{\ast }=(1-w_i){\xi }_{ij},w_{.}=\sum _{i=1}^m{w}_i$, and ${\xi }_{.j}^{\ast }=\sum _{i=1}^m{\xi }_{ij}^{\ast }$.

The EM algorithm can be easily implemented as follows. Let $({\pi }_0^{(s)},q^{(s)})$ be the current approximations to the MLE $({\hat{\pi }}_0(k),\hat{q})$ with the log-likelihood function l given in (3). The next approximation $({\pi }_0^{(s+1)},q^{(s+1)})$ is given by the EM algorithm in two steps, the so-called E-step and M-step.

E-Step: Compute the conditional expectation of the log-likelihood function

Note

Notes

We have

M-Step: In the M-step, $Q({\pi }_0,q)$ is maximized to yield the next approximation ${\pi }_0^{(s+1)}$ and $q^{(s+1)}$. Setting $\partial Q/\partial {\pi }_0=0$, we immediately have

By Lemma 2 in Appendix,

The formulas in the EM algorithm above are similar to those developed by Oluyemi (2016) and Oluyemi and Chen (2016) where there are some algebraic errors.

where