多重検定

数学いらずの医科統計学PART5 CHAPTER22の多重検定の話。
多重検定についてはぐぐると色々出てくる。
Rではp.adjustという関数でp値を各手法で補正してくれるようだ。とりあえず全部やってみる。

set.seed(123)
x <- rnorm(50, mean = c(rep(0, 25), rep(3, 25)))
p <- 2*pnorm(sort(-abs(x)))
m <- c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")

res <- NULL
for(i in seq(m)){
	res <- cbind(res, p.adjust(p, method=m[i]))
}
#順に並べた50個のp値がそれぞれどう変わるか
matplot(res, type="o", ylab="adjusted p value", pch=16, lty=1, col=rainbow(length(m)))
legend("bottomright", legend=m, bty="n", lty=1, col=rainbow(length(m)), lwd=3)
title(paste(length(p), "times tests"))
#もとのp値からどう変わるか
matplot(p, res, type="o", xlab="p value", ylab="adjusted p value", xlim=c(0, 1), pch=16, lty=1, col=rainbow(length(m)))
legend("bottomright", legend=m, bty="n", lty=1, col=rainbow(length(m)), lwd=3)
title(paste(length(p), "times tests"))

　
Benjamini-Hochberg法について話題になったので、これを確かめてみる。
p値を個計算するとき、p値は一様分布に従う。
BH法は得られたp値を小さい順に並べていって、番目のp値が、一様分布から期待されるp値()の分布と比べてどうかをみて、比べる。

k <- 50 #いくつp値を計算するか
n <- 30 #p値一個あたりの試行回数
a <- apply(matrix(runif(k*n), nr=n, nc=k), 1, sort)
b <- boxplot(t(a))

res0 <- NULL
for(i in seq(k)){
	res0 <- cbind(res0, t.test(a[i, ])$conf.int) #信頼区間を計算
}
boxplot(t(a))
for(i in seq(2)){lines(res0[i,], col=2, lwd=3)}