検定方法ではなく有意水準を調整する

　多重比較法としてTukey-Kramer法を紹介したが、これはt検定を多重比較用に改良することで多重比較に対応していた。今回紹介するBonferroniの補正やHolmの補正は、検定を改良するのではなく、有意水準を調整する、という考えに基づいている。最近は見かけることが減った気がするし、後述するように弱点もあるのであるが、幅広い状況に対応できるのが魅力の一つである。考え方もシンプルで合理的なので、検定の多重性の問題の理解にも役立ってくれるだろう。

Bonferroniの補正と証明

　Bonferroniの補正は有意水準αの検定をn回行ったときに、有意水準をα/nに補正して、帰無仮説の棄却の判定を行う。例えば、3標本A、B、Cがあるときにこれらから2標本取り出してt検定を行うとする。すると、AB、BC、CAの3ペアで検定を行うので、検定回数は3回である。もともとの有意水準を5%としていたなら、5/3 = 1.66%で帰無仮説の判定をする。AB、BC、CAの3ペアで検定で、それぞれp = 0.04、0.02、0.003だとしたら、もともとの有意水準ではすべて有意になるが、補正後はABとBCでは帰無仮説を棄却せず、CAでは棄却する、ということになる。

　Bonferroniの補正の主張は、有意水準α/nなら、n回検定を行っても危険率がαを越えない、というものだ。式で表すなら、1 - α/nが1回の検定で第1種の過誤を犯さない確率でn回の検定であれば、(1 - α/n)^nが第1種の過誤を犯さない確率である。これを1から引いた、1 - (1 - α/n)^nがn回検定した時の危険率なので、Bonferroniの主張は1 - (1 - α/n)^n < αということだ。証明しておこう。

上記のように、f(x) = (1 - α/x)^xとしたとき、f(x)は単調増加関数なので、f(x) ≧ f(1) = 1 - αとなる。両辺を-1倍して、1を足せば、1 - f(x) ≦ αとなって題意を得られる。グラフでも確認しておこう。Rを使って以下のように書ける。

------------------------------------------------------

library(plotn)

library(multcomp)

a <- 0.05

x <- seq(1,30,length = 200)

Ba <- 1 - (1 - a/x)^x #Bonferroniの補正をした多重検定の危険率

Na <- 1 - (1 - a)^x #補正なしの多重検定の危険率

as <- data.frame(x = c(x,x), a = c(Ba, Na), type = c(rep("Bonferroni", length(x)),rep("Non-adjust",length(x))))

plotn(a ~ x, data = as, group = "type", mode = "m", xlim = c(0, 20),

ylim = c(0.045, 0.055),

xlab = "検定回数", ylab = "危険率", line = T, pch = NA) #図1の描画

overdraw(abline(h = a, col = "black"))

------------------------------------------------------

無補正であると危険率が急激に上昇するが、Bonnferroniの補正によって多重検定をしても危険率は5%以下になっている。

Holmの補正

　Bonferroniの補正は検定数が少ないときは問題にならないが、検定数が増えると危険率は低いものの検出力も過剰に低くなる傾向にあることが知られている。特に全標本ペアについて差を知りたいときは注意が必要で、標本数をnとすれば2標本の組み合わせはn(n-1)/2の数になる。つまり、検定数は標本数に対して2乗のオーダーで急速に増えるということである。3標本なら3回の検定で済むが、5標本なら10回、7標本なら21回検定が必要だ。Bonferroniの補正は有意水準を一律で有意水準/検定回数に補正するので、7標本なら0.24%を下回る必要がある。これはかなり厳しい基準だ。そこで、Bonferroniの補正の改良としてHolmの補正がある。手順は以下の通りだ。

Holmの補正は、最も低いp値に関してはBonferroniの補正と同程度の有意水準で比較を行うため、危険率が補正されないときよりも低く、順次、有意水準を緩めながらp値と比較するため、Bonferroniの補正よりも検出力が高い。Bonferroniの補正がかなり有名であるが、Holmの補正も有名さのわりに、非常に有効な手段なので覚えておこう。

Bonferroni/Holmの補正の危険率

