Estimating the Proportion of True Null Hypotheses for Multiple Comparisons

Abstract

Whole genome microarray investigations (e.g. differential expression, differential methylation, ChlP-Chip) provide opportunities to test millions of features in a genome. Traditional multiple comparison procedures such as familywise error rate (FWER) controlling procedures are too conservative. Although false discovery rate (FDR) procedures have been suggested as having greater power, the control itself is not exact and depends on the proportion of true null hypotheses. Because this proportion is unknown, it has to be accurately (small bias, small variance) estimated, preferably using a simple calculation that can be made accessible to the general scientific community. We propose an easy-to-implement method and make the R code available, for estimating the proportion of true null hypotheses. This estimate has relatively small bias and small variance as demonstrated by (simulated and real data) comparing it with four existing procedures. Although presented here in the context of microarrays, this estimate is applicable for many multiple comparison situations.

Introduction

Genomic technologies are producing vast amounts of biological data that are the basis for investigations that require repetitive testing of the same hypothesis. Because the number of tests performed (e.g. differential expression) is so large, sometimes the multiple comparison procedures that control the familywise error rate are too strict for biological applications (e.g. differential methylation). In fact, many biologists would rather experience several more false positives (i.e. type I errors; false rejections of the null hypothesis) than lose important information. In an attempt to address the multiple comparison issue Benjamini and Hochberg (1995) introduced an error rate measure called False Discovery Rate (FDR). Specifically, a family of m hypothesis tests is considered, of which m₀ are true. The proportion of erroneously rejected null hypotheses among all the rejected null hypotheses can be captured by the random variable Q = V/R, where R is the number of rejected hypotheses and V is the number of false rejections (type I errors). Benjamini and Hochberg (1995) formally define the FDR to be the expected proportion of falsely rejected hypotheses among all the rejections,

$$\text{FDR=E}\left(Q\right)=E\left(V/R\right),$$

(1)

where Q = 0 when R = 0 (no rejections).

If we let p₍₁₎ ≤ p₍₂₎ ≤ … ≤ P_(m) be the ordered p-values and H_(i) be the null hypothesis corresponding to p_(i), then in Benjamini and Hochberg's (BH) FDR controlling procedure (Benjamini and Hochberg, 1995), K is considered to be the largest k such that p_k ≤ (k/m)α, where α is the pre-chosen FDR significance level. If K exists, all null hypotheses H_(i), i = 1, …, K are rejected. If no such K exists, then no hypotheses are rejected. The BH FDR controlling procedure controls the FDR at exactly the level (m₀/m)α≤ α, and hence conservatively at a, for independent test statistics and for any configuration of false null hypotheses (Benjamini and Yekutieli, 2001; Storey, et al. 2004). In 2000 Benjamini and Hochberg proposed an adaptive procedure which provides more power than the original FDR controlling procedure by comparing each p_(k) with $\left(k/{\widehat{m}}_0\right)\alpha$ where ${\widehat{m}}_0$ is an estimate of m₀. If the estimated value of m₀ is such that ${\widehat{m}}_0\ge m_0$ with probability one, then the adaptive BH FDR controlling procedure will lead to $FDR=\frac{m_0}{m}\left(\frac{m}{{\widehat{m}}_0}\alpha \right)=\frac{m}{{\widehat{m}}_0}\alpha \le \alpha$ . Because the accuracy and variation of the estimate of m₀, or π₀ = m₀/m, directly affects the performance of the adaptive FDR controlling procedure our focus is on the estimation and effect of π₀.

