« SOP、SSOPのはんだ付け | Home | 箱庭 »
Rでfalse discovery rate (FDR)を計算する
By kmgs | 1 月 27, 2009
マイクロアレイ解析で多重比較の補正を行いたいとき、定番はFDRでしょう。Rでは、
> p <- runif(100, 0, 0.3) # p値もどき > p.adjust(p, "fdr") # Benjamini-Hochberg (BH)法 > p.adjust(p, "BH") # これもBH法
あえて自分で書けば
fdr <- function(p) {
# 大きい順でインデックスを取得
o <- rev(order(p))
# 比較回数による補正
n <- length(p) / (length(p):1)
# 補正の適用
f <- p[o] * n
# 比較したもののなかで、最小のFDRを保持
# また、FDRが1を超えないよう細工
val <- sapply(1:length(p), function(x) min(f[1:x], 1))
# 元のベクトルに格納 (名前の保持)
p[o] <- val
p
}みたいにして、
> sum(p.adjust(p, "fdr") == fdr(p)) # すべて一致することを確認 [1] 100
計算イメージは以下の通りです。
ほかにも、Benjamini-Yekutieli (BY)法や、ブートストラップ法でp値の分布を推定するqvalueパッケージ、fdrtoolパッケージなどといったものがあります。