A Unified Computational Framework to Compare Direct and Sequential False Discovery Rate Algorithms for Exploratory DNA Microarray Studies

The problem of detecting differential gene expression with microarray data has led to further innovative approaches to controlling false positives in multiple testing. False discovery rate (FDR) has been widely used as a measure of error in this multiple testing context. Direct estimation of FDR was recently proposed by Storey (2002, Journal of the Royal Statistical Society, Series B 64, 479-498) as a substantially more powerful alternative to the traditional sequential FDR controlling procedure, pioneered by Benjamini and Hochberg (1995, Journal of the Royal Statistical Society, Series B 57, 289-300). Direct estimation to FDR requires fixing a rejection region of interest and then conservatively estimating the associated FDR. On the other hand, sequential FDR procedure requires fixing a FDR control level and then estimating the rejection region. Thus, sequential and direct approaches to FDR control appear very different. In this paper, we introduce a unified computational framework for sequential FDR methods and propose a class of more powerful sequential FDR algorithms, that link the direct and sequential approaches. Under the proposed unified compuational framework, both approaches simply approximate the least conservative (optimal) sequential FDR procedure. We illustrate the FDR algorithms and concepts with some numerical studies (simulations) and with two real exploratory DNA microarray studies, one on the detection of molecular signatures in BRCA-mutation breast cancer patients and another on the detection of genetic signatures during colon cancer initiation and progression in the rat.


Introduction and Background
Since the introduction of DNA microarray technology, microarray applications in research have flourished, especiallly in biomedical research.See Nguyen et al. (2002) for a more thorough description of DNA microarrays and their applications.One common application in microarray studies involves identifying differentially expressed genes between two or more biological conditions.An example, which we will consider in this paper, involves identifying differentially expressed genes in hereditary breast cancer patients with mutation in the BRCA1 gene relative to patients with BRCA2-mutation.In this context, the null hypothesis of no differential expression is tested for each gene.Because the number of genes (or gene probes) are in the thousands, controlling for errors in this multiple testing situation is very important.Thus, methods that are able to identify truly alternative hypotheses (i.e., differentially expressed genes under a defined experimental condition) with few false positives are highly desirable.The application of such methods can result in a substantial reduction in research costs and efforts during the post-analysis and follow-up validation phase of microarray experiments (Chuaqui et al., 2002).
Table 1: Notations for the possible outcomes of testing m hypotheses (Benjamini and Hochberg 1995).The proportion of true null hypotheses is π 0 ≡ m 0 /m and FDR = E( V R I {R>0} ).