We propose a simple and easy-to-implement method for estimating the proportion of true null hypotheses. The performance of this estimate is compared to existing methods via simulated and real data. Specifically, Benjamini and Hochberg (2000) estimated the number of true hypotheses from the observed p-values using the Lowest SLope (LSL) estimator. Their approach was based on a modification of the graphical method of Schweder and Spjotvoll (1982). Alternatively, Storey (2002) proposed an estimate of π₀ by assuming the p-values corresponding to true null hypotheses are uniformly distributed on the interval (0,1) and selecting a reasonable tuning parameter 0 ≤ λ < 1. Finally, Langaas, et al. (2005) derived estimators based on nonparametric maximum likelihood estimation of the p-value density, under the restriction of decreasing and convex decreasing densities. Although Benjamini and Hochberg's original and adaptive FDR controlling procedure are developed for independent statistics these procedures can also be applied to some dependence structures (Benjamini and Yekutieli, 2001). Simulations have also demonstrated that they can be used for situations where there is a weak correlation structure among the genes (Storey et al. 2004). However, because of the small number of biological replicates used in most microarray experiments, it is very difficult to measure the correlation structure among a set or family of genes. Reiner et al. (2003) proposed a procedure for the general case, but it is conservative when compared to the adaptive FDR controlling procedures.

Methods

Storey's Approach

Our approach is motivated by the work of Storey (2002), where the proportion of true null hypotheses, % is estimated by

$${\widehat{\pi }}_0(\lambda )=\frac{W(\lambda )}{(1-\lambda )m},$$

(2)

where $W\left(\lambda \right)=\#\left\{p_i:p_i>\lambda \right\},\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}0\le \lambda <1$ is a tuning parameter. The rationale for this estimate is that p-values corresponding to true null hypotheses are uniformly distributed on the interval (0,1), of which most should be close to 1. Thus, for a reasonable A, there are about m₀(1 – λ) such p-values in the interval (λ, 1] such that W(λ) ≈ m₀(1 – λ). Black (2004) pointed out that Equation (2) is an unbiased estimate of π₀ for all values of λ if all the null hypotheses are true and the p-values have a uniform distribution on the interval (0,1). However, there is an upward bias when the p-values come from both true null and true alternative hypotheses. As it turns out, choosing the tuning parameter λ in Equation (2) is very important since there is a bias-variance trade-off. When λ ↠ 0, the variance of ${\widehat{\pi }}_0\left(\lambda \right)$ becomes smaller and the bias of this estimate increases. When λ ↠ 1, the bias of ${\widehat{\pi }}_0\left(\lambda \right)$ becomes smaller, and the variance of this estimate increases. To address this point, Storey et al. (2004) proposed a bootstrap method that automatically chooses λ when estimating ${\widehat{\pi }}_0\left(\lambda \right)$ .