　ではBonferroniおよびHolmの補正について、その危険率を調べてみよう。併せて、Tukey-Kramer法とも比較をしてみる。以下にコードを示す。

------------------------------------------------------

a <- 0.05

n_samples <- 3

x <- n_samples*(n_samples - 1)/2 #検定回数

nA <- 30

nB <- 30

nC <- 30

meA <- 50

meB <- 50

meC <- 50

s <- 20

pAB <- NULL

pBC <- NULL

pCA <- NULL

ps <- NULL

for (i in 1:10000){

mA <- rnorm(nA, mean = meA, sd = s)

mB <- rnorm(nB, mean = meB, sd = s)

mC <- rnorm(nC, mean = meC, sd = s)

pAB <- c(pAB, t.test(mA, mB, var.equal = F)$p.value)

pBC <- c(pBC, t.test(mB, mC, var.equal = F)$p.value)

pCA <- c(pCA, t.test(mC, mA, var.equal = F)$p.value)

d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))

d$g <- as.factor(d$g)

fit1 <- glht(aov(d ~ g, data = d), linfct = mcp(g = "Tukey")) #Tukey-Kramer法

fit2 <- summary(fit1)

ps <- rbind(ps, fit2$test$pvalues[1:x]) #Tukey-Kramer法のp値を抽出

}

p <- data.frame(pAB, pBC, pCA)

p$pmin <- apply(p, MARGIN = 1, min) #独立に3回t検定した時のp値の最小値

#p.adjust(x, method = "xxx")で対応する補正でp値を計算しなおす

#p.adjust(x, method = "bonferroni")ならp値をp×検定回数にしている。

#その補正したp値と有意水準を比較し、有意水準を下回った数を計算している。

p$bonf <- apply(p[,-(x+1)], MARGIN = 1, FUN = function(x){sum(p.adjust(x, method = "bonferroni") < a)}) #Bonferroniの補正

p$holm <- apply(p[,-c(x+1,x+2)], MARGIN = 1, FUN = function(x){sum(p.adjust(x, method = "holm") < a)}) #Holmの補正

ps <- as.data.frame(ps) #Tukey-Kramer法のp値

colnames(ps) <- names(fit2$test$tstat)

ps$tukey <- apply(ps, MARGIN = 1, FUN = function(x){sum(x < a)}) #Tukey-Kramer法について有意水準を下回ったものの数を数える。

p <- cbind(p, ps)

head(p)

## pAB pBC pCA pmin bonf holm B - A C - A

## 1 0.7832167 0.5952187 0.7523131 0.5952187 0 0 0.9554652 0.9573710

## 2 0.6563377 0.3201086 0.5741571 0.3201086 0 0 0.9007714 0.8337667

## 3 0.6007019 0.9134639 0.4978993 0.4978993 0 0 0.8519327 0.7923190

## 4 0.7386638 0.4789175 0.3105236 0.3105236 0 0 0.9487198 0.5454544

## 5 0.2079796 0.1498576 0.8493621 0.1498576 0 0 0.4269381 0.9790919

## 6 0.1537798 0.1313544 0.8812203 0.1313544 0 0 0.3428923 0.9874012

## C - B tukey

## 1 0.8370766 0

## 2 0.5720885 0

## 3 0.9932088 0

## 4 0.7378393 0

## 5 0.3218508 0

## 6 0.2700742 0

#bonf、holm、tukeyの列に有意になった検定数が格納されている。

histn(p$pmin, xlab = "P value", ylab = "頻度") #図2の図示

#独立に3回検定した時の最小のp値の頻度

------------------------------------------------------

図2は検定の多重性問題の時と同じ図を描いた。今回は有意水準のほうを調整するのでp値の分布自体は変わらない。bonf、holm、tukey列を作ったので、ここが1以上になっている行数を計算すれば、1回以上の検定で帰無仮説を棄却していることがわかる。1回以上棄却した数を計算すれば、今回は平均値が等しい状態が正しいので第1種の過誤を犯した確率となる。

------------------------------------------------------

sum(p$pmin < 0.05) #無補正p値で1回以上帰無仮説を棄却した検定数