Accept Reject Total
Null true A particularly promising measure of error in multiple testing is the false discovery rate (FDR), the expected proportion of false discoveries among R discoveries or rejections, introduced by Benjamini and Hochberg (1995).More precisely, FDR is defined as where V is the number of erroneous rejections (Type I errors), R is the total number of rejections and I {A} denotes the indicator function for event A. Please refer to Table 1, where the outcomes for testing m hypotheses are summarized.Note in Table 1 that only R, W , and m are observable and all other quantities in the table are not observable.Also note that the well-known familywise error rate (FWER), the probability of rejecting any null hypothesis erroneously, is Pr(V > 0).Thus, FDR provides a much less strict criterion for control than the FWER in multiple testing.Hence, an obvious substantial gain in power is expected when controlling FDR compared to controlling FWER (Benjamini and Hochberg, 1995).The FWER is inappropriately strict for exploratory microarray studies, where m is in the thousands.
The traditional sequential approach to FDR, introduced by Benjamini and Hochberg (1995) (herein BH), requires fixing a FDR level of control, say α (0 < α < 1).Denote the ordered observed p-values as p (1) , . . ., p (m) .The Benjamini and Hochberg FDR (BH-FDR) controlling procedure is to find kBH = max j : p (j) ≤ j m α (1.2) and reject p (1) , . . ., p ( kBH ) , where α ∈ (0, 1) is the pre-specified target control level.BH proved that procedure (1.2) results in FDR ≤ π 0 α for 0 ≤ m 0 ≤ m, where π 0 = m 0 /m is the proportion of true null hypotheses.(It was later shown by Finner and Roters (2001) that the procedure (1.2) gives precisely FDR = π 0 α.) Since 0 ≤ π 0 ≤ 1, it follows that FDR is controlled at level α for all configuration of m 0 .However, note that the level of control is actually π 0 α, which is less than or equal to α.Thus, the BH-FDR controlling procedure (1.2) is increasingly conservative as π 0 approaches zero.This leads to a loss in power to detect true alternative hypotheses.Therefore, incorporating a less conservative, hence, more precise estimate of π 0 into the FDR controlling procedure (1.2) can improve power, as was done in Benjamini, Krieger, and Yekutieli (2001) and Benjamini and Hochberg (2000).Storey (2002;2003) and Storey and Tibshirani (2003) also recognized that estimation of π 0 is critical in the context of DNA microarray applications; however, they proposed a non-sequential, direct approach to estimate the FDR for a fixed rejection region.Storey (2002) proposed a direct estimate of FDR for a fixed rejection region [0, γ].More precisely, the proposed estimator of FDR is where Pr(P ≤ γ) = #{p j ≤ γ}/m and π0 (λ) is a conservatively biased estimator of π 0 and P is a random p-value from any test.Under independence, the main case considered in this paper, P comes from the null distribution (i.e.P ∼ Uniform[0, 1]) with probability π 0 .P is from the alternative distribution with probability 1 − π 0 .We will elaborate on the estimator π0 (λ) in Section 2.2.
We note that the benefit of estimating π 0 (or equivalently m 0 ) in multiple testing has previously been recognized, and dates back to at least Schweder and Spjøvtoll (1982).More theoretical studies of the sequential FDR procedure (1.2) and direct FDR estimation (1.3) can be found in Genovese and Wasserman (2002) and Storey, Taylor, and Siegmund (2004) respectively.
In this paper we propose a unified computational framework for numerical studies of sequential FDR methods for independent p-values based on estimation of π 0 .We also introduce a more powerful class of sequential FDR algorithms.We show that the BH-FDR procedure along with other sequential procedures, such as the two stage FDR procedure (Benjamini, Krieger, and Yekutieli, 2001), fall within this class.In addition, we illustrate that when using the same estimate of π 0 , the power of sequential FDR methods are very similar to the new, direct estimates of FDR proposed by Storey (2002).
The organization of the paper is as follows.In Section 2 we introduce a unified computational framework for comparing the FDR procedures.Under this framework, all sequential FDR methods approximate the least conservative FDR procedure through estimation of π 0 .This provides a unified framework for numerical comparisons of sequential and direct methods.We introduce a new, more powerful, family of sequential FDR algorithms in Section 2.2.Next, we illustrate these sequential FDR algorithms using two DNA microarray gene expression data sets in Section 3. In section 4 we describe a simulation framework for exploratory DNA microarray studies.Using this simulation, a comparison of the power and FDR control for the proposed sequential FDR algorithms to other sequential FDR methods is also given in Section 4. In Section 5 we compare the power of the proposed sequential FDR algorithms to the direct estimation of FDR and illustrate that the sequential and direct estimation procedures are essentially "equivalent" (in terms of power).

A Unified Computational Framework for FDR Procedures
For the original sequential FDR controlling procedure (1.2) proposed by BH, namely BH-FDR, we have FDR = π 0 α.Setting π 0 to its upper bound of one gives the desired level of FDR control, α.Thus, for the BH-FDR procedure, we can define π0 (BH) ≡ 1 because no information about π 0 was actually utilized from the distribution of observed p-values, {p j } m j=1 .Thus, the choice of π0 (BH) = 1 represents the most conservatice choice.The extreme opposite to this most conservative choice would be to use π 0 itself in the BH-FDR procedure.Clearly, this is the optimal (benchmark) choice and it is the least conservative FDR procedure.Of course, π 0 is unknown in practice.However, for numerical (or computational) studies it can be used as a benchmark.More precisely, for numerical studies of the performance (in terms of power) of the various FDR procedures, the choice of π0 (BH) = 1 and π 0 itself provide the lower and upper bound on the performance of any given FDR procedure.Thus, in this section, we introduce a unified computational view of sequential FDR procedures based on incorporation of more precise information on π 0 from the data (Section 2.1).Next, the direct estimation approach to FDR is casted into this framework via estimation of π 0 (Section 2.2).This is possible since both sequential and direct approaches to FDR approximate the optimal FDR procedure.

