一つのデータ
X
から複数の検定を行うとする.各検定の帰無仮説をH
1, H
2, . . . , H
kとし,集合
C
i(α)
を帰無仮説H
i の下での有意水準α
の棄却域とする.全ての帰無仮説H
i が同時に成り立つという帰無仮説H
を考える.ある帰無仮説H
i が棄却されればH
も棄却される(
おなじことだが,全てのH
i が棄却されなければH
も棄却されない)
よう にするのが一見自然と思われる(
これは∪
iC
i(α)
をH
の棄却域にすることと同値)
.し かし,こうした検定の有意水準は普通α
にはならず,一般的にいえることはP { X ∈ ∪
iC
i(α) } ≤ X
i
P { X ∈ C
i(α) } = kα
だけである(
検定の多重比較問題)
.ここでP { X ∈ ∪
iC
i(α/k) } ≤ k(α/k) = α
であるから
∪
iC
i(α/k)
は帰無仮説H
の有意水準α
の棄却域を与える.つまり,各検定H
i を有意水準α/k
で行えば,その結果は同時帰無仮説H
の有意水準α
の検定として使 える.これはBonferroni
の検定の多重比較補正と呼ばれる方法である.これは検定の多 重比較補正の基本となる考え方であるが,当然k
が少し大きくなるとα/k
が小さくなり すぎ,実際上無意味となる.検定の多重比較補正法としては,より精緻な方法が幾つも提案されているが,いずれも 個別の検定をある有意水準
α
i で行えば,それが結果として同時帰無仮説H
の有意水準α
の検定となるような補正有意水準α
1, α
2, . . . , α
k を与える(
一般にα
i≤ α )
.同じこ とをp
値の言葉で表現すると,本来のp
値p
1, p
2, . . . , p
k の多重比較補正値と呼ばれる 数値p
∗1, p
∗2, . . . , p
∗k(
一般にp
∗i≥ p
i)
を定め,あるi
でp
∗i≤ α
であれば,同時帰無仮説H
が有意水準α
で棄却されたことになるようにする.例えばBonferroni
の多重比較補 正法に対しては,p
値の多重比較補正をp
∗i= kp
iと取れば良い:
p
∗i≤ α ⇐⇒ p
i≤ α/k ⇐⇒ H
i を有意水準α/k
で棄却.3.3.1 p 値のベクトルの多重比較補正 p.adjust()
p.adjust()
はp
値のベクトルを与えて,幾つかの方法を用いて多重比較補正されたp
値のベクトルを返す.
書式:
p.adjust(p, method=p.adjust.methods, n=length(p)) p.adjust.methods #
可能な手法名の文字ベクトル 引数:p p
値のベクトルmethod
補正方法.以下の文字列のいずれか"holm","hochberg","hommel","bonferroni","BH","BY","fdr","none"
n
比較の個数返り値: 補正された
p
値のベクトル補正方法の一つは
Bonferroni
補正("bonferroni")
で,p
値は比較回数倍される.それほど保守的でない補正法がそれぞれ
Holm ("holm")
,Hochberg ("hochberg")
,Hommel ("hommel")
そしてBenjamini & Hochberg ("fdr")
で与えられる.無補正オ プション("none")
もある.補正法をオプションとして持ち,それを引き数p.adjust
に指定する必要がある手法の便宜用に,手法名の文字列ベクトルp.adjust.methods
が ある.最初の
4
つの方法はファミリ単位の誤差率を強力にコントロールするように作られている.
Bonferroni
補正はそのままでは,如何なる場合も使用可能なHolm
法より劣るので使用する理由は無い.
Hochberg
とHommel
の方法は,検定が独立か,非負の連関を持 つ(Sarkar, Sarkar & Chang)
という仮定下で有効である.Hommel
法はHochberg
法よ りもより強力であるが,差異は普通小さく,Hochberg
法はp
値をより高速に計算する.Benjamini & Hochberg
の"fdr"
法は棄却された仮説中の誤った発見の平均の割合であ る「誤発見率(false discovery rate)
」をコントロールする.誤発見率はファミリ単位の誤 判別率よりも厳しい条件であり,従ってBenjamini & Hochberg
法は他の手法よりもよ り強力な方法である.関連:
stats
パッケージ中のpairwise.t.test()
等のpairwise.*()
関数.> p.adjust.methods # 可能な手法名の文字ベクトル
[1] "holm" "hochberg" "hommel" "bonferroni" "fdr" "none"
> x <- rnorm(10, m=c(rep(0,5),rep(3,5))) # 平均0,3の正規乱数5個ずつ
> p <- 2*pnorm(-abs(x)) # テスト用のp値ベクトル
> round(p, 3)
[1] 0.165 0.282 0.621 0.435 0.696 0.040 0.001 0.102 0.001 0.010
> round(p.adjust(p), 3) # 既定のHolm法による補正
[1] 0.825 1.000 1.000 1.000 1.000 0.280 0.006 0.615 0.006 0.084
> round(p.adjust(p,"bonferroni"), 3) # Bonferoni法による補正 [1] 1.000 1.000 1.000 1.000 1.000 0.399 0.006 1.000 0.007 0.105
> round(p.adjust(p,"fdr"), 3) # false discovery rate法による補正 [1] 0.275 0.402 0.690 0.544 0.696 0.100 0.003 0.205 0.003 0.035
3.3.2 多重比較補正を伴う対毎の比率の比較 pairwise.prop.test()
pairwise.prop.test()
は,多重比較補正を伴う対毎の比率の比較を行う.書式:
pairwise.prop.test(x, n, p.adjust.method=p.adjust.methods, ...)
引数:
x
成功回数のベクトル.または,それぞれ成功・失敗回数を2
つの列に持つ行列n
試行回数のベクトル.もしx
が行列なら無視されるp.adjust.method p
値を補正する方法(p.adjust()
を参照) ... prop.test()
に引き渡される追加引き数返り値: クラス
"pairwise.htest"
を持つオブジェクト 関連:prop.test(), p.adjust()
.# 患者4グループ中の喫煙者数のデータ.各グループ対毎に喫煙者の割合が同一という帰無仮説の検定
# prpo.test()の参考例を参照(4グループ込みの比率の同一性の検定)
> smokers <- c(83, 90, 129, 70)
> patients <- c(86, 93, 136, 82)
# 各グループ対毎の喫煙者数比率の同一性を検定し,p値をHolm法で多重比較補正
> pairwise.prop.test(smokers, patients)
$method
[1] "Pairwise comparison of proportions"
$data.name
[1] "smokers out of patients"
# 対毎の比較の補正済みp値.第2,4グループ間の補正p値が0.1未満なので帰無仮説は10%有意
$p.value
1 2 3
2 1.0000000 NA NA
3 1.0000000 1.00000000 NA 4 0.1185648 0.09321728 0.1237680
$p.adjust.method # Holm法で多重比較補正(既定)
[1] "holm"
attr(,"class") [1] "pairwise.htest"
(以下省略)
> pairwise.prop.test(smokers, patients, p.adjust.method="fdr") # "fdr"法で多重比較補正 Pairwise comparisons using Pairwise comparison of proportions
data: smokers out of patients
1 2 3
2 1.000 - 3 0.965 0.965 -4 0.062 0.062 0.062
P value adjustment method: fdr (以下省略)
3.3.3 多重比較補正を伴うグループ水準間の対毎の t 検定による比較
pairwise.t.test()
pairwise.t.test()
は,多重比較補正を伴うグループ水準間の,対毎のt
検定による比較を行う.
書式:
pairwise.t.test(x, g, p.adjust.method = p.adjust.methods, pool.sd = TRUE, ...)
引数:
x
応答ベクトルg
グルーピング用のベクトルか因子p.adjust.method p
値を補正する方法(p.adjust()
を参照)
pool.sd
プールした標準偏差の使用を許すか許さないか... t.test()
に引き渡される追加引き数返り値: クラス
"pairwise.htest"
を持つオブジェクト関連:
t.test(), p.adjust()
.# NY市の毎日の大気観測値データairquality.各月対毎にオゾン量平均が等しいという帰無仮説のt検定
> attach(airquality) # 成分を名前で参照できるようにする
> Month <- factor(Month, labels = month.abb[5:9]) # 月数を省略文字名で因子化
> pairwise.t.test(Ozone, Month)
$method # プールした標準偏差を用いたt検定
[1] "t tests with pooled SD"
$data.name
[1] "Ozone and Month"
# 多重比較補正済みp値.(Jul,Sep),(Aug,Sep)間の補正p値が0.005より小さく,0.5%有意
$p.value
May Jun Jul Aug
Jun 1.0000000000 NA NA NA
Jul 0.0002638036 0.05112741 NA NA
Aug 0.0001949061 0.04987333 1.000000000 NA Sep 1.0000000000 1.00000000 0.004878798 0.003878108
$p.adjust.method # Holm法で多重比較補正法(既定)
[1] "holm"
attr(,"class") [1] "pairwise.htest"
# Bonferroni法で多重比較補正すると(Aug,Sep)間のみでp値が0.005以下
> pairwise.t.test(Ozone, Month, p.adj = "bonf")
$method # プールした標準偏差を用いたt検定
[1] "t tests with pooled SD"
$data.name
[1] "Ozone and Month"
$p.value
May Jun Jul Aug
Jun 1.0000000000 NA NA NA
Jul 0.0002931151 0.10225483 NA NA
Aug 0.0001949061 0.08312222 1.000000000 NA Sep 1.0000000000 1.00000000 0.006969712 0.004847635
$p.adjust.method # Bonferroni法で多重比較補正
[1] "bonferroni"
attr(,"class") [1] "pairwise.htest"
> pairwise.t.test(Ozone, Month, pool.sd = FALSE)
$method # プールした標準偏差を用いないt検定
[1] "t tests with non-pooled SD"
$data.name
[1] "Ozone and Month"
$p.value # 今度は1%有意という結果
May Jun Jul Aug
Jun 1.0000000000 NA NA NA
Jul 0.0002646679 0.01527179 NA NA
Aug 0.0019517090 0.02134659 1.000000000 NA Sep 0.8632062149 1.00000000 0.005888826 0.01720974
$p.adjust.method # Holm法で多重比較補正
[1] "holm"
attr(,"class") [1] "pairwise.htest"
> detach() # 成分を消去
3.3.4 Wilcoxon のランク和検定に対する多重比較補正
pairwise.wilcox.test()
pairwise.wilcox.test()
はグループ水準の各組合せ毎にWilcoxon
のランク和検定 を実施し,p
値に対する多重比較補正を行う.書式:
pairwise.wilcox.test(x,g,p.adjust.method=p.adjust.methods,...)
引数:x
応答ベクトルg
グルーピング用ベクトル,もしくは因子p.adjust.method p
値を補正する方法(p.adjust()
を参照) ... wilcox.test()
に引き渡される追加引き数返り値: クラス属性
"pairwise.htest"
を持つオブジェクト 関連:wilcox.test()
,p.adjust()
.# NY市の大気観測値データairquality使用.月の組合せ毎にwilcox.test()を適用し多重比較補正
> Month <- factor(Month, labels = month.abb[5:9]) # 5,6,7,8,9月名を因子に変換
> pairwise.wilcox.test(Ozone, Month)
$method
[1] "Wilcoxon rank sum test"
$data.name
[1] "Ozone and Month"
# 多重比較補正されたp値.(May,Jul)間で0.0003以下なので0.05%有意
$p.value
May Jun Jul Aug
Jun 0.5775049432 NA NA NA
Jul 0.0002996390 0.08481748 NA NA
Aug 0.0010872705 0.12953880 1.000000000 NA Sep 0.4743723336 1.00000000 0.005954083 0.02273569
$p.adjust.method # Holm法で補正(既定)
[1] "holm"
attr(,"class") [1] "pairwise.htest"
Warning messages: # データにタイがあるので警告がでる(省略)
# Bonferroni法で多重比較.再び0.05%有意という結果
> pairwise.wilcox.test(Ozone, Month, p.adj = "bonf")
$method
[1] "Wilcoxon rank sum test"
$data.name
[1] "Ozone and Month"
$p.value
May Jun Jul Aug
Jun 1.0000000000 NA NA NA
Jul 0.0002996390 0.1413625 NA NA
Aug 0.0012080783 0.2590776 1.000000000 NA Sep 1.0000000000 1.0000000 0.007442604 0.03247955
$p.adjust.method [1] "bonferroni"
attr(,"class") [1] "pairwise.htest"
Warning messages: # データにタイがあるので警告がでる(省略)
> detach()