47
$\begingroup$

My statistics program implements both the Benjamini & Hochberg (1995) and Benjamini & Yekutieli (2001) false discovery rate (FDR) procedures. I have done my best to read through the later paper, but it is quite mathematically dense and I am not reasonably certain I understand the difference between the procedures. I can see from the underlying code in my statistics program that they are indeed different and that the latter includes a quantity q that I have seen referred to in regards to FDR, but also don't quite have a grasp of.

Is there any reason to prefer the Benjamini & Hochberg (1995) procedure versus the Benjamini & Yekutieli (2001) procedure? Do they have different assumptions? What are the practical differences between these approaches?

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B, 57, 289–300.

Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 29, 1165–1188.

The 1999 paper referenced in the comments below: Yekutieli, D., & Benjamini, Y. (1999). Resampling-based false discovery rate controlling multiple test procedures for correlated test statistics. Journal of Statistical Planning and Inference, 82(1), 171-196.

Improve this question
$\endgroup$
4
  • $\begingroup$ I thought the 2001 paper establishes properties of FDR (1995) under dependence. Yekutieli and Benjamini (Journal of Statistical Planning and Inference, 1999) establishes a different FDR procedure. Any chance that is the one you are seeking? $\endgroup$
    – julieth
    Jul 6, 2013 at 1:22
  • $\begingroup$ @julieth: That was my sense of the 2001 paper upon reading the abstract alone, but the formulas in the paper (e.g. 27 - 30) do seem to involve a quantity referred to as q. Then again, so does this 1999 paper you cite. My sense though is that the 1999 paper implements a re-sampling approach which is clearly (from looking at the code) is not what my stats program is doing (R; p.adjust)... but I may be wrong. $\endgroup$
    – russellpierce
    Jul 6, 2013 at 2:43
  • 1
    $\begingroup$ At the end of the 2001 paper the 1999 paper is cited and they say "Finally, recall the resampling based procedure of Yekutieli and Benjamini (1999), which tries to cope with the above problem and at the same time utilize the information about the dependency structure derived from the sample. The resampling based procedure is more powerful, at the expense of greater complexity and only approximate FDR control." ... so I think the 2001 paper provided a closed form computational solution and that is what my stats program is implementing. $\endgroup$
    – russellpierce
    Jul 6, 2013 at 2:46
  • 1
    $\begingroup$ Ok, so you are using p.adjust. The 99 paper is different altogether as you noted. I always saw the BY option in p.adjust and didn't pay attention. That 2001 paper is usually cited in relation to the proof of FDR and 'positive regression dependence'. I never saw it as quoting a different estimator, but maybe it is in there. Looks like I need to reread it. $\endgroup$
    – julieth
    Jul 6, 2013 at 2:52

2 Answers 2

Reset to default
28
$\begingroup$

Benjamini and Hochberg (1995) introduced the false discovery rate. Benjamini and Yekutieli (2001) proved that the estimator is valid under some forms of dependence. Dependence can arise as follows. Consider the continuous variable used in a t-test and another variable correlated with it; for example, testing if BMI differs in two groups and if waist circumference differs in these two groups. Because these variables are correlated, the resulting p-values will also be correlated. Yekutieli and Benjamini (1999) developed another FDR controlling procedure, which can be used under general dependence by resampling the null distribution. Because the comparison is with respect to the null permutation distribution, as the total number of true positives increases, the method becomes more conservative. It turns out that BH 1995 is also conservative as the number of true positives increases. To improve this, Benjamini and Hochberg (2000) introduced the adaptive FDR procedure. This required estimation of a parameter, the null proportion, which is also used in Storey's pFDR estimator. Storey gives comparisons and argues that his method is more powerful and emphasizes the conservative nature of 1995 procedure. Storey also has results and simulations under dependence.

All of the above tests are valid under independence. The question is what kind of departure from independence can these estimates deal with.

My current thinking is that if you don't expect too many true positives the BY (1999) procedure is nice because it incorporates distributional features and dependence. However, I'm unaware of an implementation. Storey's method was designed for many true positives with some dependence. BH 1995 offers an alternative to the family-wise error rate and it is still conservative.

