R 3.1.1中p值的局部FDR

Question

我正在试图找到一个包来计算Efron的本地FDR来进行一系列测试。我有超过1000个协变量，所以多次修正是肯定的。

查找本地FDR将看到locfdr包在CRAN上不再可用。知道为什么会被移除吗？这似乎与当地FDR的原始出版物有着最密切的关系。

我确实找到了fdrtool，但它不能从p值计算本地FDR。我发现的其他软件包在3.1.1-LocalFDR、LocalFDR、kerfdr、twilight中是不可用的。当然，所有这些包使用的方法略有不同。即使我能找到他们，该选哪一个呢？

Answer 1

twilight包可以这样做(正如注释中所建议的那样)。它很容易使用，如下图所示，生物导体上也有相当清晰的微缩体。

首先，我们安装软件包。

# Install twilight
source("http://bioconductor.org/biocLite.R")
biocLite("twilight")

接下来，我们模拟了一些测试分数和p值.

# Simulate p p-values
set.seed(1)
p <- 10000 # Number of "genes"
prob <- 0.2 # Proportion of true alternatives

# Simulate draws from null (=0) or alternatives (=1)
null <- rbinom(p, size = 1, prob = prob) 
# Simulate some t-scores, all non-null genes have an effect of 2
t.val <- (1-null)*rt(p, df = 15) + null*rt(p, df = 15, ncp = 2)
# Compute p-values
p.val <- 2*pt(-abs(t.val), df = 15)

# Plot the results
par(mfrow = c(1,2))
hist(t.val, breaks = 70, col = "grey",
     xlim = c(-10, 10), prob = TRUE, ylim = c(0, .35))
hist(p.val, breaks = 70, prob = TRUE)

​

​

接下来，我们加载库并运行fdr分析：

library(twilight)
ans <- twilight(p.val)

我们看到，对实际替代比例的估计是相当好的：

print(1 - ans$pi0)
#[1] 0.1529
print(prob)
#[1] 0.2

包按递增顺序重新排序p值、q值和fdr值，所以我们做了一个技巧来重建原来的顺序。

o <- order(order(p.val)) 
fdr <- ans$result$fdr[o]
plot(p.val, fdr, pch = 16, col = "red", cex = .2) 

​

​

最后，我们可以把重要的与真实的交叉起来：

table(estimate = fdr < 0.5, truth = as.logical(null))
#        truth
#estimate FALSE TRUE
#   FALSE  7564 1172
#   TRUE    368  896

因此，我们有一个84.6 %的准确性在这个玩具-例子。我希望这能帮到你。twilight函数还为fdr提供了引导CIs功能，您可以在?twilight中找到它，并提供进一步的参考资料。

编辑看起来fdrtool包(它在CRAN上)实际上可以从p值计算本地fdr。在我们的例子中，我们做了以下工作：

library("fdrtool")
fdr <- fdrtool(p.val, statistic="pvalue")
fdr$lfdr # The local fdr values