Instead of choosing one specific λ, Storey and Tibshirani (2003) proposed an estimate of π₀ using $\underset{\lambda \to 1}{\lim }{\widehat{\pi }}_0\left(\lambda \right)$ so that the bias is small and there is a balance between both bias and variance. For this approach, ${\widehat{\pi }}_0\left(\lambda \right)$ is plotted over a range of λ = 0,0.05,0.10, …, 0.90 and then a natural cubic smoothing spline is fit to these data for the purpose of estimating the overall trend of ${\widehat{\pi }}_0\left(\lambda \right)$ as λ ↠ 1. In the QVALUE (http://www.faculty.washington.edu/~jstorey/) software, the predicted value of ${\widehat{\pi }}_0\left(\lambda \right)$ at λ = 0.90 is chosen as the estimate of π₀.

Average Estimate Approach

As mentioned previously, the estimate ${\widehat{\pi }}_0\left(\lambda \right)=\frac{W\left(\lambda \right)}{\left(1-\lambda \right)m}$ where 0 ≤ λ < 1, has a large bias and small variance when λ is small and a small bias and large variance when λ is big. Suppose for each λ _i , where 0 < λ₁ < λ₂ < … < λ _n < 1, we compute ${\widehat{\pi }}_0\left({\lambda }_i\right)$ as in Equation (2), then

$$E\left[{\widehat{\pi }}_0\left({\lambda }_i\right)\right]={\pi }_0+{\epsilon }_i,$$

where $E\left[{\epsilon }_i\right]={\delta }_{i,}{\delta }_i\ge {\delta }_{i+1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Var\left[{\epsilon }_i\right]={\sigma }_i^2,$ and ${\sigma }_i^2\le {\sigma }_{i+1}^2$ . Therefore, a natural choice is to consider the average of ${\widehat{\pi }}_0\left(\lambda \right)$ over the values of λ _i ,

$${\widehat{\pi }}_0=\frac{1}{n}\sum _{i=1}^{i=n}{\widehat{\pi }}_0\left({\lambda }_i\right).$$

The bias of $\hspace{0.17em}{\widehat{\pi }}_0,\hspace{0.17em}\hspace{0.17em}\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$n$}\right.{\sum }_{i=1}^{i=n}{\delta }_i$ , is smaller than δ₁ (the bias of the estimate of π₀ at λ = λ₁) and at the same time, ${\widehat{\pi }}_0$ has a smaller variance. Considering the average of ${\widehat{\pi }}_0\left(\lambda \right)$ over a range of λ to estimate π₀ reduces the problem to choosing the range of λ.

Define 0 = t₁ < t₂ < … < t_B < t_B+1 = 1 as equally spaced points in the interval [0,1] such that the interval [0,1] is divided into B small intervals with equal length 1/B. Specifically, t_i = (i–1)/B. For example, when B = 10, t₁ = 0, t₂ = 0.1, …, t₁₀ = 0.9. For each t_i (i = 1,…,B), ${\widehat{\pi }}_0\left(t_i\right)$ is an estimate of π₀ via Equation (2) with λ = t_i. The goal then becomes finding a subset of t_i's such that a new estimate of π₀ is obtained by taking the average of the corresponding values of ${\widehat{\pi }}_0\left(t_i\right)$ . Let NB_i denote the number of p-values which are greater than or equal to t_i, and let NS_i represent the number of p-values in the interval of [t_i, t_i+1). Therefore,

$$N{B}_i=\#\left\{p_k:p_k\ge t_i\right\},$$

(3)

$${\widehat{\pi }}_0\left(t_i\right)=\frac{N{B}_i}{\left(1-t_i\right)m}$$

(4)

$$N{S}_i=\#\left\{p_k:t_i\le p_k<t_{i+1}\right\},$$

(5)

where i = 1,…,B.

If the NB_i p-values come from the null distribution, then on average there are $\frac{N{B}_i}{B-i+1}$ p-values in each of the (B–i + 1) small intervals on [t_i, 1]. In other words, there are $\frac{N{B}_i}{B-i+1}$ p-values in each small interval [t_j, t_j+1) for i ≤ j ≤ B. Since the p-values corresponding to the true alternative hypotheses are smaller than those corresponding to the true null hypotheses, there are more p-values in the intervals [t_i, t_i+1) with small index i. For small i, NS_i is usually greater than $\frac{N{B}_i}{B-i+1}$ . Therefore, initiating from i = 1, we find the first i such that $N{S}_i\le \frac{N{B}_i}{B-i+1}$ . If such i exists, t_i can be considered as the change point and we assume all the p-values bigger than t_i come from the true null hypotheses. From this π₀ can be estimated by

$${\widehat{\pi }}_0\left(B\text{}\right)=\frac{1}{B-i+1}\sum _{j=i}^{j=B}{\widehat{\pi }}_0\left(t_j\right)$$

$$=\frac{1}{B-i+1}\sum _{j=i}^{j=B}\frac{N\hspace{0.17em}{B}_j}{\left(1-t_j\right)m}$$

(7)

where $i=\min \left\{i:N{S}_i\le \frac{N{B}_i}{B-i+1}\right\}$ . In order to find the range of λ, only a lower bound of λ is required. The large values of t_i are used so that it ensures the estimate of π₀ has small bias. This is equivalent to fitting a straight line with slope 0 in the right bottom part of a ${\widehat{\pi }}_0\left(t_i\right)$ versus t_i plot, such that the intercept provides the estimate of π₀. A simple modification of this approach is to estimate π₀ by taking the average of ${\widehat{\pi }}_0(t_j)$ from j = i–1 to B, that is,

$${\widehat{\pi }}_0\left(B\right)=\frac{1}{B-i+2}\sum _{j=i-1}^{j=B}\frac{N{B}_j}{\left(1-t_j\right)m}$$

(8)

where $i=\min \left\{i:N{S}_i\le \frac{N{B}_i}{B-i+1}\right\}$ . This ensures that the upward bias increases and the variance decreases, as ${\widehat{\pi }}_0\left(t_{i-1}\right)$ has smaller variance and bigger bias than ${\widehat{\pi }}_0\left(t_j\right)$ for j = i, …, B.

A remaining challenge is how to choose B. Specifically, how many λ's should be used in the interval [0,1]. Recall that a motivating factor of the proposed average estimate approach is to balance the bias and variance. The natural way to measure both the bias and variance is the mean-squared error, $E{\left[{\widehat{\pi }}_0\left(B\right)-{\pi }_0\right]}^2$ . Since the true value of π₀ is unknown and the theoretical result is intractable, we take a bootstrap approach in the following way:

1. 1.

For each B ∊ I, I = {5, 10, 20, 50,100}, compute ${\widehat{\pi }}_0\left(B\right)$ as in Equation (8).

2.

Form N bootstrap samples of the p-values, and compute the bootstrap estimates ${\widehat{\pi }}_0^{*b}\left(B\right)$ for b = 1, …, N and B ∊ {5, 10, 20, 50, 100}.

3.

For each B ∊ I, estimate its respective mean-squared error as

$$\widehat{MSE}\left(B\right)=\frac{1}{N}{\sum _{b=1}^n\left[{\widehat{{\pi }_0}}^{*b}\left(B\right)-{\overline{\pi }}_0\right]}^2,$$

Where

$${\overline{\pi }}_0={\text{average}}_{B^{\text{'}}\in I}\left\{{\widehat{\pi }}_0\left(B^{\text{'}}\right)\right\}$$

