スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

等分散検定からのt検定 vs Welch検定

独立2群間の統計学的仮説検定でよく教わる話として、「Studentのt検定は2群の母分散が等しい必要があるので、F検定などを行って等分散でなければMann-WhitneyのU検定を使う」というものがある。

一方、これが多重検定に該当するためにαエラー率が増加するとする批判もある。検索するとこの立場をとる代表的な人は生物化学研究所の井口氏や三重大学の奥村教授など。これらの方々は等分散性検定をすることは無意味であり最初からWelchのt検定を行うのが良いという立場を取っている。

調べてみると、こうした方々はモンテカルロシミュレーションでαエラー率の増加を実際に示している。しかし検定手法を評価するのにαエラーだけを示すのはフェアじゃないと思うので、検出力についてもシミュレーションしてみた。シミュレーション手法は群馬大学の青木先生の方法を参考にした。

具体的には、以下の条件の正規分布に従う乱数を発生させて2群を生成して有意水準5%で検定を100万回ずつ行い、有意と判定された割合を求めた。計算時間は64bit版R 3.1.0 on Macbook Air (Core i7, 1.7GHz) で約20時間。

・群1: 平均0、標準偏差1、サンプルサイズ10
・群2: 平均 [0.00, 0.25, 0.50, 0.75, 1.00]、標準偏差 [0.2, 0.6, 1.0, 1.4, 1.8]、サンプルサイズ [10, 20, 30]

群2で母平均が0のものの結果はαエラー率にあたり、それ以外のものは検出力にあたる(なぜならαエラー率=false positive rate、検出力=true positive rateだから)。

これら75通りについて、以下の2通りの検定を100万回ずつ行った。

1. F検定により有意水準5%での等分散性の検定
    ・有意差が出れば不等分散と判断しMann-WhitneyのU検定により二群検定
    ・有意差が出なければ等分散とみなしStudentのt検定
2. Welchのt検定のみ

次に結果を示す:
MannWhitney-Welch01 
縦軸が「実際に有意差ありと判定された回数/100万回」、横軸が群2の母平均である。
当たり前だが、二つの群の母平均が離れるほどpositive rateも増加していくことが分かる。
くどいようだが横軸が0のデータはαエラー率 (false positive rate) を、それ以外は検出力 (true positive rate) を表す。

検定に用いた有意水準 (5%) の位置に横線を引いてある。先述した青木先生のシミュレーションでは、横軸が0の場合にのみ注目してシミュレーションを行っていたわけである。

この結果を見ると、2群のサンプルサイズが等しいか母分散が等しい場合にはどちらの方法もαエラー率、検出力ともに大差はないようである。

さらに評価を進めるために、真のpositive rateとの比較を行った。群1の母平均と母分散をμ1, σ12、群2の母平均と母分散をμ2, σ22とすると、群1と群2の平均値の差はN(μ2 - μ1, σ22/n2 + σ12/n1) の正規分布に従うので、これをもとにpositive rateを計算することで真のpositive rateを算出できるはずである。こうして算出したpositive rateとシミュレーション結果との差をプロットすると次のようになる:
MannWhitney-Welch02 
今度は0の位置に横線を引いてある。この図では検出力の理論値からのずれを表しているので、この線に近いほど良い検定といえるだろう。

まだ傾向が少々分かりにくいので、横軸をeffect size(2群の平均の差を分散で割った量)にしてみる:
MannWhitney-Welch03 
この結果から、Welchのt検定のみを行う場合はeffect sizeが大きくなるに従って検出力が理論値よりも減少する傾向がある事が分かる。同時に、その挙動はサンプルサイズにはあまり影響されず安定していることも伺える。
一方等分散検定によりStudentのt検定とMann-WhitneyのU検定に振り分けた場合は、2群のサンプルサイズが等しいときはWelchのt検定とよく似た挙動を示すが、サンプルサイズに差がある場合には挙動が非常に荒ぶる事が分かる。

結論としては、等分散性検定を事前に行う従来の検定手順はサンプルサイズが等しいときや2群の母分散が等しいときには悪くないが、それ以外の時にはαエラー率、検出力が変動しやすいことが分かった。Welchのt検定もこれらの要素の影響を受けないわけではないが、従来の手順に比べると変動は少なく信頼できるように思う。

追記 (2014.07.16)

3つ目の解析(横軸をeffect sizeにしたもの)は適切な解析ではなかったかもしれない。というのも、Mann-WhitneyのU検定がそもそも「二群がほぼ等分散である」ことを仮定しているため、不等分散な検定に用いること自体が不適切だから。とはいえこの手法が一般的に使われているようなのでそれを批判するつもりでシミュレーションしたら思いのほか悪くなかった(後述)ので驚いたというのが正直なところ。

今回のシミュレーションで意外だったのは、「等分散検定による事前検定は不適切」という批判が必ずしも妥当とは言えないことが判明したことだった。二群のサンプルサイズが等しい場合のαエラー率と検出力はどちらもWelch検定単体の場合とほとんど変わらなかった。母分散に差がない場合も同様の結果だった。むしろサンプルサイズの不均等が不等分散時の弊害を顕在化させると考えることができる。

母分散の差(比)は厳密には知りえないことを考えると、二群検定をする場合には検定方法を吟味することも大切だが、それ以前に実験計画段階でサンプルサイズを等しくする努力の重要性を示す結果となったと言える。



今回のシミュレーションのRコード:

