Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084 www.elsevier.com/locate/jspi
Two-stage stepup procedures controlling FDR夡 Sanat K. Sarkar Department of Statistics, Temple University, Philadelphia, USA Received 10 June 2005; received in revised form 6 August 2006; accepted 23 March 2007 Available online 17 July 2007
Abstract A two-stage stepup procedure is defined and an explicit formula for the FDR of this procedure is derived under any distributional setting. Sets of critical values are determined that provide a control of the FDR of a two-stage stepup procedure under iid mixture model. A class of two-stage FDR procedures modifying the Benjamini–Hochberg (BH) procedure and containing the one given in Storey et al. [2004. Strong control, conservative point estimation and simultanaeous conservative consistency of false discovery rates: a unified approach. J. Roy. Statist. Soc. Ser. B 66, 187–205] is obtained. The FDR controlling property of the Storey–Taylor–Siegmund procedure is proved only under the independence, which is different from that presented by these authors. A single-stage stepup procedure controlling the FDR under any form of dependence, which is different from and in some situations performs better than the Benjamini–Yekutieli (BY) procedure, is given before discussing how to obtain two-stage versions of the BY and this new procedures. Simulations reveal that procedures proposed in this article under the mixture model can perform quite well in terms of improving the FDR control of the BH procedure. However, the similar idea of improving the FDR control of a stepup procedure under any form dependence does not seem to work. © 2007 Elsevier B.V. All rights reserved. Keywords: Modified BH procedure; Mixture model; Average power
1. Introduction Since its introduction by Benjamini and Hochberg (1995) in multiple testing, the concept of false discovery rate (FDR) as a measure of error rate has been receiving considerable attention due to its relevance in many scientific investigations, particularly in the context of DNA microarray analysis. The FDR is the expected proportion of type I errors among the rejected hypotheses. There are many procedures that control the FDR. Among them, the stepup procedures of Benjamini and Hochberg (1995) has received the most attention. Consider n null hypotheses H1 , . . . , Hn that are being tested simultaneously using the corresponding p-values P1 , . . . , Pn . Let these Pi ’s be ordered as 0 = P0:n < P1:n · · · Pn:n . Then, given a set of ordered critical values 0 = 0 < 1 · · · n < 1, a stepup procedure rejects Hi for all i KSU , where KSU = max{0 i n : Pi:n i }. Assuming that Pi ∼ U(0, 1), for each i = 1, . . . , n, under the null hypotheses, the critical values of the Benjamini–Hochberg (BH) procedure controlling the FDR at are given by i = i/n,
i = 1, . . . , n.
夡 The
research is supported by NSF Grants DMS-0306366 and DMS-0603868. E-mail address:
[email protected].
0378-3758/$ - see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2007.03.058
(1.1)
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
1073
One of the most important results of the BH procedure is that its FDR is exactly equal to (n0 /n), where n0 /n is the proportion of true null hypotheses, when the p-values are independent and is less than or equal to (n0 /n) when they are positive regression dependent on subset (PRDS) of those corresponding to the null hypotheses, a positive dependence property shared by many of the multivariate distributions arising in multiple testing (Benjamini and Yekutieli, 2001; Sarkar, 2002). Finner and Roters (2001) and Storey et al. (2004) offered different proofs of the above equality result under independence. Zhang (2005) contains some of these results. In the absence of a knowledge of any particular form of dependence among the Pi ’s, the procedure of Benjamini and Yekutieli (2001) is generally recommended. The critical values of this procedure (referred to as the BY procedure) are given by ⎧ n ⎫ ⎨ 1⎬ i = i , n ⎩ j⎭
i = 1, . . . , n,
(1.2)
j =1
with the FDR being less than or equal to (n0 /n). Since n0 is unknown, both BH and BY procedures could potentially be improved in the sense of providing a better control of the FDR if they are modified using an appropriately chosen estimate of n0 . This idea, nevertheless, has so far been implemented successfully only in modifying the BH procedure under the independence. Benjamini and Hochberg (2000) first proposed a data-adaptive modification of the BH procedure, modifying further a similar approach in Hochberg and Benjamini (1990) aimed at improving the familywise error rate (FWER) control of the Bonferroni single-step and related stepwise procedures. They borrowed the idea of estimating n0 based on the number of relatively large p-values outlined in Schweder and Spjotvoll (1982). This has also been used in Benjamini et al. (2006) and Storey et al. (2004) who provided alternative modifications to the BH procedure. While Benjamini et al. (2006) and Storey et al. (2004) offer proofs of the FDR controlling property of their modified BH procedures under independence of the p-values, no proof is yet available for the adaptive BH procedure of Benjamini and Hochberg (2000). These are all basically two-stage procedures that estimate n0 at the first stage before applying the BH procedure at the second stage with its critical values modified incorporating this estimate. Sarkar (2006) has recently modified the single-step Bonferroni procedure using the estimate of n0 in Storey et al. (2004) through a two-stage approach and proved assuming the independence that it controls the FDR. In this article, we obtain a wider class of two-stage modifications of the BH procedure under mixture model (Genovese and Wasserman, 2002; Storey, 2002) and obtain similar two-stage versions of the BY procedure and an alternative to it proposed for the first time in this article under any form of dependence. We start with a general two-stage stepup procedure in Section 2 and obtain in the same section an explicit expression for its FDR under any distributional setting. We then determine in Section 3, assuming the mixture model, sets of critical values of a two-stage stepup procedure that provide a control of the FDR at a pre-specified value ∈ (0, 1). The resulting two-stage FDR procedures represent a newer class of modified BH procedures, including the one given in Storey et al. (2004). For the Storey–Taylor–Siegmund procedure (to be referred to as the STS procedure), we also give a slightly more general result providing a proof of its FDR controlling property only under the independence, that is, without having the additional assumption on the probability distribution of different configurations of true and false null hypotheses. This is different from Storey et al. (2004). We choose one in our proposed class of procedures and study it numerically in comparison with the STS, the adaptive BH and the original BH procedures in terms of their FDR and average power. While our procedure appears to have the best performance for small difference between true and false nulls, the STS procedure seems to work the best when this difference is large. We obtain an inequality in Section 4 for the FDR of a two-stage stepup procedure under any form of dependence among the p-values. This provides an idea of developing the newer FDR procedure and leads to the developments of two-stage versions of this and the BY procedures. To set the stage for a numerical investigation of the question if a two-stage stepup procedure really improves the FDR control of its respective single-stage version in the dependence case, we decide to numerically assess how comparable our new FDR procedure is to the commonly used BY procedure. Interestingly, our procedure is seen to exhibit better control of the FDR when there is a high dependence among the p-values and not many of the null hypotheses are false. Unfortunately, however, none of these two-stage procedures does seem to work in terms of improving the FDR control of the single-stage procedures.
1074
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
2. Two-stage stepup procedure and its FDR In this section, we will define a two-stage stepup procedure and obtain an explicit formula for its FDR under any distributional setting. Definition 2.1 (Two-stage stepup procedure). Given 0 = 0 < 1 · · · n < n+1 = 1 and a set of other constants j k , 0 k j n, satisfying 0 = j 0 < j 1 · · · jj j , for each j, let J = max{0 i n : Pi:n i }, and K(J ) = max{0 i J : Pi:n J i }. Reject Hi for all i K(J ). Remark 2.1. A two-stage stepup procedure starts with a stepup test with the critical values 1 · · · n at the first stage and continues testing as long as no hypothesis is rejected. Once a rejection occurs, testing at the first stage stops, the n − J hypotheses accepted are declared true, and those which are not accepted are deferred to the second stage to be further tested using another stepup test with a smaller set of critical values. Given J = j , K(j ) is the number of hypotheses that are finally rejected. It reduces to a single-stage stepup procedure with the critical values 1 · · · n if j k = j ,1k j n. The critical values j ’s at the first stage are chosen to be relatively large compared to those at the second stage. The reason behind this is that the hypotheses that are accepted based on larger critical values (in terms of p-values) should not be further tested. Moreover, the p-values that are insignificant compared to larger critical values are better indicators of the true null hypotheses, providing a reliable information about the number of true null hypotheses that may be used to modify the critical values before testing the remaining hypotheses. Our definition of two-stage stepup procedure does not cover the type considered by Benjamini et al. (2006) who introduced a two-stage version of linear stepup procedure. Unlike ours, their procedure allows the first-stage critical values to be smaller than those at the second stage. This, however, presents an inconsistency in the procedure in that some hypotheses that are insignificant at the first stage may be found significant at the second stage. Our definition avoids such inconsistencies. It is actually more in the spirit of Storey et al. (2004) who restrict the significance threshold for their modified BH procedure to the region which is found significant at the first stage. Notice that we accept the hypotheses found insignificant at the first stage before proceeding to the second stage; whereas, Benjamini et al. (2006) make no decision at the first stage unless all the hypotheses are accepted or rejected and Storey et al. (2004) make no decision unless all the hypotheses are accepted. As seen from the following proposition, under certain conditions on the critical values of our procedure, it makes no difference whether or not we make any decision at the first stage. Proposition 2.1. With additional set of constant j k , 1 j < k n, satisfying jj j,j +1 · · · j n , for each j, let K ∗ (J )=max{0 i n : Pi:n J i } and replace K(J ) by K ∗ (J ) in Definition 2.1. This new two-stage stepup procedure will remain equivalent to the original one as long as j k k , 1 j k n. Towards deriving a formula for the FDR of a two-stage stepup procedure, first note that if J = 0, that is, if all the hypotheses are declared true in the first stage, the FDR is 0. Therefore, by letting {Hi , i ∈ I0 } the set of true null hypotheses, we can write a formula for the FDR of a two-stage stepup procedure as follows: FDR = =
j n 1 i∈I0 j =1 k=1 j n i∈I0 j =1 k=1
k
Pr{J = j, K(j ) = k, Hi is rejected}
1 Pr{J = j, K(j ) = k, Pi j k }. k
Note that, for j = 2, . . . , n, j 1 k=1
k
=
Pr{J = j, K(j ) = k, Pi j k }
j 1 k=1
k
Pr{J = j, K(j ) k, Pi j k } −
j 1 k=2
k
Pr{J = j, K(j ) k − 1, Pi j k }
(2.1)
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
j −1
I (Pi j k ) I (Pi j,k+1 ) 1 − . = Pr{J = j, Pi jj } + E Pr{J = j, K(j ) k|Pi } j k k+1
1075
(2.2)
k=1
(−i)
(−i)
Let P1:n−1 · · · Pn−1:n−1 denote the ordered components of the set obtained by removing Pi from (P1 , . . . , Pn ). Then, we have, for any j = 1, . . . , n, {J = j, Pi jj } = {Pj :n j , Pj +1:n > j +1 , . . . , Pn:n > n , Pi jj } (−i)
(−i)
(−i)
= {Pj −1:n−1 j , Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n , Pi jj }.
(2.3)
And, for any 1 k < j n and j,k+1 , {J = j, K(j ) k, Pi } =
k
{J = j, K(j ) = l, Pi }
l=0
=
k
{Pl:n j l , Pl+1:n > j,l+1 , . . . , jj Pj :n < j , Pj +1:n > j +1 , . . . , Pn:n > n , Pi }
l=0
= {Pk+1:n > j,k+1 , . . . , Pj −1,n > j,j −1 , jj Pj :n < j , Pj +1:n > j +1 , . . . , Pn:n > n , Pi } (−i)
(−i)
(−i)
(−i)
= {Pk:n−1 > j,k+1 , . . . , Pj −2,n−1 > j,j −1 , jj Pj −1:n−1 < j , Pj :n−1 > j +1 (−i)
, . . . , Pn−1:n−1 > n , Pi } (−i)
(2.4)
(−i)
(with P0:n−1 = 0, and Pn:n−1 = 1). Thus, applying (2.3) and (2.4), respectively, to the probability in the first term and the conditional probabilities inside the summation in the second term of (2.2) and going back to (2.1), we have the following lemma providing an explicit expression for the FDR of a two-stage stepup procedure under any distributional setting. Lemma 2.1. The FDR of a two-stage stepup procedure with the first-stage critical values 0 = 0 < 1 · · · n < n+1 = 1 and the second-stage critical values 0 = j 0 < j 1 · · · jj , with jj j , j = 1, . . . , n, is given by FDR =
n 1 (−i) (−i) (−i) Pr{Pj −1:n−1 j , Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n , Pi jj } j
i∈I0 j =1
+
j −1 n i∈I0 j =2 k=1
(−i)
(−i)
(−i)
E Pr{Pk:n−1 > j,k+1 , . . . , Pj −2,n−1 > j,j −1 , jj Pj −1:n−1 < j ,
(−i) (−i) Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n
I (Pi j k ) I (Pi j,k+1 ) − | Pi } k k+1
.
(2.5)
Remark 2.2. Lemma 2.1 generalizes the formula for the FDR given in Lemma 3.2 (with r = 1) of Sarkar (2002) from a single-stage stepup to a two-stage stepup procedure. It reduces to the formula for a single-stage stepup procedure with the critical values 0 = 0 < 1 · · · n < n+1 = 1, when j k = j , 1 k j n. 3. FDR-controlling two-stage stepup procedures under independence We determine in this section the first-stage critical values j , j = 1, . . . , n, and the second-stage critical values j k , 1 k j n, in the above two-stage stepup procedure providing a control of the FDR at ∈ (0, 1), assuming that the Pi ’s are independent and, when the null hypotheses are true, they are distributed as U(0, 1).
1076
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
Theorem 3.1. Consider a two-stage stepup procedure with the second-stage critical values j k ’s satisfying j k = min{j , kj }, for k = 1, . . . , j , given the first-stage critical values j ’s and some other constants j ’s. The FDR of this procedure under the independence is given by FDR =
n
min
i∈I0 j =1
j (−i) (−i) (−i) , j Pr{Pj −1:n−1 j , Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n }. j
(3.1)
Proof. The FDR of a two-stage stepup procedure, under independence, as we see from Lemma 2.1, is FDR =
n jj i∈I0 j =1
+
j
(−i)
(−i)
(−i)
Pr{Pj −1:n−1 j , Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n }
j −1 n i∈I0 j =2 k=1
(−i)
(−i)
(−i)
Pr{Pk:n−1 > j,k+1 , . . . , Pj −2,n−1 > j,j −1 , jj Pj −1:n−1 < j , (−i)
(−i)
Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n }
j,k+1 j k − . k k+1
(3.2)
For every fixed j = 2, . . . , n in the above triple summation, consider the two possible situations: (i) j j j and (ii) j j > j . In situation (i), since kj j j j , we have j k = min{j , kj } = kj , for all k = 1, . . . , j . Whereas, in situation (ii), we have j = min{j , j j } = jj , whatever be k. In each of these situations, the summation over k = 1 to j − 1 for that particular j will be zero. In other words, the expression (3.2) reduces to the right-hand side of (3.1) if we choose the j k ’s as stated in the theorem. Thus, the theorem is proved. Remark 3.1. Theorem 3.1 with the j ’s and j ’s satisfying the inequality FDR provides a class of two-stage stepup FDR procedures under the independence. The single-stage BH procedure belongs to this class. To see this, let j =j /n and choose any kj j /n, k = 1, . . . , j . Since min{j /j, j } = /n in this case, the FDR is equal to n0 /n, as n j =1
(−i)
(−i)
(−i)
Pr{Pj −1:n−1 j , Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n } = 1,
(3.3)
for any 1 · · · n . However, as the j k ’s satisfy j k =j /n, 1 k j n, this is no different from the BH procedure (see Remark 2.1). Remark 3.2. To obtain a proper two-stage FDR procedure, first note that if we let j = (1 − j )/(n − j + 1), j = 1, . . . , n, then the FDR in (3.1) becomes less than or equal to
n 1 − j (−i) (−i) (−i) Pr{Pj −1:n−1 j , Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n } n−j +1
i∈I0 j =1
= n0 (1 − n ) +
n−1 1 − j (−i) (−i) Pr{Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n } n−j +1
i∈I0 j =1
−
n 1 − j (−i) (−i) Pr{Pj −1:n−1 > j , . . . , Pn−1:n−1 > n } n−j +1
i∈I0 j =2
= n0 (1 − n ) +
n−1 i∈I0 j =1
(−i)
(−i)
Pr{Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n }
1 − j +1 1 − j − . n−j +1 n−j
(3.4)
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
1077
We choose the j ’s as follows. For any fixed 1 m n and 0 < < 1, let 1 − j = (n − j + 1)(1 − )/n for j = 1, . . . , m, and j = m for j = m + 1, . . . , n. For this choice of the j ’s, the expression (3.4) reduces to n0 (1 − m ) − (1 − m )
n−1 Pr{P (−i) > } m j :n−1 i∈I0 j =m
(n − j )(n − j + 1)
.
(3.5)
We now consider the following: Mixture model: With i ≡ (Hi ) defined to be 0 or 1 according as Hi is true or false, let (Pi , i ), i = 1, . . . , n, be iid with Pr{Pi u|i } = (1 − i )F0 (u) + i F1 (u), i ∼ Bernoulli(1 − 0 ),
(3.6)
where F0 (u) = u and F1 is stochastically smaller than F0 . Under this model, the Pi ’s are now iid with F (u) = 0 u + (1 − 0 )F1 (u). Therefore, (3.5) can be expressed as ⎡ ⎤ n−1 Pr{Pj :n−1 > m } ⎦ n0 (1 − m ) ⎣1 − (n − j )(n − j + 1) j =m
⎡
= n0 (1 − m ) ⎣1 − ⎡ n0 (1 − m ) ⎣1 −
⎤ Pr{Pn−j :n−1 > m } ⎦ j (j + 1)
n−m j =1
n−m j =1
⎤
Pr{Pn−m+1−j :n−m > m } ⎦ j (j + 1)
⎡ ⎤ n−m Pr{Pn−m+1−j :n−m > m } 0 (1 − m ) ⎦ = n {1 − F (m )} ⎣1 − 1 − F (m ) j (j + 1) j =1
n0 }, {1 − n−m+1 m n−m+1
(3.7)
where P1:n−m . . . Pn−m:n−m are the ordered components of a subset of n − m of the original n p-values. The first inequality in (3.7) follows from the fact that, since the jth ordered component of a set of n values is more than the (j −1)th ordered component of any of its subset of n−1 values, we have Pn−j :n−1 Pn−j −1:n−2 · · · Pn−m+1−j :n−m . Whereas, the second inequality follows, first from that fact that 0 (1 − m ) 0 (1 − m ) = 1, 1 − F (m ) 0 (1 − m ) + (1 − 0 )(1 − F1 (m )) and second knowing that ⎡ (n − m + 1){1 − F (m )} ⎣1 −
n−m j =1
⎤ Pr{Pn−m+1−j :n−m > m } ⎦ = Pr{Pn−m+1:n−m+1 > m }, j (j + 1)
(3.8)
from Sarkar (1998, p. 499), which is less than or equal to Pr F0 {Pn−m+1:n−m+1 > m }, as the Pi ’s are stochastically smaller under F than under F0 . If we choose 0 < < 1 and 1 m n is such a way that m − 1 1/(n−m+1) n−m+1 m = 1 − (1 − ) , n n
1078
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
that is,
((m − 1)/n)1/(n−m+1) − (m − 1)/n , 1 − (m − 1)/n
(3.9)
then the FDR in (3.7) is controlled at 0 . Thus, we have the following theorem. Theorem 3.2. The following class of two-stage stepup procedures controls the FDR at under the mixture model (3.6): (i) At stage one, consider the set of constants j = 1 − ((n − j + 1)/n)(1 − ), j = 1, . . . , m, and j = m , j = m + 1, . . . , n, for any fixed 0 < < 1 and 1 m n satisfying (3.9). Find j = max{0 i n: Pi:n i }. (ii) Accept Hi for all i > j and go to the second stage. (iii) At stage two, determine the constants 0 = j 0 j 1 · · · jj satisfying j l = min{j , l(1 − j )/(n − j + 1)}, l = 0, 1, . . . , j , and find k = max{0 l j : Pl:n j l }. (iv) Reject Hi for all i k and accept the rest. Remark 3.3. The FDR of a procedure in Theorem 3.2 under the mixture model, as we see from (3.1), is given by
FDR =
n i∈I0 j =1
min
j (1 − j ) (−i) (−i) (−i) , Pr{Pj −1:n−1 j , Pj :n−1 > j +1 , . . . , Pn−1:n−1 > n }. j n−j +1
(3.10)
Theorem 3.2 provides a newer class of FDR procedures under the mixture model modifying the BH procedure in the spirit of Benjamini et al. (2006) and Storey et al. (2004). Each modifies the critical values of the BH procedure, from the j k = k/n to those satisfying j k = min{j , k/nˆ 0 (j )} at the second stage using an estimate nˆ 0 (j ) of n0 based on the number of hypotheses accepted at the first stage. The rationale behind such a modification, as given in Benjamini and Hochberg (2000), Benjamini et al. (2006) and Storey et al. (2004), is that if in the BH procedure is replaced by min{1, n/n0 }, the FDR, at least under the independence, can be controlled less conservatively. The idea of estimating n0 based on the number of insignificant hypotheses, first outlined in Schweder and Spjotvoll (1982), was used in the aforementioned papers and Sarkar (2006). Notice that we have chosen the estimate nˆ 0 (j ) = (n − j + 1)/(1 − j ), with n − j + 1 in the numerator, and restricted the second-stage critical values to the region found significant at the first stage, as in Storey et al. (2004) and Sarkar (2006). This is in contrast with Benjamini et al. (2006) who chose n − j , the number of hypotheses accepted at the first stage, and allowed the second-stage critical values to fall in the insignificant region found at the first stage. Also notice that in the above theorem j k k for all 1 j k n. So, according to Proposition 2.1, it does not make any difference whether or not we make any decision at stage one in the above procedure unless all the hypotheses are accepted. In other words, the procedures in Theorem 3.2 can be described equivalently as in the following corollary. Corollary 3.1. The following class of two-stage stepup procedures controls the FDR at under the mixture model: At stage one, consider the set of constants as in Theorem 3.2 and find j = min{0 i n: Pi:n i }. If j = 0, accept all hypotheses and stop; otherwise, go to stage two and apply the stepup procedure to all the hypotheses using the modified BH critical values j k = min{j , k(1 − j )/(n − j + 1)}, k = 1, . . . , n. The above theorem or its corollary provides a number of two-stage modified BH procedures, for different choices of m, each controlling the FDR under the mixture model. It is important to note, however, that for m = 1 and n, we do not need to assume the mixture model; just independence under any fixed configuration, yet unknown, of true and false null hypotheses will suffice. More precisely, when m = n, since 1 − j = (n − j + 1)(1 − )/n, i = 1, . . . , n,
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
1079
the FDR in (3.4) equals n0 (1 − )/n, for any 0 < < 1. When m = 1, (3.4) is equal to n0 (1 − ) − (1 − )
(−i)
n−1 i∈I0 j =1
Pr{Pj :n−1 > } (n − j )(n − j + 1) (−i)
n−1 Pr{Pj :n−1 } n0 = (1 − ) + (1 − ) n (n − j )(n − j + 1) i∈I0 j =1
(1 − ) + (1 − )
(−i)
n n−1 i=1 j =1
Pr{Pj :n−1 } (n − j )(n − j + 1)
= Pr{Pn:n > },
(3.11)
which is less than or equal to for any 0 < < 1. The last line in (3.11) follows from Sarkar (1998). In fact, our procedure with m = 1, that is, the two stage procedure with j k ’s given by k(1 − ) j k = min , , k = 1, . . . , n, n−j +1 for any arbitrary 0 < < 1, is same as the one developed by Storey et al. (2004). The FDR of this procedure under independence, as seen from (3.1), is given by n (1 − ) (−i) (−i) FDR = min (3.12) , Pr{Pj −1:n−1 < Pj :n−1 }. j n−j +1 i∈I0 j =1
which is much more explicit than given in Storey et al. (2004), and we have given here an alternative proof that it controls the FDR. Some other interesting results related to this procedure are worth noting. For instance, FDR
n (−i) (−i) Pr{Pj −1:n−1 < Pj :n−1 } j
i∈I0 j =1
= FDR1 ,
(3.13)
with the equality holding when /(n + ), where FDR1 is the FDR of the first-stage single-step procedure. Let FNR1 and pFNR1 , respectively, be the FNR and pFNR of the this single-step procedure. Then, we also have FDR
n i∈I0 j =1
1− (−i) (−i) Pr{Pj −1:n−1 < Pj :n−1 } n−j +1
⎡
= ⎣Pr{Pn:n > } −
n i∈I \I0 j =1
⎤ 1 (−i) (−i) P {Pj −1:n−1 < Pj :n−1 , Pi > }⎦ n−j +1
= [Pr{Pn:n > } − FNR1 ] = Pr{Pn:n > }(1 − pFNR1 ),
(3.14)
with the equality holding when n/(1 + n). See Sarkar (2006) for formulas of FDR and FNR of a single-step procedure. Again, the inequality (3.14) provides an alternative proof of the FDR controlling property of this procedure under the independence. Storey et al. (2004) used the inequality FDR Pr{max1 i n0 Pi > } to prove this property. 4. Numerical study Given 1000 independent random variables Xi ∼ N(i , 1), i = 1, . . . , 1000, we consider simultaneous testing of the null hypotheses Hi : i = 0, i = 1, . . . , 1000, against their respective alternatives Ki : i > 0, i = 1, . . . , 1000. Among
1080
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
Adaptive BH BH
Sarkar m=1 0.2
δ=1
0.4
0.6
0.8
1.0
δ=2
0.05
FDR
0.04
0.03
0.02
0.01
0.2
0.4
0.6
0.8
1.0 π0
Fig. 1. Comparison of different procedures in terms of the FDR.
several possible choices within the proposed class of FDR procedures, we consider the one with m = 10 and = 0.9953 and numerically compare it with the STS procedure (m = 1) with = 0.5, the adaptive BH procedure and the original BH procedure. Notice that the choice of in our procedure is dictated by the constrained in (3.9) given m = 10. We examine the amount of improvement our procedure (lebeled Sarkar) offers over the BH procedure relative to the STS and the adaptive BH procedures. We also compare these procedures in terms of the average power, defined as the expected proportion false null hypotheses that are rejected (Dudoit et al., 2003; Shaffer, 2002). Having generated values of n = 1000 independent random variables Xi ∼ N(i , 1), i = 1, . . . , 1000, we perform 1000 Z tests of = 0 against > 0 using each of the four procedures at the FDR level = 0.05. We randomly identify each hypothesis to be true ( = 0) or false ( = ), for some fixed > 0, according to a Bernoulli model with 0 = 1 − 1 as the probability of a null hypothesis being true. The values of Q = V /R, the ratio of the number of rejected true null hypotheses (V) to the total number of rejected null hypotheses (R), and S/n1 , the proportion of the number of rejected false null hypotheses (S) to the total number of false null hypotheses (n1 ), are then calculated for each procedure. Repeating this 30,000 times, we estimate the FDR by averaging the 30,000 Q values and the average power by averaging the 30,000 S/n1 -values. The value of Q is considered to be zero if R = 0 in a particular repetition. Whereas, since n1 is random in a mixture model and can be zero with positive probability (even though it is extremely small), we also consider S/n1 = 0 if n1 = 0. Figs. 1 and 2 compare the FDR and the average power, respectively, of these procedures for = 1 and 2. As seen from Fig. 1, for improving the FDR control of the BH procedure, ours seems to have the best performance when in a false null hypothesis the value of the mean departs slightly from its specified null hypothetical value; whereas, the STS procedure seems to work the best when this departure is large. Similar picture is seen in Fig. 2 that compares the average power, even though, for each of these procedures, the power of detecting small departures of the means from their null hypothetical values is generally low, could be even lower than unless many of the null hypotheses are false.
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
Adaptive BH
1081
Sarkar
BH
m=1 0.2
δ=1
0.4
0.6
0.8
1.0
δ=2
1.0
AvePower
0.8
0.6
0.4
0.2
0.0 0.2
0.4
0.6
0.8
1.0 π1
Fig. 2. Comparison of different procedures in terms of the average power.
5. FDR-controlling two-stage stepup procedures under dependence In this section, we discuss how to determine critical values of a two-stage stepup procedure providing a control of the FDR under any form of dependence among the test statistics. To this end, we first have the following lemma. Lemma 5.1. The FDR of a two-stage stepup procedure satisfies the following inequality: FDR
j n j,k − j,k−1 . k
(5.1)
i∈I0 j =1 k=1
Proof. Note that in the expression (2.1) for the FDR of a two-stage stepup procedure, j 1 k=1
k
P {J = j, K(j ) = k, Pi j k } =
j k 1 k=1
=
k
j j 1 l=1 k=l
k
j j 1 l=1
l
P {J = j, K(j ) = k, j,l−1 < Pi j l }
l=1
k=1
P {J = j, K(j ) = k, j,l−1 Pi < j,l }
P {J = j, K(j ) = k, j,l−1 Pi < j l }
1082
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
j 1 l=1
l
P {J = j, j,l−1 Pi < j l }
j j l − j,l−1 . l
(5.2)
l=1
Using (5.2) in (2.1), we get the inequality (5.1).
Let us choose the j k ’s satisfying kj , j k = min j , j k=1 (1/k)
(5.3)
1k j n, for some j ’s, modifying the BY critical values at the second stage using the number of hypotheses accepted (which is n − j ) at the first stage and truncating each at j , as in the case of modifying the BH procedure (Section 3). Then, we have from (5.1) that FDR n0 nj=1 j , which can be made less than or equal to by choosing the j ’s in several different ways. For instance, one might choose the j ’s as follows: j = /{n(n − j + 1) nk=1 (1/k)}. Thus, we have the following as a two-stage version of the BY procedure: Procedure 5.1. At stage one, consider any set of constants 0 = 0 < 1 · · · n < n+1 = 1. At stage two, consider the critical values j k ’s satisfying k , k = 1, . . . , n. (5.4) j k = min j , j n(n − j + 1) k=1 k1 nk=1 k1 Considering Lemma 5.1 for a single-stage procedure, we notice that one can develop a single-stage procedure different from the BY procedure that will also control the FDR under any form dependence among the p-values. We present this single-stage procedure in the following before constructing a two-stage version of this. Lemma 5.2. Consider a single-stage stepup procedure with the critical values k ’s satisfying k = k(k + 1)/{2n2 }. It has the FDR less than or equal to n0 /n. Proof. From Lemma 5.1, we see that the FDR of this single-stage stepup procedure is given by FDR
n k − k−1 k
i∈I0 k=1
n 1 2 [k(k + 1) − (k − 1)k] 2n k i∈I0
n0 = . n
k=1
(5.5)
The following is a two-stage version of the alternative single-stage procedure in Lemma 5.2. Procedure 5.2. At stage one, consider any set of constants 0 = 0 < 1 · · · n < n+1 = 1. At stage two, consider the critical values dj k ’s satisfying k(k + 1) j k = min j , , k = 1, . . . , n. (5.6) 2nj (n − j + 1) nk=1 k1
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
BY
1083
Sarkar 200
ρ = 0.5
400
600
800
1000
ρ = 0.75 0.006
0.004
0.002
FDR
0.000 ρ=0
ρ = 0.25
0.006
0.004
0.002
0.000 200
400
600
800
1000 n0
Fig. 3. Comparison of the BY and the new procedures in terms of the FDR.
Before numerically investigating how Procedures 5.1 and 5.2 perform as two-stage modifications of their respective single-stage versions, we assess how good the new single-stage procedure proposed in Lemma 5.2 is compared to the well-known BY procedure. Fig. 3 compares the simulated vales of the FDR for the BY procedure and the new procedure (labeled Sarkar) in simultaneous testing of the means of 1000 normal random variables, each mean being tested at 0 against a positive value using the Z test. The correlations are assumed to be common, which is , and each variance is 1. The simulated values are based on 20,000 repetitions and the alternative mean in each test is chosen to be 2. Interestingly, the new procedure seems to perform much better than the BY procedure when the p-values are highly dependent and the number of true null hypotheses is neither extremely small nor extremely large. However, when we brought Procedures 5.1 and 5.2 into our comparisons, we noticed, unfortunately, that none of these two-stage procedures provide any significant improvement over the BY procedure or its alternative in terms of the FDR control. Thus, it appears that the way we estimate n0 to improve the BH procedure under the independence may not work under dependence for improving the BY or its alternative procedure. 6. Concluding remarks Starting with a general two-stage stepup procedure, we have made an attempt in this article to present a wider class of FDR procedures. In particular, we have offered a theory of modifying the BH procedure through a two-stage approach incorporating an estimate of the number of true null hypotheses and maintaining the FDR control, which is more general than Storey et al. (2004). Although our general theory relies on the iid mixture model for the p-values, for the STS procedure, it not only reduces to the one requiring only the independence condition of the p-values, as
1084
S.K. Sarkar / Journal of Statistical Planning and Inference 138 (2008) 1072 – 1084
they have assumed, but also it provides an alternative proof of its FDR control property. We have chosen only one from our proposed class of procedures and shown numerically that it works quite well in terms of improving the FDR control of the BH procedure. We believe that there are other procedures in this class that would also perform well and more studies are needed to investigate it. Another positive contribution we have made in this article is that we have introduced a new single-stage stepup procedure that might work better than the commonly used BY procedure in terms of controlling the FDR under any form of dependence among the p-values in some instances. Acknowledgment The author is grateful to an Associate Editor and two referees for making useful and constructive comments, and to Wenge Guo for giving some additional comments and pointing out an error in an earlier version of the manuscript, all of which led to a much improved form of the paper. He also thanks Zijiang Yang for doing the numerical calculations. Appendix Proof of Proposition 2.1. The two events {J = j, K ∗ (j ) = j } and {J = j, K(j ) = j } can be seen to be identical if j k k , for 1j k n. This is because, for 1 k < j n, {J = j, K ∗ (j ) = k} = {Pj :n j , Pj +1:n > j +1 , . . . , Pn:n > n } ∩ {Pk:n j k , Pk+1,n > j,k+1 , . . . , Pn:n > j,n } = {Pk:n j k , Pk+1,n > j,k+1 , . . . , jj < Pj :n j , Pj +1:n > j +1 , . . . , Pn:n > n },
(A.1)
which is same as {J = j, K(j ) = k}. For 1 k = j n, {J = j, K ∗ (j ) = j } = {Pj :n jj , Pj +1:n > j +1 , . . . , Pn:n > n } = {J = j, K(j ) = j }. And, for 1j < k n, both of these events are null. This proves the proposition.
(A.2)
References Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. Roy. Statist. Soc. Ser. B 57, 289–300. Benjamini, Y., Hochberg, Y., 2000. On the adaptive control of the false discovery rate in multiple testing with independent statistics. J. Educ. Behav. Statist. 25, 60–83. Benjamini, Y., Yekutieli, D., 2001. The control of the false discovery rate in multiple testing under dependence. Ann. Statist. 29, 1165–1188. Benjamini, Y., Krieger, A.M., Yekutieli, D., 2006. Adaptive linear step-up procedures that control the false discovery rate. Biometrika 93, 491–507. Dudoit, S., Shaffer, J.P., Boldrick, J.C., 2003. Multiple hypothesis testing in microarray experiments. Statist. Sci. 18, 71–103. Finner, H., Roters, M., 2001. On the false discovery rate and expected number of type I errors. Biometrical J. 43, 985–1005. Genovese, C., Wasserman, L., 2002. Operarting characteristics and extensions of the false discovery rate procedure. J. Roy. Statist. Soc. Ser. B 64, 499–517. Hochberg, Y., Benjamini, Y., 1990. More powerful procedures for multiple significance testing. Statist. Med. 9, 811–818. Sarkar, S.K., 1998. Some probability inequalities for ordered MTP2 random variables: a proof of the simes conjecture. Ann. Statist. 26, 494–504. Sarkar, S.K., 2002. Some results on false discovery rate in stepwise multiple testing procedures. Ann. Statist. 30, 239–257. Sarkar, S.K., 2006. False discovery and false non-discovery rates in single-step multiple testing procedures. Ann. Statist. 34, 394–415. Schweder, T., SpjZtvoll, E., 1982. Plots of p-values to evaluate many tests simultaneously. Biometrika 69, 493–502. Shaffer, J.P., 2002. Multiplicity, directional (type III) errors, and the null hypothesis. Psychol. Methods 7, 356–369. Storey, J.D., 2002. A direct approach to false discovery rates. J. Roy. Statist. Soc. Ser. B 64, 479–498. Storey, J.D., Taylor, J.E., Siegmund, D., 2004. Strong control, conservative point estimation and simultanaeous conservative consistency of false discovery rates: a unified approach. J. Roy. Statist. Soc. Ser. B 66, 187–205. Zhang, Z., 2005. Multiple hypothesis testing for finite and infinite number of hypotheses. Ph.D. Thesis, Case Western Reserve University.