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
当我使用 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