Approximating the least conservative FDR controlling procedure
Since the original BH-FDR provides FDR = π 0 α ≤ α, it is conservative by a factor of π 0 = m 0 /m.If π 0 (m 0 ) is known, then the conservativeness can be corrected by applying the BH-FDR procedure at level α = α/π 0 , instead of α.This correction provides FDR control at level α, since FDR = π 0 α = α.Thus, we define the optimal or least conservative FDR (LC-FDR) procedure, assuming that α is known, as finding kLC = max j : and then rejecting the hypotheses corresponding to the p-values p (1) , . . ., p ( kLC ) .
The LC-FDR procedure, although useless in practice because π 0 is unknown, provides the benchmark for studying the precision of estimating FDR and the power to detect true positives (or true alternative hypotheses).Thus, we can define a class of sequential FDR controlling procedures that approximate the LC-FDR procedure as finding k = max j : and then rejecting p (1) , . . ., p ( k) such that FDR ≤ α, where π0 is a conservative estimate of π 0 , α * ∈ (0, 1), and α ∈ (0, 1) is the target FDR level of control.Under this framework, the least conservative data-based estimate of π0 such that FDR ≤ α is desirable.For example, the original sequential FDR controlling procedure, namely BH-FDR, trivially falls into class (2.2) with π0 given by π0 (BH) ≡ 1 and α * = α: .
As pointed out earlier, it is the most conservative FDR controlling procedure according to (2.2).
Since the choice of α * = α is more strict than α * = α (i.e.α < α), BKY also proposed the following two-stage modified FDR (2SM-FDR) procedure kBKY-M = max j : , where π0 (BKY-M) is the first stage estimate of π 0 at level α instead of α .Note that under this unified view, given by (2.2), FDR methods, such as the BH-FDR, 2S-FDR, and 2SM-FDR, essentially "mimic" (or approximate) the least conservative FDR process.More precisely, they use estimates of the form α * /π 0 (•) as a plug-in for α/π 0 , where α * is some pre-chosen level.Thus, we can view these FDR controlling procedures as approximations to the LC-FDR procedure.In this setting, the least conservative data-based estimate of π 0 such that FDR ≤ α is desirable.Next, we describe a new family of sequential FDR algorithms that provide a better approximation of the LC-FDR procedure.This procedure uses the estimator of π 0 proposed by Storey (2002) in his direct estimation approach to FDR.

