1. はじめに
多重検定について、2群間の検定であるt-検定等をそのまま繰り返すと落とし穴があります。
これから、具体的にRでシミュレーションしてみたいと思います。
まず、小学生６年生（全国平均）は、平均 149cm 標準偏差4cmのデータの確率密度関数(正規分布)に従うと仮定します（100万人の母集団とします）。
library(tidyverse)
rm(list=ls())
#小学6年生の母集団を設定する
set.seed(123)
syogakusei = rnorm(1000000, 149, 4)
2. 2群間の検定
まず最初に同じ母集団から、無作為に10人を非復元抽出しA、Bの2群として2群間の平均を比較することを試みます。
#A,Bの非復元抽出を指定
A = sample(syogakusei,10, replace=F)
B = sample(syogakusei,10, replace=F)
A,Bの2群の平均の違いについてt検定をします。
有意水準5%とすると、A,Bは同じ母集団からのサンプリングですが、この検定によってH0:帰無仮説が棄却される確率は5%存在します。
t.test(A,B)
>Welch Two Sample t-test
>
>data: A and B
>t = -0.26971, df = 17.943, p-value = 0.7905
>alternative hypothesis: true difference in means is not equal to 0
>95 percent confidence interval:
> -4.256628 3.288267
>sample estimates:
>mean of x mean of y
> 149.0848 149.5689
上記のp値が0.05以下であれば、H0：帰無仮説は棄却されるとなります。
この場合、同じ母集団からのサンプリングなので帰無仮説が正しいため、95％の確率で棄却されませんが、5%の確率で棄却される可能性があります。
それでは1000回繰り返し、検定が棄却される状況をシミュレーションします
set.seed(12345)
keka = rep(NA,10000) #記録するためのベクトル
for(i in 1:10000){
A = sample(syogakusei,10, replace=F)
B = sample(syogakusei,10, replace=F)
k=as.numeric(t.test(A,B)[3][1])
keka[i] = k #p値を記録する
}
さてp値が0.05以下で帰無仮説が棄却された回数はの割合は？
keka = keka <= 0.05
mean(keka)
>[1] 54
0.054とほぼ期待どおり５%の確率で棄却されたことが確認できました。
3. 3群間の検定 t-検定の繰り返し
次に同じことを3群間でt検定を繰り返してみます。
３群間（A,B,C）の場合、組み合わせは次のとおりとなります。
- A-B
- A-C
- B-C
set.seed(12345)
keka=rep(NA,1000)
i=1
for(i in 1:10000){
A = sample(syogakusei,10, replace=F)
B = sample(syogakusei,10, replace=F)
C = sample(syogakusei,10, replace=F)
k1 = as.numeric(t.test(A,B)[3][1])
k2 = as.numeric(t.test(A,C)[3][1])
k3 = as.numeric(t.test(B,C)[3][1])
kk =c (k1,k2,k3)
kt = sum(kk<=0.05) #3回の内、1回でも0.05以下があれば棄却される
if (kt>0) {
keka[i] = 1} else{keka[i] = 0
}
}
さて帰無仮説が棄却された回数の割合は？
mean(keka)
>[1]0.132
0.132と約13%となりました。
棄却される確率は5%なのはずなのに、なぜ上昇したのでしょうか？
3群間の検定の組み合わせをみると
2回連続で検定するということは同時確率になるので乗算します。
棄却しない確率95% x 95% =90% となるため、棄却される確率は10%になり、5%より上昇します。
このことから、本来5%の有意水準なので、有意とならない事象を有意と判定するタイプIのエラー率が上昇します。
これが多重検定時における罠になります。
4. 多重検定の有意水準の調整
多重検定では有意水準は、そのままでは有意水準が適用できないため、調整する必要がでてきます。
1つ目は、全グループのデータのばらつきを元に各群同士に差があるかを検定していく方法です。Tukey-Kramer法がこれに該当します。
２つ目は、先述の「本当は差がないものの中から誤って一つでも差がある判定としてしまう確率（Familywise error rate）」を5%に調整する方法です。もとの有意水準を検定の回数で割って新たな有意水準とするBonferroni法などがこれに該当します。
まずBonferroni法 を使ってみます。
５％有意水準で検定を2回繰り返すため、調整済みの有意水準は0.05/2 = 0.025となります。
先ほどと同じように1000回繰り返し、棄却状況についてシミュレーションします。
set.seed(12345)
keka=rep(NA,1000)
i=1
for(i in 1:1000){
A = sample(syogakusei,10, replace=F)
B = sample(syogakusei,10, replace=F)
C = sample(syogakusei,10, replace=F)
k1 = as.numeric(t.test(A,B)[3][1])
k2 = as.numeric(t.test(A,C)[3][1])
k3 = as.numeric(t.test(A,C)[3][1])
kk = c(k1,k2,k3)
kt = sum(kk<=(0.05/2)) #Bonferroni修正を実施
if (kt>0) {
keka[i] = 1} else{keka[i] = 0
}
}
mean(keka)
>[1] 0.055
Bonferroni修正前の棄却された割合が0.132だったのに対し、Bonferroni修正後は0.055になり、期待値である5%に近くなりました。
次に別の調整法であるTukey-Kramer法を使って調整を試みます。
まずデータフレーム化してtidyデータ化します。
data = data.frame(A=A,B=B,C=C)
data_tidy = data %>% pivot_longer(A:C, names_to = "treatment",values_to = "value")
多重検定には次の関数を使います。
TukeyHSD( aov( ) )
TukeyHSD(aov(value ~ treatment,data = data_tidy))
> Tukey multiple comparisons of means
> 95% family-wise confidence level
>
>Fit: aov(formula = value ~ treatment, data = data_tidy)
>
>$treatment
> diff lwr upr p adj
>B-A 1.3581244 -2.403252 5.119500 0.6478809
>C-A 0.9400231 -2.821353 4.701399 0.8106926
>C-B -0.4181013 -4.179477 3.343275 0.9590564
調整されたp値( p adj)が一番右に表示されています。
それでは今までと同様にTukey-Kramer法を1000回繰り返し、シミュレーションします。
set.seed(12345)
keka=rep(NA,1000)
i=1
for(i in 1:1000){
A = sample(syogakusei,10, replace=F)
B = sample(syogakusei,10, replace=F)
C = sample(syogakusei,10, replace=F)
data = data.frame(A=A,B=B,C=C)
data_tidy = data %>% pivot_longer(A:C, names_to = "treatment",values_to = "value")
k = TukeyHSD(aov(value ~ treatment,data = data_tidy))$treatment
kk = k[,4]
kt = sum(kk<=0.05)
if (kt>0) {
keka[i] = 1} else{keka[i] = 0
}
}
さて棄却された割合は？
mean(keka)
>[1] 0.051
修正前の棄却された割合0.132に対し、修正後は0.051になり、期待値である5%に近くなりました。
5.まとめ
今回、ご紹介した方法以外にも調整方法はありますが、今回の方法についてまとめます。
〇多重検定の使い分け方（パラメトリック）
1.比較する各群の分散が同じであると仮定する場合
コントロールがある場合 : Dunnett検定 (対照群と処置群との比較)
パッケージ multcompやSimComp等を使います。
コントロールがない場合 : Tukey-Kramer検定 (全ての群間比較)
上記のTukeyHSD( )関数や、パッケージ multcompを使います。
2.比較する各群の分散が異なると仮定する場合
Bonferroni補正を使います。
間違いがあれば、ご指摘いただけると大変ありがたいです。
6.参考