R如何计算错误发现率

How Does R Calculate the False Discovery Rate

当我使用 R 的 p.adjust 函数计算错误发现率时,我似乎得到了不一致的结果。基于 documentation 中引用的论文 调整后的 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 对可以解释这两个结果的数据做了什么?

原因是 FDR 计算确保 FDR 永远不会随着 p 值的减小而增加。这是因为您始终可以选择为拒绝规则设置更高的阈值,前提是更高的阈值会使您的 FDR 更低。

在您的案例中,您的第二个假设的 p 值为 0.0006,FDR 为 0.010971585,但第三个假设的 p 值更大,FDR 更小。如果您可以通过将 p 值阈值设置为 0.0019 来实现 0.009120228 的 FDR,那么没有理由仅仅为了获得更高的 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的顺序倒退。

您可以在您 link 的 Benjamini-Hochberg paper 中看到这一点,包括在第 293 页的过程定义中,其中指出(强调我的):

let k be the largest i for which P(i) <= i / m q*;

then reject all H_(i) i = 1, 2, ..., k