A "new" family of sequential FDR algorithms: Linking direct FDR estimation to sequential FDR procedures
A gain in power to detect true alternative hypotheses will require a less conservative estimate of π 0 .This gain in power was first shown by Storey (2002) for independent p-values and illustrated with extensive simulation by Nguyen (2004), including cases involving violation of key assumptions of the FDR framework.In this section we propose a family of sequential "FDR algorithms" by utilizing the estimator of π 0 proposed by Storey (2002) and further implemented by Storey and Tibshirani (2003) in their direct estimation approach to FDR.
The estimator of π 0 proposed by Storey (2002) is as follows.Since it is much more likely that very large p-values correspond to true null hypotheses, consider the set of large p-values falling into the upper interval (λ, 1] to estimate π 0 (for some chosen 0 < λ < 1).Furthermore, note that if no genes are differentially expressed, then the null p-values are uniformly distributed, denoted P ∼ U(0, 1).Hence, Pr{P ∈ (λ, 1]} = Pr(P > λ) = 1 − λ.It follows that the expected number of null p-values that would fall into the interval (λ, 1] is (1 − λ)m 0 .In addition, if we know the number of null p-values in (λ, 1], #{Null p j > λ}, then an unbiased estimate of π 0 is π0 The numerator of (2.4), #{Null p j > λ}, is not observable in practice.However, replacing the numerator with #{p j > λ}, an observable quantity, leads to a conservatively biased estimate of π 0 : It can be seen that the estimate, π0 (λ), is conservatively biased from the following simple inequality Thus, E[π 0 (λ)] ≥ E[π 0 (UB)] = π 0 .Note that the index parameter, λ, is actually a tuning parameter which balances bias and variance.More precisely, as λ approaches 1, #{p j > λ} consists mostly of truly null p-values; therefore, the bias decreases.However, the interval used to estimate π0 (λ), specifically (λ, 1], shrinks to zero as λ → 1; hence, the variance increases.Thus, Storey and Tibshirani (2003) proposed an automatic algorithm for choosing the optimal λ to minimize the mean squared error of π0 (λ).The optimal estimator of π 0 is denoted as π0 (OPT) and is reproduced in the Appendix for convenience.
Using the estimator π0 (λ), we propose the following family of sequential FDR algorithms (indexed by λ) to better approximate the LC-FDR controlling procedure (2.1) kλ = max j : p (j) ≤ j m α π0 (λ) . (2.6) The proposed FDR algorithm simply replaces π0 (BH) ≡ 1 in the ordinal BH-FDR procedure with π0 (λ), a less conservative estimate of π 0 .We emphasize that the FDR algorithm given by (2.6) is sequential and it utilizes the exact information (estimator) of π 0 as with the direct estimation approach.Therefore, one might expect that the power of the two approaches would be equivalent.We will elaborate further on this in Section 5, using numerical studies.
We illustrate in the next section the proposed sequential FDR algorithms with two exploratory DNA microarray studies involving: (1) the detection of molecular signatures in BRCA-mutation breast cancer patients and ( 2) the detection of genomic signatures during colon cancer initiation and progression in the rat.

Detection of molecular signatures in BRCA-mutation breast cancer patients
One common application of DNA microarray technologies in biomedical research is the detection of differentially expressed genes between two or more biological conditions (groups).For example, Hendenfalk et al. (2001) applied cDNA microarray to the study of hereditary breast cancer.In particular, one goal of the microarray study was to identify genes, among m = 3, 226 genes, that are differentially expressed in breast cancer patients with mutations in the BRCA1 gene relative to patients with BRCA2-mutations.There were n 1 = 7 patients with BRCA1-mutation and n 2 = 8 patients with BRCA2-mutation.For each gene, we computed the two sample t-statistic and the corresponding p-value for testing the null hypothesis of no differential gene expression.Figure 3A displays the density histogram of the observed p-values, {p j } m j=1 .For this distribution of p-values, it is estimated that the proportion of true null hypotheses, π 0 , is π0 (λ) = 0.678 using λ = 1/2.The optimal estimate, π0 (OPT), gave a similar estimate of 0.668.We applied sequential FDR methods, namely BH-FDR, 2S-FDR, 2SM-FDR, and the proposed FDR algorithm (2.6) using π0 (λ).For illustration, we applied the methods to control false discoveries among the m tests at 10% (α = .10).

Detection of genomic signatures during colon cancer initiation and progression in the rat
In a study of the chemopreventive n-3 polyunsaturated fatty acids (PUFA) as it relates to colon cancer initiation and progression in the rat, Davidson, Nguyen et al. (2004) used CodeLink TM oligonucleotide microarrays (Ramakrishnan et al., 2002) to monitor gene expression profiles from 93 rat tissue samples.Each microarray contains 9,685 gene probes.The study utilized a 3 × 2 × 2 factorial design with three types of dietary fat (rich in n-6 PUFA, n-3 PUFA or n-9 monounsaturated fatty acid (MUFA)), two treatments (injection of carcinogen azoxymethane (AOM) or saline) and two time points (12 hours (initiation) and 10 weeks (progression) post first injection).93 RNA samples from 59 rats (including 34 replicates) were randomly assigned into the three-factor combination.It is of interest in the study to examine the molecular mechanisms by which dietary fat composition exerts tumor enhancing or inhibiting effects.More specifically, it is of interest to identify gene expressions that are significantly different between n-3 PUFA and n-6 PUFA enriched diets.For illustration, Figure 2A displays the distribution (density histogram) of the 9,685 observed p-values contrasting n-3 PUFA and n-6 PUFA enriched diets from a linear mixed model for repeated measures.Applying Storey and Tibshirani (2003) estimator for π 0 , the the proportion of truly non-differentially expressed genes between n-3 PUFA and n-6 PUFA enriched diets, gives π0 (OPT) = 0.758.Using this estimate, the proposed FDR method results in 31 differetially expressed genes at FDR level α = 0.05.The original BH-FDR and the two-stage FDR procedures both yield 22 genes declared differetially expressed at the same FDR level.