## [1] 1198

sum(p$bonf > 0) #Bonferroniの補正で1回以上帰無仮説を棄却した検定数

## [1] 441

sum(p$holm > 0) #Holmの補正で1回以上帰無仮説を棄却した検定数

## [1] 441

sum(p$tukey > 0) #Tukey-Kramer法で1回以上帰無仮説を棄却した検定数

## [1] 490

mean(p$bonf) #Bonferroniの補正で帰無仮説を棄却した検定数の平均数

## [1] 0.0498

mean(p$holm) #Holmの補正で帰無仮説を棄却した検定数の平均数

## [1] 0.0535

mean(p$tukey) #Tukey-Kramer法で帰無仮説を棄却した検定数の平均数

## [1] 0.0569

------------------------------------------------------

以上のようにp値の補正をしない場合は危険率が12%近くになっているが、Bonferroniの補正、Holmの補正、Tukey-Kramer法では5%程度の危険率となっている。Bonferroniの補正とHolmの補正は最も小さいp値に対しての有意水準は同じなので、危険率は同じとなる。また、棄却した検定数の平均値はいずれも0に近い値となっており、大差はない。続いて、棄却した検定数の分布を確認してみよう。これの確認は特に検出力の確認の時に重要になる。

------------------------------------------------------

histn(p$bonf, xlab = "棄却された検定数") #図3の図示

histn(p$holm, xlab = "棄却された検定数") #図4の図示

histn(p$tukey, xlab = "棄却された検定数") #図5の図示

------------------------------------------------------

図3～5でその分布は大きく変わらない。この結果から、3標本の比較においてはBonferroniの補正、Holmの補正およびTukey-Kramer法で使い勝手はそれほど変わらない。続いて、5標本の比較をした時の様子を確認してみよう。次のコードは非常に解析時間がかかるので、手元でやるときは100回程度の繰り返し程度にとどめておくことをお勧めする。ここでは10000回行う。

------------------------------------------------------

a <- 0.05

n_samples <- 5

x <- n_samples*(n_samples - 1)/2

n <- 30

me <- 50

s <- 20

pAB <- NULL

pAC <- NULL

pAD <- NULL

pAE <- NULL

pBC <- NULL

pBD <- NULL

pBE <- NULL

pCD <- NULL

pCE <- NULL

pDE <- NULL

ps <- NULL

for (i in 1:10000){ #手元で行うときは100回程度にしておこう

mA <- rnorm(n, mean = me, sd = s)

mB <- rnorm(n, mean = me, sd = s)

mC <- rnorm(n, mean = me, sd = s)

mD <- rnorm(n, mean = me, sd = s)

mE <- rnorm(n, mean = me, sd = s)

pAB <- c(pAB, t.test(mA, mB, var.equal = F)$p.value)

pAC <- c(pAC, t.test(mA, mC, var.equal = F)$p.value)

pAD <- c(pAD, t.test(mA, mD, var.equal = F)$p.value)

pAE <- c(pAE, t.test(mA, mE, var.equal = F)$p.value)

pBC <- c(pBC, t.test(mB, mC, var.equal = F)$p.value)

pBD <- c(pBD, t.test(mB, mD, var.equal = F)$p.value)

pBE <- c(pBE, t.test(mB, mE, var.equal = F)$p.value)

pCD <- c(pCD, t.test(mC, mD, var.equal = F)$p.value)

pCE <- c(pCE, t.test(mC, mE, var.equal = F)$p.value)

pDE <- c(pDE, t.test(mD, mE, var.equal = F)$p.value)

d <- data.frame(d = c(mA, mB, mC, mD, mE),

g = c(rep("A", n), rep("B", n), rep("C", n),

rep("D", n), rep("E", n)))

d$g <- as.factor(d$g)

fit1 <- glht(aov(d ~ g, data = d), linfct = mcp(g = "Tukey")) #Tukey-Kramer法

fit2 <- summary(fit1)

ps <- rbind(ps, fit2$test$pvalues[1:x])

}

p <- data.frame(pAB,pAC,pAD,pAE,pBC,

pBD,pBE,pCD,pCE,pDE)