Benjamini, Y and Y Hochberg. On the Adaptive Control of the False Discovery Rate in Multiple Testing with Independent Statistics. Journal of Educational and Behavioral Statistics, 2000.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Thanks a lot! Could you revise your question to clarify the following points/issues: "resampling the null distribution" is the 1999 paper? Would you please provide the citation for the 2000 paper? You seemed familiar with p.adjust, is it actually implementing the BY procedure? Must one use BH when the hypothesis tests aren't dependent? What causes hypothesis tests to be considered dependent? - Please let me know if any of these questions goes beyond the present scope and requires a new question be asked. $\endgroup$
    – russellpierce
    Jul 6, 2013 at 12:17
  • $\begingroup$ p.adjust has options for both (BH and BY). However, I thought these were the same, so I missed something. $\endgroup$
    – julieth
    Jul 6, 2013 at 14:12
  • $\begingroup$ And the underlying code is different too (I checked) so they will produce different numbers. $\endgroup$
    – russellpierce
    Jul 6, 2013 at 15:39
  • $\begingroup$ So what procedure is it that you think p.adjust is performing with the BY argument? I don't think it is the 1999 procedure. The underlying code is pmin(1, cummin(q * n/i * p[o]))[ro]. BH is pmin(1, cummin(n/i * p[o]))[ro]. So they only differ in q which is sum(1/(1:n)) where n = the number of pvalues. o and ro just serve to put the p values in decreasing numeric order for the function and then spit them back out in the same order the user inputed them in. $\endgroup$
    – russellpierce
    Jul 7, 2013 at 1:09
  • 1
    $\begingroup$ So, since no new answers are coming in, I'll accept this answer and summarize my understanding. p.adjust may be misciting for BY. What is performed is not resampling. BH, 2000 introduced the adaptive FDR procedure, and this involves the estimation of the null proportion, which may be the q that appears in the BY code. In the interim, it seems the sensible thing to do is to cite p.adjust directly as that reflects the actual procedure used when you use the option "BY" and to just be aware that "BY" may actually be implementing Benjamini & Hochberg, 2000. $\endgroup$
    – russellpierce
    Aug 28, 2013 at 16:52
7
$\begingroup$

p.adjust is not misciting for BY. The reference is to Theorem 1.3 (proof in Section 5 on p.1182) in the paper:

Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 29, 1165–1188.

As this paper discusses several different adjustments, the reference on the help page (at the time of writing) for p.adjust() is somewhat obscure. The method is guaranteed to control FDR, at the stated rate, under the most general dependence structure. There are informative comments in Christopher Genovese's slides at: www.stat.cmu.edu/~genovese/talks/hannover1-04.pdf Note the comment on slide 37, referring to the method of Theorem 1.3 in the BY 2001 paper [method='BY' with p.adjust()] that: "Unfortunately, this is typically very conservative, sometimes even more so than Bonferroni."

Numerical example: method='BY' vs method='BH'

The following compares method='BY' with method='BH', using R's p.adjust() function, for the p-values from column 2 of Table 2 in the Benjamini and Hochberg (2000) paper:

> p <-    c(0.85628,0.60282,0.44008,0.41998,0.3864,0.3689,0.31162,0.23522,0.20964,
0.19388,0.15872,0.14374,0.10026,0.08226,0.07912,0.0659,0.05802,0.05572,
0.0549,0.04678,0.0465,0.04104,0.02036,0.00964,0.00904,0.00748,0.00404,
0.00282,0.002,0.0018,2e-05,2e-05,2e-05,0)
> pmat <- rbind(p,p.adjust(p, method='BH'),p.adjust(p, method='BY'))
> rownames(pmat)<-c("pval","adj='BH","adj='BY'")
> round(pmat,4)


 
           [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
pval     0.8563 0.6028 0.4401 0.4200 0.3864 0.3689 0.3116 0.2352 0.2096
adj='BH  0.8563 0.6211 0.4676 0.4606 0.4379 0.4325 0.3784 0.2962 0.2741
adj='BY' 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
          [,10]  [,11]  [,12]  [,13]  [,14]  [,15]  [,16]  [,17]  [,18]
pval     0.1939 0.1587 0.1437 0.1003 0.0823 0.0791 0.0659 0.0580 0.0557
adj='BH  0.2637 0.2249 0.2125 0.1549 0.1332 0.1332 0.1179 0.1096 0.1096
adj='BY' 1.0000 0.9260 0.8751 0.6381 0.5485 0.5485 0.4856 0.4513 0.4513
          [,19]  [,20]  [,21]  [,22]  [,23]  [,24]  [,25]  [,26]  [,27]
pval     0.0549 0.0468 0.0465 0.0410 0.0204 0.0096 0.0090 0.0075 0.0040
adj='BH  0.1096 0.1060 0.1060 0.1060 0.0577 0.0298 0.0298 0.0283 0.0172
adj='BY' 0.4513 0.4367 0.4367 0.4367 0.2376 0.1227 0.1227 0.1164 0.0707
          [,28]  [,29]  [,30] [,31] [,32] [,33] [,34]
pval     0.0028 0.0020 0.0018 0e+00 0e+00 0e+00     0
adj='BH  0.0137 0.0113 0.0113 2e-04 2e-04 2e-04     0
adj='BY' 0.0564 0.0467 0.0467 7e-04 7e-04 7e-04     0

Note: The multiplier that relates the BY values to the BH values is $\sum_{i=1}^m (1/i)$, where $m$ is the number of p-values. Multipliers are, for example values m = 30, 34, 226, 1674, 12365:

> mult <- sapply(c(11, 30, 34, 226, 1674, 12365), function(i)sum(1/(1:i)))

setNames(mult, paste(c('m =',rep('',5)), c(11, 30, 34, 226, 1674, 12365)))

 m = 11     30     34    226   1674  12365 
  3.020  3.995  4.118  6.000  8.000 10.000

Check that for the example above, where $m$=34, the multiplier is 4.118

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.