Simulation Framework for Exploratory DNA Studies: Some Computational Experiments and Results
In this section we describe the simulation framework to compare FDR algorithms for exploratory DNA microarray studies.Using this framework, we compare the power of the proposed sequential FDR algorithm (2.6) to detect true alternative hypotheses.We illustrate via simulation the gain in power from using the proposed sequential FDR algorithm relative to the BH-FDR, 2S-FDR, 2SM-FDR, and, more importantly, to the least conservative FDR (LC-FDR) procedure.Recall that the LC-FDR procedure has optimal power and FDR control.We also examine the FDR control for the proposed FDR algorithm (2.6).

Simulation design considerations
The simulation was designed in the context of a two-sample comparison, anal-ogous to the BRCA-mutation breast cancer example in Section 3.More specifically, we generated an m × n gene expression data matrix,

X = X
(01) Of the total m genes, m 0 are truly null and the remaining m 1 = m − m 0 are truly alternative.
Under the null setting, we generated the expression value for gene j in both groups i = 1 and 2, independently from a N(µ 0 , σ 2 0 ) distribution: x ji ∼ N(µ 0 , σ 2 0 ) for j = 1, . . ., m 0 and i = 1, 2. This is the first m 0 rows of the data matrix X above.Under the null hypothesis, there is no difference in gene expression between groups 1 and 2 for all m 0 genes.
However, there is a difference in gene expression between groups 1 and 2 under the alternative setting.Thus, expression values for gene j in group 2, x j2 (j = m 0 + 1, . . ., m), were independently generated from a N(µ 1 , σ 2 1 ) distribution.For group 1, x j1 ∼ N(µ 0 , σ 2 0 ).Rather than setting µ 1 to a fixed value, the true alternative mean expression for group 2 was allowed to vary above the null mean value of µ 0 .Also, in order to more reflect the heterogeneity of variance often encountered with microarray data in practice, we similarly allowed the alternative variance parameter, σ 2 1 to vary.Table 2 summarizes the simulation model and parameters.
In this simulation framework, our intent is to test the null hypothesis that there is no differential expression for each gene between groups 1 and 2. Thus, after the data generation we computed the two sample t-statistic and the corresponding p-value for each gene, as was done in the BRCA-mutation breast cancer example.Because the data was generated with heterogeneity of variance, we examined both sets of p-values obtained from (1) erroneously assuming a t distribution for the null distribution and (2) using a permutation method to approximate the null distribution.For (1) the p-value for gene j was computed as p j = Pr(|t(n − 2)| > t j ), where t j is the observed t-statistics and t(n − 2) denotes the t distribution with n − 2 degrees of freedom.For (2), we obtained B permutations of the n sample labels.For the bth permutation, denote the re-calculated t-statistics based on the permuted data by {t 0b j } m j=1 (b = 1, . . ., B).The permutation-based p-values can be computed as (Storey and Tibshirani, 2003).
For example, Figure 3B displays the density histogram of the observed pvalues for a simulated data set with m = 1, 000 genes, group sample size n 1 = n 2 = 8, and π 0 = m 0 /m = 0.70.Note that the simulated data resembles closely the BRCA-mutation hereditary breast cancer data described earlier in Section 3. The estimate of the proportion of true null hypotheses is π0 (λ) = 0.710 for the choice of λ = 1/2.We summarized below some basic points of the simulation model which we believe more accurately reflect real microarray data.
• The number of genes, m, is large.
• The range of the sample size is small to moderate and reflects the range afforded in real microarray studies.
• The model encompasses heterogeneity of variance, often encountered with real microarrray data.This applies to both the null and alternative settings.
• The mean expression for genes differentially expressed is allowed to vary.Again, this applies to both the null and alternative settings.
• The above two items (alone and/or together) imply that the sampling distribution of the test statistic is not exact.Therefore, the p-values generated are only approximations.
• Biological theory and assumptions imply dependence in gene expression, so violation of independence should be considered in the simulation model.This was done, for example, in Nguyen (2004) and we refer the reader there for details.
• Measurement error are well recognized with gene expression measurements and can be incorporated into the simulation model.(For example, see Zien et al., 2002).Although measurement error is an important issue with microarray data, we will not address this point further and will address it in future work.However, the framework described in Section 2 can be applied to any simulation model, including models with measurement errors and dependent data.
Under these more realistic simulation conditions, it will be interesting to examine the power and FDR control.For illustration, we will describe the results for the model given in Table 2 with Both µ 1 and σ 2 1 are randomly selected for each iteration (data set simulated).Also, we only describe the case where n i = 8 per group, similar to the BRCA-mutation data set.Details and the patterns of the other cases are similar.However, some of the above cases, including violation of independence, non-normal models for gene expression, and different configuration of sample sizes can be found in Nguyen (2004), although not under the unified computational framework presented here.
Table 2: Simulation model and parameters.Data were generated such that there is no differential expression in m 0 genes from a total of m genes.