p$pmin <- apply(p, MARGIN = 1, min)

p$bonf <- apply(p[,-(x+1)], MARGIN = 1, FUN = function(x){sum(p.adjust(x, method = "bonferroni") < a)})

p$holm <- apply(p[,-c(x+1,x+2)], MARGIN = 1, FUN = function(x){sum(p.adjust(x, method = "holm") < a)})

ps <- as.data.frame(ps)

colnames(ps) <- names(fit2$test$tstat)

ps$tukey <- apply(ps, MARGIN = 1, FUN = function(x){sum(x < a)})

p <- cbind(p, ps)

head(p)

## pAB pAC pAD pAE pBC pBD pBE

## 1 0.4494687 0.4619952 0.3965028 0.4802990 0.9802798 0.9691175 0.9721925

## 2 0.8821492 0.9658862 0.6138184 0.5059831 0.9316642 0.7352591 0.4654855

## 3 0.3862097 0.7865140 0.8287744 0.2710762 0.2706819 0.5070805 0.9272519

## 4 0.2317006 0.9818460 0.9567691 0.8330420 0.2339299 0.2753950 0.1743149

## 5 0.5325003 0.4588261 0.4183853 0.3301437 0.8632075 0.8624147 0.7459694

## 6 0.4016304 0.9743053 0.1712597 0.9122776 0.3995514 0.0226176 0.3340582

## pCD pCE pDE pmin bonf holm B - A C - A

## 1 0.9478896 0.9545697 0.9978236 0.3965028 0 0 0.9479638 0.9537527

## 2 0.6893521 0.5443744 0.3164450 0.3164450 0 0 0.9999381 0.9999994

## 3 0.6302274 0.1734927 0.3870740 0.1734927 0 0 0.8862120 0.9987497

## 4 0.9403339 0.8538000 0.7990113 0.1743149 0 0 0.7768633 0.9999999

## 5 0.9847419 0.9114985 0.8808277 0.3301437 0 0 0.9686875 0.9228611

## 6 0.1449025 0.8828668 0.2021596 0.0226176 0 0 0.9160841 0.9999997

## D - A E - A C - B D - B E - B D - C E - C

## 1 0.9386114 0.9378204 0.9999999 0.9999996 0.9999995 0.9999970 0.9999964

## 2 0.9889601 0.9731606 0.9999861 0.9966422 0.9468752 0.9920669 0.9661963

## 3 0.9995157 0.8468374 0.7510244 0.9536922 0.9999841 0.9882258 0.6966395

## 4 0.9999980 0.9995518 0.7643115 0.8061848 0.6502433 0.9999919 0.9997129

## 5 0.9290940 0.8816297 0.9997267 0.9998266 0.9982062 1.0000000 0.9999614

## 6 0.6142851 0.9999594 0.9266176 0.1616338 0.8723830 0.5931031 0.9998899

## E - D tukey

## 1 1.0000000 0

## 2 0.8102809 0

## 3 0.9295374 0

## 4 0.9988892 0

## 5 0.9999262 0

## 6 0.6873573 0

histn(p$pmin, xlab = "P value", ylab = "頻度") #図6の図示

------------------------------------------------------

続いて危険率を計算しよう。

------------------------------------------------------

sum(p$pmin < 0.05) #無補正p値で1回以上帰無仮説を棄却した検定数

## [1] 2794

sum(p$bonf > 0) #Bonferroniの補正で1回以上帰無仮説を棄却した検定数

## [1] 369

sum(p$holm > 0) #Holmの補正で1回以上帰無仮説を棄却した検定数

## [1] 369

sum(p$tukey > 0) #Tukey-Kramer法で1回以上帰無仮説を棄却した検定数

## [1] 468

mean(p$bonf) #Bonferroniの補正で帰無仮説を棄却した検定数の平均数

## [1] 0.0455

mean(p$holm) #Holmの補正で帰無仮説を棄却した検定数の平均数

## [1] 0.0471

mean(p$tukey) #Tukey-Kramer法で帰無仮説を棄却した検定数の平均数

## [1] 0.0617

