R如何计算错误发现率答案

【问题标题】：How Does R Calculate the False Discovery RateR如何计算错误发现率
【发布时间】：2015-05-01 18:38:15
【问题描述】：

当我使用 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 对可以解释这两个结果的数据做了什么？

【问题讨论】：

我认为这是一个更好的问题，但它本质上是 stackoverflow.com/questions/10323817/… 的副本

标签： r fdr

【解决方案1】：

原因是 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 的顺序倒退。

您可以在链接到的Benjamini-Hochberg paper 中看到这一点，包括在第 293 页的过程定义中，其中指出（强调我的）：

让 k 成为最大的 i P(i)

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

【讨论】：