Benjamini & Hochberg (1995) 和 Benjamini & Yekutieli (2001) 错误发现率程序之间的实际区别是什么?

机器算法验证 事后 错误发现率
2022-01-20 13:34:32

我的统计程序实现了 Benjamini & Hochberg (1995) 和 Benjamini & Yekutieli (2001) 错误发现率 (FDR) 程序。我已尽力通读后面的论文,但它在数学上相当密集,我不能合理地确定我理解这些程序之间的区别。我可以从我的统计程序中的底层代码中看到它们确实是不同的,并且后者包括一个我见过的关于 FDR 的数量 q,但也不太了解。

是否有任何理由更喜欢 Benjamini & Hochberg (1995) 程序而不是 Benjamini & Yekutieli (2001) 程序?他们有不同的假设吗?这些方法之间的实际区别是什么?

Benjamini, Y. 和 Hochberg, Y. (1995)。控制错误发现率:一种实用且强大的多重测试方法。皇家统计学会杂志 B 系列,57,289–300。

Benjamini, Y. 和 Yekutieli, D. (2001)。依赖关系下多次测试中错误发现率的控制。统计年鉴 29, 1165–1188。

以下评论中引用了 1999 年的论文:Yekutieli, D., & Benjamini, Y. (1999)。基于重采样的错误发现率控制相关测试统计的多个测试过程。统计规划与推理杂志,82(1),171-196。

2个回答

Benjamini and Hochberg (1995) 介绍了错误发现率。Benjamini and Yekutieli (2001) 证明了估计量在某些形式的依赖下是有效的。依赖可以如下产生。考虑 t 检验中使用的连续变量和与之相关的另一个变量;例如,测试两组的 BMI 是否不同,以及这两组的腰围是否不同。因为这些变量是相关的,所以得到的 p 值也将是相关的。Yekutieli 和 Benjamini (1999) 开发了另一种 FDR 控制程序,可以通过重新采样零分布在一般依赖下使用。因为比较是针对零排列分布的,所以随着真阳性总数的增加,该方法变得更加保守。事实证明,随着真阳性数量的增加,BH 1995 也是保守的。为了改进这一点,Benjamini 和 Hochberg (2000) 引入了自适应 FDR 程序。这需要估计一个参数,即空比例,该参数也用于 Storey 的 pFDR 估计器。Storey 进行了比较,并认为他的方法更强大,并强调了 1995 程序的保守性。Storey 也有依赖的结果和模拟。

上述所有测试在独立下都是有效的。问题是这些估计可以处理什么样的偏离独立性。

我目前的想法是,如果你不期望有太多的真阳性,那么 BY (1999) 程序会很好,因为它结合了分布特征和依赖性。但是,我不知道实施。Storey 的方法是为许多具有一定依赖性的真阳性而设计的。BH 1995 提供了家庭错误率的替代方案,它仍然是保守的。

本杰明尼,Y 和 Y 霍赫伯格。关于具有独立统计的多重测试中错误发现率的自适应控制。教育与行为统计杂志,2000 年。

p.adjust 没有混淆 BY。参考论文中的定理 1.3(第 1182 页第 5 节中的证明):

Benjamini, Y. 和 Yekutieli, D. (2001)。依赖关系下多次测试中错误发现率的控制。统计年鉴 29, 1165–1188。

由于本文讨论了几种不同的调整,p.adjust() 帮助页面上的参考(在撰写本文时)有些模糊。该方法保证在最一般的依赖结构下以规定的速率控制 FDR。Christopher Genovese 的幻灯片中有丰富的评论,网址为:www.stat.cmu.edu/~genovese/talks/hannover1-04.pdf 请注意幻灯片 37 上的评论,参考 BY 2001 论文中定理 1.3 的方法 [method= 'BY' 与 p.adjust()] 表示:“不幸的是,这通常非常保守,有时甚至比 Bonferroni 还要保守。”

数值示例: method='BY' vsmethod='BH'

下面使用 R 的 p.adjust() 函数对 Benjamini 和 Hochberg (2000) 论文中表 2 第 2 列的 p 值比较 method='BY' 和 method='BH':

> 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

注意:将 BY 值与 BH 值相关联的乘数是i=1m(1/i), 在哪里m是 p 值的数量。例如,乘数是 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

检查上面的示例,其中m=34,乘数为 4.118