------------------------------------------------------

p値の補正をしない場合は危険率が28%近くになっており、3~4回の複数回検定で1回の間違いを犯すことがわかる。一方で、Bonferroniの補正、Holmの補正、Tukey-Kramer法では5%程度の危険率となっている。危険率に関しては補正が有効に働いている。続いて、棄却した検定数の分布を確認してみよう。

------------------------------------------------------

histn(p$bonf, xlab = "棄却された検定数") #図7の図示

histn(p$holm, xlab = "棄却された検定数") #図8の図示

histn(p$tukey, xlab = "棄却された検定数") #図9の図示

------------------------------------------------------

以上のように5標本になっても、各補正とTukey-Kramer法の性能に大きな違いはない。3標本でも、5標本でも、Tukey-Kramer法よりも各補正のほうが危険率が低いくらいであり、これは図1の様子とも一致した結果となっている。

Bonferroni/Holmの補正の検出力

　ではBonferroniおよびHolmの補正について、その検出力を調べてみよう。3標本では、標本A、B、Cのうち、Aだけが20程度大きい平均を持つ母集団とする。このとき、3回の検定のうちAとB、AとCの間で有意になり、BとCの間では有意にならない、つまり帰無仮説の棄却をされる回数は2回である。

------------------------------------------------------

a <- 0.05

n_samples <- 3

x <- n_samples*(n_samples - 1)/2 #検定回数

nA <- 30

nB <- 30

nC <- 30

meA <- 70

meB <- 50

meC <- 50

s <- 20

pAB <- NULL

pBC <- NULL

pCA <- NULL

ps <- NULL

for (i in 1:10000){

mA <- rnorm(nA, mean = meA, sd = s)

mB <- rnorm(nB, mean = meB, sd = s)

mC <- rnorm(nC, mean = meC, sd = s)

pAB <- c(pAB, t.test(mA, mB, var.equal = F)$p.value)

pBC <- c(pBC, t.test(mB, mC, var.equal = F)$p.value)

pCA <- c(pCA, t.test(mC, mA, var.equal = F)$p.value)

d <- data.frame(d = c(mA, mB, mC), g = c(rep("A", nA), rep("B", nB), rep("C", nC)))

d$g <- as.factor(d$g)

fit1 <- glht(aov(d ~ g, data = d), linfct = mcp(g = "Tukey")) #Tukey-Kramer法

fit2 <- summary(fit1)

ps <- rbind(ps, fit2$test$pvalues[1:x]) #Tukey-Kramer法のp値を抽出

}

p <- data.frame(pAB, pBC, pCA)

p$pmin <- apply(p, MARGIN = 1, min)

p$bonf <- apply(p[,-(x+1)], MARGIN = 1, FUN = function(x){sum(p.adjust(x, method = "bonferroni") < a)}) #Bonferroniの補正

p$holm <- apply(p[,-c(x+1,x+2)], MARGIN = 1, FUN = function(x){sum(p.adjust(x, method = "holm") < a)}) #Holmの補正

ps <- as.data.frame(ps) #Tukey-Kramer法のp値

colnames(ps) <- names(fit2$test$tstat)

ps$tukey <- apply(ps, MARGIN = 1, FUN = function(x){sum(x < a)})

p <- cbind(p, ps)

head(p)

## pAB pBC pCA pmin bonf holm B - A

## 1 7.447779e-07 0.3068460 1.544818e-04 7.447779e-07 2 2 3.180969e-06

## 2 6.632140e-03 0.5957989 1.241446e-03 1.241446e-03 2 2 1.828742e-02

## 3 5.096932e-04 0.2964528 2.022780e-05 2.022780e-05 2 2 1.825286e-03

## 4 2.196696e-04 0.4693331 3.208982e-05 3.208982e-05 2 2 2.781116e-04

## 5 1.631154e-05 0.2634409 1.222334e-03 1.631154e-05 2 2 5.695799e-05

## 6 2.496160e-06 0.7249096 1.570326e-06 1.570326e-06 2 2 6.350351e-07

## C - A C - B tukey

## 1 1.666753e-04 0.5535358 2

