Rで解析：多重量検定のFDR調整比較「BonEV」パッケージ

多重比較で有意水準を調整していないt検定の結果で主張している内容を見かけます。FDRをBenjamini-Hochbergなどできちんと制御しましょう。

本パッケージではBenjamini-Hochberg,Storey_adjp,Bon-EVの方法が利用できます。コマンドのコードも合わせて紹介しますので参考にしてください。

パッケージバージョンは1.0。実行コマンドはR version 3.2.3で確認しています。

パッケージのインストール

下記、コマンドを実行してください。

#パッケージの利用には「qvalue」パッケージが必要です
install.packages("BiocManager")
BiocManager::install("qvalue")

#パッケージのインストール
install.packages("BonEV")

実行コマンド

詳細はコメント、パッケージのヘルプを確認してください。

#パッケージの読み込み
library("qvalue")
library("BonEV")

###仮想のp値を用意#####
set.seed(1111)
n <- 5
PValue <- sample(runif(n*10, min = 0, max = 1), n, replace = TRUE)
#####

#Benjamini-Hochberg,Storey_adjp,Bon-EVでFDRを制御:Bon_EVコマンド
Result <- Bon_EV(pvalue = PValue, alpha = 0.05)

#結果をdata.frame化
FDRResult <- data.frame(raw_P_value = Result$raw_P_value,
                        BH_adjp = Result$BH_adjp,
                        Storey_adjp = Result$Storey_adjp,
                        Bon_EV_adjp = Result$Bon_EV_adjp)
#表示
FDRResult
raw_P_value   BH_adjp Storey_adjp Bon_EV_adjp
1  0.79980081 0.7998008  0.11212878           1
2  0.07420536 0.1814442  0.02543773           1
3  0.12096279 0.1814442  0.02543773           1

#おまけ:Bon_EVコマンドのコード
#statパッケージ:p.adjustコマンド
#qvalueパッケージ:Storey JDらの方法
#論文:
#A direct approach to false discovery rates. Journal of the Royal Statistical Society, Series B, 64: 479-498.
#False discovery rates. In International Encyclopedia of Statistical Science. 
#http://genomine.org/papers/Storey_FDR_2011.pdf 
#Bon_EV_adjpはコードを参照
function (pvalue, alpha) 
{
  new_MTP_adjp <- rep(length(pvalue))
  BH <- p.adjust(pvalue, "BH")
  qobj <- qvalue(p = pvalue)
  qvalues <- qobj$qvalues
  ngene <- length(pvalue)
  for (i in 1:ngene) {
    adjpv <- ngene * (pi0est(pvalue)$pi0)/sum(BH &lt;= alpha, 
                                              na.rm = TRUE) * pvalue
    new_MTP_adjp[i] <- min(adjpv[i], 1)
  }
  mylist <- list(raw_P_value = pvalue, BH_adjp = BH, Storey_adjp = qvalues, 
                 Bon_EV_adjp = new_MTP_adjp)
  return(mylist)
}

少しでも、あなたの解析が楽になりますように！！