4.

Let $\widehat{B}=\arg {\min }_{B\in I}\widehat{MSE}\left(B\right)$ , then the estimate of π₀ is ${\pi }_0\hspace{0.17em}is\hspace{0.17em}{\widehat{\pi }}_0=\hspace{0.17em}{\widehat{\pi }}_0\left(\widehat{B}\right)$ .

Notice that in step three the value of π₀ is estimated by the average of the $\hspace{0.17em}{\widehat{\pi }}_0\left(B\right)$ over arrange of B.

Results

Simulation Studies

To investigate the performance of the proposed average estimate approach, a simulation study was performed. Taking m = 1,000 (i.e. 1,000 genes are tested for differential expression), let π₀ vary over a wide range, say π₀ = 0.50,0.60,…,0.90 which are reasonable for microarray experiments. Hypotheses, H₀: μ = 0 versus H_a. μ > 0, are tested for independent random variables Z_i(i = 1,…,m) from null distribution N(0,1) and alternative distribution N(2,1). Specifically, mπ₀ and m(1 –π₀) random variables have mean 0 and 2, respectively. For each test, the p-value is computed as p_i = P(Z > z_i), where Z is a random variable from a standard normal distribution N(0,1) and z_i is the observed value of Z_i. For each value of π₀, l = 1,000 data sets were simulated.

For the choice of B, B is either fixed (i.e. B = 5, 10, 20, 50, and 100) or chosen by the proposed bootstrap approach. For each of the l = 1,000 simulated data sets, when B is fixed, the estimate of π₀ is computed via Equation (8), that is, ${\widehat{\pi }}_0=\frac{1}{B-i+2}{\sum }_{j=i-1}^{j=B}{\widehat{\pi }}_0\left(t_j\right)$ where $i=\min \left\{i:N{S}_i\le \frac{N{B}_i}{B-i+1}\right\}$ . If such i does not exist, π₀ is estimated by the average of ${\widehat{\pi }}_0\left(t_{B-1}\right)$ and ${\widehat{\pi }}_0\left(t_B\right)$ . For the bootstrap approach to automatically choose B, the range of B is 5, 10,20,50, 100.

For completion the performance of the proposed average estimate approach is compared with several existing procedures:

1. 1.

Benjamini and Hochberg's lowest slope estimate (LSL) (Benjamini and Hochberg, 2000),

2.