## 2 3.604675e-03 0.8459836 2

## 3 3.501682e-05 0.5208166 2

## 4 2.241782e-05 0.7860305 2

## 5 4.445240e-03 0.4386492 2

## 6 2.198281e-06 0.9402612 2

histn(p$pmin, xlab = "P value", ylab = "頻度") #図10の図示

#独立に3回検定した時の最小のp値の頻度

------------------------------------------------------

次に検出力の計算を今まで通り、以下のように計算してみる。

------------------------------------------------------

sum(p$pmin < 0.05) #無補正p値で1回以上帰無仮説を棄却した検定数

## [1] 9933

sum(p$bonf > 0) #Bonferroniの補正で1回以上帰無仮説を棄却した検定数

## [1] 9736

sum(p$holm > 0) #Holmの補正で1回以上帰無仮説を棄却した検定数

## [1] 9736

sum(p$tukey > 0) #Tukey-Kramer法で1回以上帰無仮説を棄却した検定数

## [1] 9784

mean(p$bonf) #Bonferroniの補正で帰無仮説を棄却した検定数の平均数

## [1] 1.8474

mean(p$holm) #Holmの補正で帰無仮説を棄却した検定数の平均数

## [1] 1.9105

mean(p$tukey) #Tukey-Kramer法で帰無仮説を棄却した検定数の平均数

## [1] 1.876

------------------------------------------------------

ただし、多重比較における検出力は今までと比べて少し難しい。というのも、正しい状態は、標本AとB、AとCで帰無仮説が棄却され、BとCでは棄却されない、というものである。そうなら、少なくとも棄却された数が2以外であれば、間違いを犯していることになる（前回の項では分散分析との比較のために計算しなかった）。したがって上記の計算は多重比較において検出力を過大評価している。適切な検出力の確認のためには、棄却された検定数を分布を確認する必要がある。

------------------------------------------------------

histn(p$bonf, xlab = "棄却された検定数") #図11の図示

histn(p$holm, xlab = "棄却された検定数") #図12の図示

histn(p$tukey, xlab = "棄却された検定数") #図13の図示

------------------------------------------------------

図11～13を見るとわかるが、必ずしも棄却された検定数は2ではない。今回の場合は、いずれの方法でも大差はなさそうだ。以下のようにtable()関数を使って計算して、より正しい検出力を計算する（より正確にはAB、ACで有意になっているもののみを選ぶ必要があるが、煩雑になるのでここでは計算しない。興味がある人はやってみよう）。

------------------------------------------------------

table(p$bonf)

##

## 0 1 2 3

## 264 1088 8558 90

table(p$holm)

##

## 0 1 2 3

## 264 726 8651 359

table(p$tukey)

##

## 0 1 2 3

## 216 923 8746 115

------------------------------------------------------

今回の場合は、いずれの方法でも86～87％程度で、上記のように1回以上棄却した検定数と比べると、10%程度低い値だ。方法間で大差はなく、3標本の比較においては、いずれの方法であっても判断に大きな違いは生まれない可能性が高いだろう。

　では、次に5標本の場合を考える。今回はAとBがC、D、Eに対して平均値が20大きい状況を考える。すると、10回の2標本検定のうち、2（A or B）× 3（C or D or E）= 6回は有意差が出ることが期待される。

------------------------------------------------------

a <- 0.05

n_samples <- 5

x <- n_samples*(n_samples - 1)/2

n <- 30

me1 <- 50

me2 <- 70

s <- 20

pAB <- NULL

pAC <- NULL

pAD <- NULL

pAE <- NULL

pBC <- NULL

pBD <- NULL

pBE <- NULL

pCD <- NULL

pCE <- NULL

pDE <- NULL

ps <- NULL