Power
Figure 4A displays the true proportion of null hypotheses, π 0 , and its various conservative estimates, namely π0 (BH), π0 (BKY), π0 (BKY-M), and π0 (λ) based on 10,000 simulations for each π 0 ∈ {0.1, 0.2, . . ., 0.9}.Thus, a total of 90,000 data sets were generated.Also, for this simulation experiment we used a sample size of n i = 8 per group, similar to the BRCA microarray data.It is clear from Figure 4A that π0 (λ) is the least conservative estimate of π 0 ; hence, much closer to the target, π 0 .The resulting power, for the same simulation experiment, is given in Figure 4B.It is evident that the proposed FDR algorithm using π0 (λ) is substantially more powerful than the other sequential FDR methods.In addition, it is much closer to the optimal power, given by the LC-FDR procedure.
We note that there is little difference between the two stage procedure, 2S-FDR, and the modified version, 2SM-FDR.

FDR control
For the original BH-FDR controlling procedure, FDR = π 0 α; hence, it is conservative because it actually controls FDR at a lower level of π 0 α, rather at the target level of α.It follows that the actual level of FDR control for the BH-FDR procedure is linearly decreasing from α, as π 0 → 0. On the other extreme, the least conservative or optimal procedure, namely LC-FDR, completely corrects for this conservativeness by utilizing π 0 itself.Thus, as expected and confirmed by the results in Figure 5, the LC-FDR control is at the target level of α and the BH-FDR control is lower, at precisely π 0 α.(For illustration, the results in Figure 5 is for a target FDR control level of α = 0.05.)Also, the two-stage procedures control FDR below the target α; however, they are nearly as conservative as the BH-FDR control, the most conservative case.As can be seen from Figure 5, our proposed FDR algorithm (2.6), which uses Story's less conservative estimate, π0 (λ), also controls FDR at the target level of α = 0.05.In other words, FDR ≤ α for the proposed FDR algorithm (2.6).However, the FDR control is substantially less conservative than the other procedures, thus affording a large gain in power, as was described earlier.Note that he BH-FDR control is lower, at level π 0 α.The proposed FDR algorithm using π0 (λ) is less conservative and, hence, closest to the target control level.Note that he BH-FDR control is lower, at level π 0 α.The proposed FDR algorithm using π0 (λ) is less conservative and, hence, closest to the target control level.