Storey's bootstrap estimate (Storey_boot) (Storey et al. 2004),

3.

Storey and Tibshirani's smoother estimate (ST_smoother) (Storey and Tibshirani, 2003),

4.

Langass et al.'s nonparametric maximum likelihood estimate (convest) (Langaas et al. 2005).

For procedures 2 and 3, the QVALUE software (http://www.faculty.washington.edu/~jstorey/) was employed. For procedure 4, the R function ‘con-vest’ was downloaded from the R library ‘limma’ as part of the Bioconductor project (http://www.bioconductor.org).

Table 1 summarizes the simulation results. Bias and the standard deviation of the estimates are estimated by

$$\begin{array}{l}\widehat{Bias}=\frac{1}{l}{\displaystyle \sum _{i=1}^{i=l}}\left({\widehat{\pi }}_{0i}-{\pi }_0\right)\hfill \\ \widehat{Std}=\sqrt{\frac{1}{l-1}{{\displaystyle \sum _{i=1}^{i=l}}\left({\widehat{\pi }}_{0i}-\frac{1}{l}{\displaystyle \sum _{i=1}^{i=l}}{\widehat{\pi }}_{0i}\right)}^2,}\hfill \end{array}\hspace{0.17em}$$

Table 1

The estimate of the proportion of true null hypotheses is compared for: Benjamini and Hochberg's lowest slope approach (LSL), Storey's ${\widehat{\pi }}_0\left(\lambda \right)$ estimate with λ selected via bootstrapping (Storey_boot), Storey and Tibshirani's smoother method (ST_smoother), Langass's nonparametric maximum likelihood approach (convest), and the proposed average estimate approach with fixed values of B = 5, 10, 20, 50, 100 and with B chosen via the bootstrapping procedure (B_boot). There are 1,000 simulated data sets, each with a total of m = 1, 000 hypothesis tests, for each value of π₀.

π0	0.5	0.6	0.7	0.8	0.9
	Estimates of π0
LSL	0.7151	0.7889	0.8561	0.9184	0.9683
Storeyboot	0.4814	0.5789	0.6765	0.7728	0.8660
STsmoother	0.4951	0.5939	0.6980	0.7993	0.8973
convest	0.4963	0.5938	0.6947	0.7921	0.8882
B = 5	0.5132	0.6113	0.7136	0.8086	0.9058
B = 10	0.5082	0.6084	0.7083	0.8045	0.9052
B = 20	0.5141	0.6128	0.7115	0.8076	0.9064
B = 50	0.5196	0.6175	0.7156	0.8106	0.9078
B = 100	0.5243	0.6210	0.7180	0.8122	0.9085
B boot	0.5195	0.6175	0.7148	0.8113	0.9082
	Standard deviation of π0 estimates
LSL	0.0323	0.0269	0.0225	0.0155	0.0092
Storeyboot	0.0467	0.0491	0.0513	0.0522	0.0549
STsmoother	0.0513	0.0570	0.0608	0.0654	0.0656
convest	0.0331	0.0364	0.0337	0.0321	0.0328
B = 5	0.0335	0.0356	0.0420	0.0428	0.0382
B = 10	0.0391	0.0390	0.0402	0.0412	0.0366
B = 20	0.0331	0.0343	0.0358	0.0371	0.0331
B = 50	0.0293	0.0309	0.0321	0.0334	0.0315
B = 100	0.0272	0.0291	0.0307	0.0321	0.0312
B boot	0.0301	0.0301	0.0313	0.0313	0.0311

where ${\widehat{\pi }}_{0i}$ estimates π₀ for the ith simulation, and π₀ is the true value. As demonstrated, the LSL approach has the largest upward bias which guarantees that Benjamini and Hochberg's adaptive FDR controlling procedure controls the FDR below a pre-chosen FDR level. However, the FDR can be much lower than the pre-chosen FDR level. The LSL approach also has the smallest variation. The last three approaches [2–4] all underestimate the proportion of true null hypotheses. The proposed average estimate approach provides estimates of π₀ that have upward but relatively small bias and relatively small variance regardless of whether B is fixed or automatically chosen via bootstrap procedure. When B increases, the bias increases and the variation decreases. Both the small upward bias and small variance provide evidence that the proposed average estimate approach has better properties when compared to the other approaches.

The average of the true false discovery rate (FDR) from 1000 simulations is also compared in this simulation study by applying Benjamini and Hochberg's adaptive FDR controlling procedure (Benjamini and Hochberg, 2000) with π₀ estimated using the above mentioned five methods (Fig. 1). The FDR significance level was chosen as α = 0.05. For the purpose of comparison, the original BH FDR controlling procedure (Benjamini and Hochberg, 1995) and the adaptive FDR controlling procedure with the incorporation of the true value of π₀ were also applied to the p-values. It can be seen that the original BH FDR controlling procedure has the lowest FDR as expected. Because Benjamini and Hochberg's lowest slope approach overestimates π₀, the FDR is below, but much lower than, the pre-chosen α, although this approach has a bigger FDR than that of the BH procedure. Storey's bootstrap estimate, the smoother estimate and convest estimate produce higher FDRs than the pre-chosen level because all three methods underestimate π₀. Our proposed average estimate approach overestimates π₀, its FDR is below but very close to the pre-chosen significance level α = 0.05. Table 1 also demonstrates that the FDR for the proposed average estimate has the relatively small variation.

Figure 1

Simulation results of the False Discovery Rate (FDR) at significance level α = 0.05 for seven procedures: Benjamini and Hochberg's FDR controlling procedure with incorporation of the true π₀ (BHπ₀), Benjamini and Hochberg's FDR controlling procedure (BH), Benjamini and Hochberg's adaptive approach with incorporation of the estimate of π₀ which is estimated by the proposed average estimate procedure where B is chosen via bootstrapping (B_boot), Benjamini and Hochberg's lowest slope approach (LSL), Storey's bootstrapping approach (Storey_boot), Storey and Tibshirani's smoother method (ST_smoother), and Langass et al.'s nonparametric maximum likelihood estimate (convest), respectively. The black straight line represents FDR = 0.05. The total number of hypotheses tests is m = 1, 000 and the size of simulation study 1,000 for each value of π₀.

The power of the five adaptive FDR controlling procedures is compared (Fig. 2). The power of a procedure is measured by average power which is defined to be the ratio of average number of correct rejections of true alternative hypotheses to the total number of true alternative hypotheses. Formally, average power = E(S)/(m – m₀). As illustrated, the power decreases when π₀ increases for all of the FDR controlling procedures. The original BH procedure has the lowest power, while Benjamini and Hochberg's adaptive procedure has the second lowest power. It is not surprising that Storey_boot procedure has the largest statistical power, because the FDR of this procedure exceeds the pre-chosen FDR significance level (Fig. 1).

Figure 2

Simulation results for the evaluation of statistical power at significance level α = 0.05 for seven procedures: Benjamini and Hochberg's FDR controlling procedure with incorporation of the true π₀ (BHπ₀), Benjamini and Hochberg's FDR controlling procedure (BH), Benjamini and Hochberg's adaptive approach with incorporation of the estimate of π₀ which is estimated by the proposed average estimate procedure where B is chosen via bootstrapping (B_boot), Benjamini and Hochberg's lowest slope approach (LSL), Storey's bootstrapping approach (Storey_boot), Storey and Tibshirani's smoother method (ST_smoother) and Langass et al.'s nonparametric maximum likelihood estimate (convest), respectively. The total number of hypotheses tests is m = 1, 000, and the size of simulation study is 1,000 for each value of π₀.

Microarray Data Application

The same five estimating π₀ methods were also applied to the training samples of the leukemia data of Golub et al. (1999), which consist of 27 patients with acute lymphoblastic leukemia (ALL) and 11 patients with acute myeloid leukemia (AML). The samples were assayed using Affymetrix Hgu6800 chips and the gene expression data of 7129 genes (Affymetrix probes) are available from R library golubEsets. For each gene, a simple two-sample t-test was employed for testing differential gene expression and the p-value was computed. Table 2 gives the estimate of the proportion of true null hypotheses and the number of statistically significant genes.

Table 2

The estimate of the proportion of true null hypotheses and the number of statistically significant genes for the leukemai data (Golub et al. 1999) at significance level α = 0.05 after applying Benjamni and Hochberg's adaptive FDR controlling procedure with π₀ estimated using five methods: Benjamini and Hochberg's lowest slope approach (LSL), Storey's ${\widehat{\pi }}_0\left(\lambda \right)$ estimate with λ selected via bootstrapping (Storey_boot), Storey and Tibshirani's smoother method (ST_smoother), Langass's convest approach (convest), and the proposed average approach with B chosen via the bootstrapping procedure (B_boot). A two-sample t-test was used to compute the p-values.

Method	Estimate of π0	Number of Signicant genes
LSL	0.899	584
Storeyboot	0.595	787
STsmoother	0.583	791
convest	0.595	787
B boot	0.604	776

From this real data analysis, it can be seen that the Benjamni and Hochberg's LSL approach conservatively overestimates π₀, hence it leads to lowest power in terms of the number of rejections. Our proposed average approach provides a slightly larger estimate than Storey's bootstrap approach, the smoother estimate and the nonparametric maximum likelihood approach (convest), even though they end up with a similar number of rejections.

Summary

As array technology improves, it is anticipated that the number of features per array will only increase, hence multiple comparisons will continue to be a challenging problem. Specific to microarrays, the false discovery rate (FDR) is preferred to familywise error rate (FWER) because the FDR controlling procedures have more statistical power than the FWER controlling procedures, even at the cost of a few more type I errors (i.e. false positives). Since Benjamini and Hochberg (1995) proposed their FDR controlling procedure, a variety of methods have been proposed to estimate π₀ the proportion of true null hypotheses. As seen here, overestimating π₀ controls the FDR below the specified rate. When our and others, estimate of π₀ is incorporated into the Benjamini and Hochberg's FDR controlling procedure, the adaptive FDR controlling procedure has more power and an FDR close to the pre-chosen level.

In this work, we have compared several methods for estimating the proportion of true null hypotheses (π₀). Benjamini Hochberg's lowest slope approach (Benjamini and Hochberg, 2000) overestimates % Storey's estimate ${\widehat{\pi }}_0\left(\lambda \right)$ (Storey, 2002) also overestimates π₀ for any fixed value 0 ≤ λ < 1. When λ ↠ 1, the bias becomes smaller, and the variance becomes bigger. In order to find the optimal λ such that ${\widehat{\pi }}_0\left(\lambda \right)$ has small variation, Storey proposed a bootstrapping method (Storey et al. 2004). However, this method underestimates π₀ and the downward bias increases as the true value π₀ gets bigger. Storey and Tibshirani (2003) proposed a smoother method to estimate $\underset{\lambda \to 1}{\lim }{\widehat{\pi }}_0\left(\lambda \right)$ such that this estimate has small bias. Unfortunately, this method also underestimates π₀, although the bias is very small. Furthermore, the variation of this estimate is relatively large, which makes the adaptive FDR controlling procedure unstable. More recently, Langaas et al. (2005) proposed an estimate based on the nonparametric maximum likelihood function of the p-value density restricted to convex decreasing densities. However, this method also underestimates π₀, most likely because the distribution of the p-values is not decreasing for large p-values and tends to be flat.

Using the limitations of the existing approaches for estimating π₀ as the motivation, we propose the average estimate approach by taking average of the estimates of π₀ over a range of equally spaced points on the interval [0,1]. While our average estimate approach has a slightly larger bias, it also has smaller variation than any of the other methods. Furthermore, when compared to the other methods it is easy to implement (e.g. Excel) when the number of points used in approach is fixed (say, B = 10), and can be automated to choose B via a bootstrap procedure (R code available: http://www.stat.purdue.edu/~doerge). When our proposed estimated value of π₀ is incorporated into Benjamini and Hochberg's adaptive FDR controlling procedure, more statistical power is gained such that the FDR can be controlled below, yet extremely close to a desired level α.