tests.power.core <- function(dist, mean.dif, sd.ratio, n1, n2, sim.size) {
x <- matrix(rnorm(n1*sim.size, mean = 0, sd = 1), nrow = n1)
y <- matrix(rnorm(n2*sim.size, mean = mean.dif, sd = sd.ratio), nrow = n2)
welch <- sapply(seq_len(sim.size), function(i) {t.test(x[,i], y[,i])$p.value})
v2m <- sapply(seq_len(sim.size), function(i) {if(var.test(x[,i], y[,i])$p.value > 0.05) t.test(x[,i], y[,i], var.equal = TRUE)$p.value else wilcox.test(x[,i], y[,i])$p.value})
t <- sapply(seq_len(sim.size), function(i) {t.test(x[,i], y[,i], var.equal = TRUE)$p.value})
return(c(sum(welch < 0.05),
sum(v2m < 0.05)))
}

tests.power.proc <- function(dist, mean.difs, sd.ratios, n1, n2s, sim.size) {
powers <- numeric()
for(n2 in n2s) {
for(sd.ratio in sd.ratios) {
for(mean.dif in mean.difs) {
powers <- c(powers, tests.power.core(dist, mean.dif, sd.ratio, n1, n2, sim.size))
}
}
}
return(powers)
}

tests.power <- function(dist, mean.difs, sd.ratios, n1, n2s) {
mean.size <- length(mean.difs); sd.size <- length(sd.ratios); n.size <- length(n2s)
powers <- data.frame(n1 = rep.int(n1, mean.size * sd.size * n.size * 2),
n2 = sort(rep.int(n2s, mean.size * sd.size * 2)),
sd.ratio = rep.int(sort(rep.int(sd.ratios, mean.size * 2)), n.size),
mean.dif = rep.int(sort(rep.int(mean.difs, 2)), sd.size * n.size),
type = rep.int(c("welch.only", "mann.after.f"), mean.size * sd.size * n.size),
power = rep.int(0, mean.size * sd.size * n.size * 2),
trial = rep.int(0, mean.size * sd.size * n.size * 2))
proc <- function(sim.size) {
powers$power <<- powers$power + tests.power.proc(dist, mean.difs, sd.ratios, n1, n2s, sim.size)
powers$trial <<- powers$trial + sim.size
}
return(list(power = function() { powers },
proc = proc))
}

labeller <- function(var, value){
if(var == "mean.dif") {
value <- paste("Mean diff =", value)
} else if(var == "sd.ratio") {
value <- paste("SD ratio =", value)
} else if(var == "n2") {
value <- paste("n2 =", value)
}
return(value)
}

test <- tests.power(rnorm,
mean.difs = seq(0, 1, len = 5),
sd.ratios = seq(0.2, 2.0, by=0.4),
n1 = 10,
n2 = seq(10, 30, by = 10))
test$proc(1e+6)

pp <- test$power()
sds<-sqrt(1/pp$n1 + pp$sd.ratio^2/pp$n2)
powers.ideal <- pnorm(qnorm(0.975, mean = 0, sd = sds), mean = pp$mean.dif, sd = sds, lower.tail = FALSE) + pnorm(qnorm(0.025, mean = 0, sd = sds), mean = pp$mean.dif, sd = sds)
pp$ideal <- powers.ideal
pp$sd.stat <- sds

ggplot(pp, aes(x = mean.dif, y = power/trial, colour = type, shape = type)) +
geom_hline(aes(yintercept = 0.05)) + geom_point() + geom_line() +
facet_grid(n2 ~ sd.ratio, labeller = labeller) +
xlab("Difference of population means of the two samples") +
ylab(expression(paste("Simulated power (",10^6," trials)", sep=""))) +
scale_colour_discrete(name = "Test Method",
breaks = c("mann.after.f", "welch.only"),
labels = c("Mann-Whitney U/Student's t", "Welch's t test")) +
scale_shape_discrete(name = "Test Method",
breaks = c("mann.after.f", "welch.only"),
labels = c("Mann-Whitney U/Student's t", "Welch's t test"))

ggplot(pp, aes(x = mean.dif, y = power/trial - ideal, colour = type, shape = type)) +
geom_hline(aes(yintercept = 0)) + geom_point() + geom_line() +
facet_grid(n2 ~ sd.ratio, labeller = labeller) +
xlab("Difference of means") + ylab("Difference from the theoretical true positive rate") +
scale_colour_discrete(name = "Test Method",
breaks = c("mann.after.f", "welch.only"),
labels = c("Mann-Whitney U/Student's t", "Welch's t test")) +
scale_shape_discrete(name = "Test Method",
breaks = c("mann.after.f", "welch.only"),
labels = c("Mann-Whitney U/Student's t", "Welch's t test"))

ggplot(pp, aes(x = mean.dif/sd.stat, y = power/trial - ideal, colour = type, shape = type)) +
geom_hline(aes(yintercept = 0)) + geom_point() + geom_line() +
facet_grid(n2 ~ ., labeller = labeller) +
xlab(expression(paste("Effect size (", Delta, mu/sigma, ")", sep = ""))) +
ylab("Difference between the simulated power and the theoretical true positive rate") +
scale_colour_discrete(name = "Test Method",
breaks = c("mann.after.f", "welch.only"),
labels = c("Mann-Whitney U/Student's t", "Welch's t test")) +
scale_shape_discrete(name = "Test Method",
breaks = c("mann.after.f", "welch.only"),
labels = c("Mann-Whitney U/Student's t", "Welch's t test"))

スポンサーサイト

コメントの投稿

管理者にだけ表示を許可する

プロフィール

前世

Author:前世
Tipsとかてきとうに
趣味スクリプト書きだよ。
リグルが大好きだよ。
(元の生息地はここ)

Rユーザー
・Python/Javascript
・Applescriptもたまに書く

りぐる数
最新記事
月別アーカイブ
最新コメント
最新トラックバック
カテゴリ
検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

QRコード
QR
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。