A "Reconciliation" Between Sequential and Direct FDR Approaches: Diminishing the Power Gap
Multiple testing in microarray applications often involves conducting thousands of tests simultaneously to screen for differentially expressed genes, as illustrated by the BRCA breast cancer data in Section 3. The widely popular sequential method, namely BH-FDR, completely ignores important information regarding π 0 .Hence, it critically lacks power to detect truly alternative hypotheses, especially in microarray applications where m is extremely large.In fact, increasing information regarding π 0 becomes available when m → ∞, and this is precisely the case in micrarray applications.However, the landmark paper introducing FDR by Benjamini and Hochberg (1995) did not envision applications with this scale.This is clearly apparent from their numerical study with m ≤ 64.When m is extremely large, the investigator often must accept a high level of false discovery rate (α) in order to obtain some rejections or significant genes for follow-up studies under the BH-FDR procedure.Storey (2002) recognizes some of the shortcomings of the BH-FDR procedure in microarray applications.He proposed estimating π 0 and FDR directly for a fixed rejection region.For example, if the investigator decides, a priori, to reject all genes with p-values less or equal to γ = 0.005, what is her/his expected FDR?
The direct approach provides a conservatively biased estimator, FDR(γ).Thus, the investigator fixes, before hand, the rejection; in this example, the rejection region is [0, γ].The direct estimation approach to FDR (Storey, 2002) is an important and fundamental departure from the tradition sequential FDR approach in the following sense.Essentially, traditional sequential methods conservatively estimate the rejection region for a fixed FDR level, and the direct estimation approach conservatively estimate FDR directly for a fixed rejection region.
Thus, the interpretation and implementation of the two approaches are different.Differing views, advocating the practice of one approach over the other remains (Storey, 2002).However, based on numerical studies, there is a substantial gain in power from using direct estimation relative to using the sequential BH-FDR procedure.(See Figure 1 and Table 2 of Storey (2002).)This is due to the fact that the direct estimate of FDR uses a better estimate of π 0 , and not to whether one fixes the rejection region (direct method) or the FDR level (sequential method).As we will illustrate, the key lies in utilizing a less conservative, hence more precise, estimate of π 0 .More precisely, if we use exactly the same information (estimate of π 0 ) in both the sequential and direct FDR methods, will the differences in power of the two approaches diminish?Storey (2002) proposed a direct estimate of FDR for a fixed rejection region [0, γ], which uses π0 (λ).As described in the Introduction section, the direct estimator of FDR is FDR λ (γ) = π0 (λ) γ Pr(P ≤ γ) , where Pr(P ≤ γ) = #{p j ≤ γ}/m.We compare the power curve of FDR λ (γ) to the proposed sequential FDR algorithm (2.6), which uses the same estimate of π 0 , namely π0 (λ).We note that sequential FDR methods cannot be directly compared with the direct estimation approach (Storey, 2002).This is because, as pointed out earlier, the former estimates the rejection region, whereas the later estimates FDR directly.However, this can be circumvented by using the sequential methods to control FDR at level α = FDR λ (γ) for each iteration (simulation run).The results, based on 10,000 simulation runs and using the same setup as described in Section 4, are given in Figure 6.Sequential methods using poor estimates of π 0 , namely BH-FDR, 2S-FDR, and 2SM-FDR, have very low power; hence their power curves fall far below the optimal power curve (LC-FDR).The power for the proposed sequential FDR method (2.6) is close to the power for the direct estimator, FDR λ (γ), because they both use the same estimate for π 0 .This result is not surprising in light of the connection between the BH-FDR procedure and FDR λ (γ) shown by Storey (2002) for independent p-values; BH-FDR is a special case of FDR λ (γ) when π0 (λ) = 1.Similarly, it is not difficult to show analytically that the proposed sequential FDR algorithm is equivalent to FDR λ (γ) for independent p-values.