for (i in 1:10000){ #手元で行うときは100回程度にしておこう

mA <- rnorm(n, mean = me2, sd = s)

mB <- rnorm(n, mean = me2, sd = s)

mC <- rnorm(n, mean = me1, sd = s)

mD <- rnorm(n, mean = me1, sd = s)

mE <- rnorm(n, mean = me1, sd = s)

pAB <- c(pAB, t.test(mA, mB, var.equal = F)$p.value)

pAC <- c(pAC, t.test(mA, mC, var.equal = F)$p.value)

pAD <- c(pAD, t.test(mA, mD, var.equal = F)$p.value)

pAE <- c(pAE, t.test(mA, mE, var.equal = F)$p.value)

pBC <- c(pBC, t.test(mB, mC, var.equal = F)$p.value)

pBD <- c(pBD, t.test(mB, mD, var.equal = F)$p.value)

pBE <- c(pBE, t.test(mB, mE, var.equal = F)$p.value)

pCD <- c(pCD, t.test(mC, mD, var.equal = F)$p.value)

pCE <- c(pCE, t.test(mC, mE, var.equal = F)$p.value)

pDE <- c(pDE, t.test(mD, mE, var.equal = F)$p.value)

d <- data.frame(d = c(mA, mB, mC, mD, mE),

g = c(rep("A", n), rep("B", n), rep("C", n),

rep("D", n), rep("E", n)))

d$g <- as.factor(d$g)

fit1 <- glht(aov(d ~ g, data = d), linfct = mcp(g = "Tukey")) #Tukey-Kramer法

fit2 <- summary(fit1)

ps <- rbind(ps, fit2$test$pvalues[1:x])

}

p <- data.frame(pAB,pAC,pAD,pAE,pBC,

pBD,pBE,pCD,pCE,pDE)

p$pmin <- apply(p, MARGIN = 1, min)

p$bonf <- apply(p[,-(x+1)], MARGIN = 1, FUN = function(x){sum(p.adjust(x, method = "bonferroni") < a)})

p$holm <- apply(p[,-c(x+1,x+2)], MARGIN = 1, FUN = function(x){sum(p.adjust(x, method = "holm") < a)})

ps <- as.data.frame(ps)

colnames(ps) <- names(fit2$test$tstat)

ps$tukey <- apply(ps, MARGIN = 1, FUN = function(x){sum(x < a)})

p <- cbind(p, ps)

head(p)

## pAB pAC pAD pAE pBC pBD

## 1 0.39455900 3.215918e-03 8.461087e-04 8.415939e-06 1.517574e-02 4.340646e-03

## 2 0.83797900 6.797058e-05 1.993153e-04 4.084134e-05 7.631945e-05 2.622943e-04

## 3 0.93276244 2.709786e-05 5.024432e-04 1.476192e-03 1.761931e-05 4.521463e-04

## 4 0.58260741 1.471369e-07 9.644995e-05 2.293575e-10 4.956545e-05 1.862616e-03

## 5 0.49425980 3.499058e-07 1.567346e-03 6.506127e-05 3.398514e-07 4.129370e-03

## 6 0.01981968 1.229728e-03 1.977115e-05 2.150835e-03 2.733201e-06 4.396183e-08

## pBE pCD pCE pDE pmin bonf holm

## 1 3.765940e-05 0.92319101 0.09464733 0.07804042 8.415939e-06 5 5

## 2 5.379815e-05 0.34595682 0.81982150 0.44117624 4.084134e-05 6 6

## 3 1.288731e-03 0.91791677 0.39013493 0.52439536 1.761931e-05 6 6

## 4 8.655813e-07 0.58327336 0.27416792 0.14756408 2.293575e-10 6 6

## 5 1.276428e-04 0.04504682 0.14870123 0.48656108 3.398514e-07 6 6

## 6 5.218374e-06 0.59194511 0.89828880 0.50066431 4.396183e-08 6 6

## B - A C - A D - A E - A C - B D - B

## 1 0.9201255 9.446741e-03 0.0069629236 6.540831e-06 9.617701e-02 7.614488e-02

## 2 0.9996365 1.368023e-04 0.0051217639 3.534024e-04 5.714462e-05 2.577572e-03

## 3 0.9999928 5.002292e-04 0.0008036849 1.036845e-02 7.778039e-04 1.058150e-03

## 4 0.9833666 1.091757e-05 0.0001488644 7.744512e-07 1.075739e-04 1.154698e-03

