來自R中的p值的本地FDR 3.1.1

我試圖找到一個包來計算Efron的一系列測試的本地FDR。我有超過1000個協變量，所以多重修正肯定是有序的。來自R中的p值的本地FDR 3.1.1

尋找本地FDR是看到在CRAN上不再可用的locfdr軟件包。任何想法爲什麼它被刪除？這似乎與當地的FDR的原始出版物最密切相關。

我確實找到了fdrtool，但它無法從p值計算本地FDR。我發現的其他軟件包不適用於3.1.1 - LocalFDR，localFDR，kerfdr，twilight。當然，所有這些軟件包使用的方法稍有不同。即使我能找到他們，哪一個可以選擇？

來源

2014-10-06 KirkDCO

存檔版本都可以在這裏http://cran.r-project.org/src/contrib/Archive/locfdr/你可以下載和'install.packages（文件， repos = NULL，type =「source」）'。你可能想要調查爲什麼包被刪除。 – 2014-10-06 16:18:53

感謝您的提示。我確實得到了存檔版本，安裝了它，並且它正在工作。然而，它的日期是2011年2月，這讓我想知道爲什麼它過時了，爲什麼它被拉出來等等。還有一長串的其他軟件包不再在CRAN上 - nomi，LocalFDR，kerfdr， localFDR ......看起來奇怪的如此多的工具，做同樣的事情都丟失了。 – KirkDCO 2014-10-06 18:28:29

我對方法不熟悉，所以不能評論，但它似乎不是一個簡單的巧合。 – 2014-10-06 19:56:40

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)

Imgur

接下來我們加載庫和運行FDR分析：

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

我們看到，真正的替代比例的估計是相當不錯：

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

該軟件包將p值，q值和fdr值重新排序以增加orde r，所以我們做一些技巧來重建原來的訂單。

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

Imgur

最後，我們可以在交叉製表對顯著真：

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

因此，我們有84.6％在這個玩具，例如精確度。我希望這有幫助。 twilight函數還具有fdr的引導配置項，您可以在?twilight中找到它以及其他參考。

編輯看來，fdrtool包（這是CRAN）實際上可以計算從p值的本地FDR。在我們的例子中，我們做到以下幾點：

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

來源

2014-10-30 11:00:19

感謝您的評論和偉大的例子。這非常有幫助！ – KirkDCO 2014-10-30 15:06:09

我更感興趣「任何想法爲什麼它被刪除？」 – daj 2014-10-31 03:03:42

@daj沒有idaa。我找不到任何信息。但'fdrtool'軟件包也是一個選項（正如我編輯的答案中所見）。 – 2014-10-31 11:08:38

來自R中的p值的本地FDR 3.1.1

回答

相關問題