R如何计算错误发现率

Question

当我使用R的p.adjust函数来计算错误发现率时,我似乎得到了不一致的结果.基于所引用的纸张文档经调整的p值应该这样来计算:

adjusted_p_at_index_i= p_at_index_i*(total_number_of_tests/i).

现在,当我跑步时,p.adjust(c(0.0001, 0.0004, 0.0019),"fdr")我得到了预期的结果

c(0.0003, 0.0006, 0.0019)

但是当我跑步时,p.adjust(c(0.517479039, 0.003657195, 0.006080152),"fdr")我得到了这个

c(0.517479039, 0.009120228, 0.009120228)

而不是我计算的结果:

c(0.517479039, 0.010971585, 0.009120228)

R对可以解释这两种结果的数据做了什么？

Answer 1

原因是FDR计算确保FDR不会随着p值的减小而增加.这是因为如果更高的阈值会使您的FDR更低,您总是可以选择为拒绝规则设置更高的阈值.

在您的情况下,您的第二个假设的p值为0.0006FDR 0.010971585,但第三个假设的p值较大且FDR较小.如果0.009120228通过将p值阈值设置为可以达到FDR 0.0019,则永远不会有理由设置较低的阈值以获得更高的FDR.

您可以在代码中输入p.adjust以下内容来查看:

...
}, BH = {
    i <- lp:1L
    o <- order(p, decreasing = TRUE)
    ro <- order(o)
    pmin(1, cummin(n/i * p[o]))[ro]

该cummin函数采用向量的累积最小值,按顺序向后p.

您可以在链接到的Benjamini-Hochberg论文中看到这一点,包括第293页的过程定义,其中说明了(强调我的):

令k 为 P(i)<= i/mq*的最大i ;

然后拒绝所有H_(i)i = 1,2,...,k