Discussion
We have provided a unified computational framework for comparing FDR procedures in numerical studies, encompassing traditional sequential methods and the new direct approach to FDR.We have demonstrated that when using the same information regarding π 0 in both sequential FDR methods and in direct estimation of FDR, the power curves are equivalent.In particular, we introduced a class of sequential FDR of the form k = {j : p (j) ≤ (j/m)(α * /π 0 )} and showed that various sequential FDR methods fall within this class.In addition, when a less conservative estimate of π 0 is used, specifically π0 = π0 (λ), the power is equivalent to the direct estimator of FDR.
These conclusions also hold for other simulation configurations, in addition to those described earlier.For example, note that in the simulation we used a sample size of n i = 8 per group, similar to the BRCA data.In the context of two sample comparison in microarray experiments, this is a moderate sample size.However, the results hold for even smaller sample sizes.Of course, as n i grows larger, the power for all methods converge to the optimal power of the LC-FDR procedure.
We emphasize that the family of sequential FDR algorithms given by (2.6) was motivated and derived at by examining the estimation of π 0 in the sequential FDR procedures of Benjamimi and colleages.Thus, along this same line, utilizing the more precise estimator π0 (λ) of Storey's to improve power is sensible.However, it is important to note that the later algorithms do not guarantee the FDR control, in the sense of Bejamini and Hochberg (1955) for the BH-FDR procedure.Thus, one needs to take care to examine the FDR control in numerical studies (on average) as described in Section 4.3.This is particularly true also for numerical studies involving violation of the assumptions in the FDR framework.In these cases, even the procedures with proven FDR control no longer hold because of the violation of assumptions.Furthermore, it is important to note here that Storey (2002), more or less, already recognized the link between his direct estimation approach to the sequential approach early on.(See also Storey et al. (2004).)However, their numerical comparisons aims to demonstrate the improvement in power within the assumptions of the FDR framework.As described in Section 4 the statistical assumptions do not hold under the minimal conditions of real microarray studies, including inhomogenous mean and variance, small sample sizes, dependence in gene expression, and so on.Thus, it is important to examine the performance of both FDR approaches under these more realistic conditions and also gauge their performance relative the the (optimal) least conservative FDR procedure.The unified computation framework proposed here, base on estimation of π 0 and approximating the least conservative FDR procedure, is one way to make such comparisons.
Finally, we note that the real data analyses use the optimal choice of λ for direct estimation, namely π 0 (OPT).(See also the appendix section.)However, to reduce the computational cost in the large simulation we used the approximation π0 (λ = 1/2).A prelimary simulation suggests that the power of the direct estimation method using π 0 (OPT) has slightly higher power than the approximation π0 (λ = 1/2), as expected.Therefore, when comparing the power of direct estimation to sequential FDR methods, the results reported here remain valid.

Figure 1 :
Figure1: Sequential FDR methods applied to the BRCA data.The circles plotted are j/m versus p (j) .A sequential FDR method is equivalent to finding the first time in the (reversed) sequence of ordered p-values, {p (m) , p (m−1) , . . ., p (1) }, that crosses the line with slope α * /π 0 .The right solid vertical line indicates this crossing, at j = 281, for the proposed FDR algorithm (dotted line).The left vertical line indicate the crossing for both the BH-FDR and the 2S-FDR procedure, which occurred at the same location (j = 162).It is also the same for 2SM-FDR (not shown).

Figure 2 :
Figure 2: Sequential FDR methods applied to the colon cancer data.(A, top) Displayed is the density histogram of the 9,685 observed p-values comparing n-3 PUFA and n-6 PUFA enriched diets.(B, bottom) Displayed is a graphical summary of the sequential FDR analyses.Estimate of π 0 , the proportion of truly non-differentially expressed genes, is π0 ((OPT)) = 0.758.

Figure 5 :
Figure5: FDR control for sequential methods.Displayed are the observed levels of FDR control, averaged over the 10,000 simulation runs (α = 0.05).Note that he BH-FDR control is lower, at level π 0 α.The proposed FDR algorithm using π0 (λ) is less conservative and, hence, closest to the target control level.

Figure 6 :
Figure6: FDR control for sequential methods.Displayed are the observed levels of FDR control, averaged over the 10,000 simulation runs (α = 0.05).Note that he BH-FDR control is lower, at level π 0 α.The proposed FDR algorithm using π0 (λ) is less conservative and, hence, closest to the target control level.