## 5 0.9625378 2.725950e-07 0.0030195625 1.823790e-04 7.968850e-06 2.496447e-02

## 6 0.1133601 9.821159e-03 0.0017727985 1.515406e-02 3.765275e-07 2.470213e-08

## E - B D - C E - C E - D tukey

## 1 2.584369e-04 0.9999799 0.3541997 0.4103442 4

## 2 1.625886e-04 0.8668620 0.9991559 0.9503834 6

## 3 1.307842e-02 0.9999643 0.9154472 0.9480117 6

## 4 1.675607e-06 0.9729076 0.8743446 0.5201660 6

## 5 2.306159e-03 0.2121871 0.6341068 0.9465185 6

## 6 8.613412e-07 0.9868075 0.9999145 0.9675499 6

histn(p$pmin, xlab = "P value", ylab = "頻度") #図14の図示

------------------------------------------------------

棄却された平均検定数だけ計算してみよう。

------------------------------------------------------

mean(p$bonf) #Bonferroniの補正で帰無仮説を棄却した検定数の平均数

## [1] 5.0117

mean(p$holm) #Holmの補正で帰無仮説を棄却した検定数の平均数

## [1] 5.2891

mean(p$tukey) #Tukey-Kramer法で帰無仮説を棄却した検定数の平均数

## [1] 5.2471

------------------------------------------------------

この様子をみると、Bonnferroniの補正の検出力が低い様子が見て取れる。それを分布でも確認してみよう。

------------------------------------------------------

histn(p$bonf, xlab = "棄却された検定数") #図15の図示

histn(p$holm, xlab = "棄却された検定数") #図16の図示

histn(p$tukey, xlab = "棄却された検定数") #図17の図示

------------------------------------------------------

今回正しい検出力を計算するためには、棄却された検定数が6回の部分を抽出する必要がある。図を見るとBonferroniの補正では、6回棄却された回数が5000回程度に対し、Holmの補正やTukey-Kramer法では6000回程度と、約10％も検出力に差がついている。では以下の計算で検出力を求めよう。

------------------------------------------------------

table(p$bonf)

##

## 0 1 2 3 4 5 6 7

## 41 97 342 737 1767 2065 4925 26

table(p$holm)

##

## 0 1 2 3 4 5 6 7 8

## 41 74 258 503 1298 1549 6105 151 21

table(p$tukey)

##

## 0 1 2 3 4 5 6 7 8

## 28 56 216 534 1382 1897 5842 43 2

------------------------------------------------------

上記の結果から、Bonferroniの補正では検出力は49%程度に対し、Holmの補正では61%、Tukey-Kramer法では58%であることがわかる。このように5標本で対比較を行うとなると、Bonferroniの補正では検出力が小さくなることが分かる。これは標本数が増えるほど顕著になる。その点で、Holmの補正はTukey-Kramer法並みの検出力を持ち、多標本の比較でも効果を発揮してくれる。

Bonferroni/Holmの補正の使いどころ

　BonferroniおよびHolmの補正の使いどころはどんな時だろう。もし、Tukey-Kramer法ができる前提条件が整っているのであれば、こちらの方がプログラムを実行するだけで計算できるし、間違いも少ない。特に標本数が多くて、比較数が多いときには手作業で行うとエラーが乗ってしまう可能性もある（一部を検定し忘れた、など。まあこれも自分でプログラムできるようになれば解消していく問題であるが）。一方で、どんな検定であれ、柔軟に対応できるのがBonferroniやHolmの補正だ。分散分析やTukey-Kramer法が使えないときに、ノンパラメトリックな多重比較、例えば、Kruskal-Wallis検定（分散分析のノンパラメトリック版）やSteel-Dwass検定（Tukey-Kramer法のノンパラメトリック版）を使えるならまだよいが、これらの検定の条件にも合わないときは、もっと仮定が弱いノンパラメトリック検定（例えば、Brunner-Munzel検定【外部リンクも参照】）をする必要も出てくる。しかし、こういうノンパラメトリック検定は、多重比較に拡張されておらず、2標本までしか取り扱えないことも多い。そんな時は、BonferroniやHolmの補正の出番となる。