<< トップページへ

<< 目次へ

最終更新日: 2008.11.24

分散分析法における多重比較



Contents
  1. p値と有意水準・危険率
  2. 第1種の誤りと第2種の誤り
  3. p値の多重性
  4. 多重性を考慮した検定法
  5. 交互作用項が有意である場合の多重比較
  6. 一般線形モデルにおける多重比較



■ p値と有意水準・危険率

 統計的仮説検定においてはp値によって帰無仮説(H0)が棄却されるかどうかを判断することになっている。そして帰無仮説を棄却するかどうかという判断の基準として有意水準(α)をあらかじめ設定しておき、得られたp値が有意水準よりも小さければ帰無仮説を棄却するのである。有意水準は一般的に0.05と設定されることが多いが、医療統計学の分野ではしばしば0.01や0.001という、より厳しい有意水準が設定されることもある。

 ここで改めてp値とは何かということを考えてみよう。p値とは「帰無仮説が正しいという条件の下で偶然にその検定統計量が得られる確率(条件つき確率)のこと」なのである。つまり、得られたp値が非常に小さい場合(例えばp = 0.0001である場合)は偶然に検定統計量が得られたという確率は非常に小さいということである。もう少しくだけた言い方をすれば「たまたま得られた結果であるとは考えられない」ということである。

 2群の平均値の差の検定(独立2標本のt検定)の例を通して、以上のことを考えてみよう。

# 第1群のデータ
> G1 <- c(78, 70, 66, 76, 78, 76, 88, 76)

# 第2群のデータ
> G2 <- c(55, 60, 63, 50, 69, 55, 62, 65)

# t検定を行う
> t.test(G1, G2)

        Welch Two Sample t-test

data:  G1 and G2 
t = 5.1135, df = 13.984, p-value = 0.0001583
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
  9.360791 22.889209 
sample estimates:
mean of x mean of y 
   76.000    59.875 

 ここで得られた検定統計量はt = 5.1135である。この値が偶然によって得られた値である確率がp値であるわけで、今回の検定で得られたp値はp = 0.0001583と非常に小さい。これはt = 5.1135という値が偶然によって得られる確率が0.0001583である、ということを意味している。

 それではp値がどれほどの大きさであれば「偶然でない」といえるのだろうか。その判断の基準となるのが有意水準である。前述したように、一般的に有意水準はα = 0.05と設定されることが多い。これは得られたp値が0.05よりも小さければ、検定統計量が偶然によって得られたものではないだろう、とするものである。

 ところでこの有意水準は危険率とも呼ばれている。それは、統計的仮説検定は常に誤った結論を導く確率も有しているからである。例えば、先ほどの検定で得られたp = 0.0001583という値は、0.0001583という確率で誤った結論を導く危険を有しているということができるのである。言い換えれば約10000回に1回は誤った結論を得るということである。

 したがって、有意水準をα = 0.05と設定するということは、同じ検定を100回行ったら、そのうちに5回は誤った結論を得る確率であるといえる。つまり、100回に5回の失敗なら許容しようというものである。先に述べたように、医療統計学の分野でα = 0.001とするのは1000回に1回の失敗なら許容するが、100回に5回の失敗では失敗の確率が高すぎるので許容できないという理由からそのようにするのである。もし風邪薬を投与した場合に100回に5回は副作用が現われるというよりは、1000回に1回の副作用が現われることのほうがより安心であろう、ということである。

<< このページのトップへ



■ 第1種の誤りと第2種の誤り

 α = 0.05という値は100回中5回は誤った結果を得てしまう確率であると述べた。では「誤った結果」とはなんのことであろうか。それは「帰無仮説が正しいのにもかかわらず、帰無仮説を棄却してしまう」という誤りのことである。この誤りのことを特に第1種の誤りといっている(第1種の過誤とかType I errorとも表現される)。ややこしいが、結局はαは有意水準、危険率、第1種の誤り、いずれの意味も含んでいるのである。

 これに対して「帰無仮説が誤っているにもかかわらず、帰無仮説を採択してしまう」という誤りのことを第2種の誤りといい、これをβで表すのが慣わしである。

 ここで統計的仮説検定において、起こりうる4つの結果を以下にまとめておく。

 第1種の誤りがαならば、「帰無仮説が正しく、帰無仮説を採択する」確率は1 - αで表せる。また、一方で第2種の誤りがβならば、「帰無仮説が誤りで、帰無仮説を採択する」確率は1 - βで表すことができる。

 統計的仮説検定において、多くは帰無仮説を棄却することが分析者にとって望ましいことであることが多いだろう。したがって、αを最小に抑えつつ(1 - β)が大きくなるように標本設計を行うことが好ましい。一般にサンプルサイズが大きくなるほど(1 - β)は大きくなるが、サンプルサイズを必要以上に大きくしすぎると、今度はかえってαを大きくしてしまうので注意が必要である。

 このようにαと(1 - β)の関係を考慮して適切なサンプルサイズを決めることをパワーアナリシスという。これは日本語で検定力検定とか検出力検定と呼ばれることがあるが、それは(1 - β)のことを検定力、あるいは検出力と呼んでいるからである。

 この内容については永田「サンプルサイズの決め方」朝倉書店に詳しい解説があるので、興味のある方はぜひとも参考にしてもらいたい。また、Rにはパワーアナリシスを行うための関数として、例えば、power.t.test()などといった関数がいくつか用意されている。英語が苦にならなければG*Powerというパワーアナリシスを行うためのフリーウェアが用意されている。

<< このページのトップへ



■ p値の多重性

 ある同一の検定を繰り返し行ったとすると、そこには多重性の問題が生じることになる。例えば、A、B、Cという3つのグループの平均値を2グループずつ対比較してみたいとしよう。その場合、考えらえる組み合わせは{A vs B, A vs C, B vs C}という3通りの組み合わせが考えられる。つまり、3回の検定を繰り返し行うことになる。

 実際にその様子を以下に示してみよう(Rのt.test()が返す出力の1部は省略してある)。

> A <- c(12, 18, 18)
> B <- c(20, 19, 26)
> C <- c(26, 32, 31)
> t.test(A, B)

        Welch Two Sample t-test

data:  A and B 
t = -1.9126, df = 3.969, p-value = 0.1289
mean of x mean of y 
 16.00000  21.66667 

> t.test(A, C)

        Welch Two Sample t-test

data:  A and C 
t = -5.0089, df = 3.978, p-value = 0.007552
mean of x mean of y 
 16.00000  29.66667 

> t.test(B, C)

        Welch Two Sample t-test

data:  B and C 
t = -2.7899, df = 3.898, p-value = 0.05078
mean of x mean of y 
 21.66667  29.66667

 検定の結果、α = 0.05という有意水準のもとではA vs Cの対比較のみ有意であることが認められた。しかし、このやり方は正しくない。

 しつこいようだが、α = 0.05という値は100回に5回は誤った結果(第1種の誤り)を得る確率のことである(言い換えると100回に5回の失敗なら許容するということ)。つまり、5%の確率で"ミスる"可能性があるというわけだ。ということは、同じ検定を3回繰り返したら、5%の確率で"ミスる"という確率は5%よりも大きくなると考えられるだろう。

 もう少し分かりやすい例をあげてみよう。例えば、赤信号の横断歩道を渡ると5%の確率で事故に会うとしよう。度胸のある人なら「100回に5回しか事故に会わないのだから」といって渡るかもしれない。実際、5%くらいの確率ならそれほど危険に感じないかもしれない。

 ところが、事故に会う確率が5%であったとしても、2回、3回、4回と渡る回数を重ねていくほど事故に会う確率は高くなっていくことが直感的に分かるだろう。これが統計的仮説検定における多重性の問題というものである。

<< このページのトップへ



■ 多重性を考慮した検定法

 統計学の初学者にとって悩まされる問題の1つは「どの多重比較法を用いたらよいか」ということであろう。たしかに、多重性を考慮した検定法は数多く存在するが、ここではボンフェローニ法とそれに関連する多重比較法のみを紹介する。なぜなら、ボンフェローニ法はどのような場合においても適用することができ、かつその理論が分かりやすく検定の際に特別な計算をする(プログラムを用いる)必要がないからである。とはいっても、この方法にもいくつかの問題点もあるので、それを含めて説明していくことにしよう。


ボンフェローニ法

 ボンフェローニ法とはボンフェローニの不等式に基づいて行われる多重比較法の1つである。ボンフェローニ法はk個の帰無仮説を検定する際に、全体の有意水準をαとしたとき、個々の検定の有意水準をα / kとする方法である。したがって、この方法を採用した場合、有意水準をコントロールするだけなので、通常通りのt検定によって繰り返し検定を行うだけでよい。

 先ほどのデータセットを用いて、ボンフェローニ法による多重比較法を適用してみよう。

A 12, 18, 18
B 20, 19, 26
C 26, 32, 31

 既に述べたように、考えられる全ての組み合わせは{A vs B, A vs C, B vs C}の3通りである。すなわち、ここには3つの帰無仮説が存在しているので、検定を3回繰り返すことになる。

HA vs B: A群とB群の平均値に差はない

HA vs C: A群とC群の平均値に差はない

HB vs C: B群とC群の平均値に差はない

 全体の有意水準をα = 0.05とすると、k = 3なので、個々の検定における有意水準はα' = α / k = 0.05 / 3 = 0.01666667となる。

 検定統計量、およびそれに対するp値は通常通りのt検定を行えばよい。以下は先ほどの解析と全く同じである。

> A <- c(12, 18, 18)
> B <- c(20, 19, 26)
> C <- c(26, 32, 31)
> t.test(A, B)

        Welch Two Sample t-test

data:  A and B 
t = -1.9126, df = 3.969, p-value = 0.1289
mean of x mean of y 
 16.00000  21.66667 

> t.test(A, C)

        Welch Two Sample t-test

data:  A and C 
t = -5.0089, df = 3.978, p-value = 0.007552
mean of x mean of y 
 16.00000  29.66667 

> t.test(B, C)

        Welch Two Sample t-test

data:  B and C 
t = -2.7899, df = 3.898, p-value = 0.05078
mean of x mean of y 
 21.66667  29.66667

 ボンフェローニの修正を加えた有意水準はα = 0.02(ここでは分かりやすいように四捨五入した)であるから、HA vs Bについてはp = 0.1289 > α = 0.02となるので有意ではないといえる。同様にHA vs Cについてはp = 0.007552 < α = 0.02となり有意であり、HB vs Cについてはp = 0.05078 < α = 0.02となり有意ではないといえる。

 以上に示したように、ボンフェローニ法では個々の検定に対する有意水準をコントロールするだけなので非常に簡単な方法である。しかもどのような場合においても適用することができる。しかしながら、このボンフェローニ法ではあまりにも検定に対してあまりにも保守的であるという場合が生じるのが欠点である。

 仮に全体の有意水準をα = 0.05として、帰無仮説の数が10個あった場合は0.05 / 10 = 0.005といったように個々の検定に対して非常に保守的になってしまうことがわかる。要するになかなか有意差が認められなくなってしまうということである。


ホルムの方法

 前述したボンフェローニ法を改良した方法としてホルムの方法というものがある。これもボンフェローニ法と同じく全体の有意水準αから個々の有意水準α'をコントロールするというものである。しかし、ホルムの方法ではそのコントロールの仕方が少し異なる。

 帰無仮説の数がk個であり、全体の有意水準をα = 0.05とした場合、個々の検定における有意水準α'iを次のように定める。

α'1 = α/k, α'2 = α/(k - 1), α'3 = α/(k - 2), ..., α'n = α

 個々の検定によって得られたp値をpiとして、piを小さい順に並べたとき、それをp(1), p(2), p(3), ..., p(n)とすれば、それに対する有意水準は

p(1): α'1, p(2): α'2, p(3): α'3, ..., p(n): α'n

となる。

 再び以下のデータセットを用いて、ホルムの方法を適用して解析してみよう。

A 12, 18, 18
B 20, 19, 26
C 26, 32, 31

 まずは全体の有意水準をα = 0.05としたとき、個々の検定における有意水準α'iを定める。

α'1 = α/k = 0.05/3 = 0.01666667
α'2 = α/(k-1) = 0.05/(3-1) = 0.025
α'3 = α = 0.05

 つぎに個々の検定によって得られたp値を小さい順に並べる。

p(1): {A vs C} = 0.007552, p(2): {B vs C} = 0.05078, p(3): {A vs B} = 0.1289

 結局、

p(1) = 0.007552 < α'1 = 0.01666667
p(2) = 0.05078 > α'2 = 0.025
p(3) = 0.1289 > α'3 = 0.05

となり、有意であるのは{A vs C}のみであることが分かる。

 以上のように、ボンフェローニ法では全ての検定に対する有意水準が0.02であったのに対して、ホルムの方法では個々の検定によって有意水準を工夫してコントロールしている。そのため、ボンフェローニ法では帰無仮説kの個数が大きくなった場合に保守的でありすぎるという問題が改善されている。

<< このページのトップへ



■ 交互作用項が有意である場合の多重比較

 2要因の分散分析において、A要因とB要因の交互作用が認められた場合にはどのような処理を行えばよいだろうか。それを考えるために、まず交互作用とはなんであったかを簡単に復習しておこう。交互作用とはA要因のある水準の下ではB要因のある水準の値が大きくなる(または小さくなる)ということであった。

 例えば、A要因に醤油味、塩味、味噌味という3水準があり、B要因に(麺のゆで加減が)かため、ふつう、やわらかめという3水準があったとして、被験者にラーメンの美味しさを100点満点で評価してもらったとしよう。この場合、A*Bという交互作用項が有意であるという結果が得られた場合どのような解釈ができるだろうか。これは醤油味という水準の下で麺がかためだと評価点が高くなるというようなことである。

 交互作用項が有意である場合は以下の手順に従って多重比較法を適用する必要がある。

手順1 B要因の水準b1, b2, ..., bjそれぞれの下でのA要因の水準a1, a2, ..., ajの効果を調べる。

手順2 A要因の水準a1, a2, ..., ajそれぞれの下でのB要因の水準b1, b2, ..., bjの効果を調べる。

b1 b2 b3
a1 a1b1 = {*, *, *} a1b2 = {*, *, *} a1b3 = {*, *, *}
a2 a2b1 = {*, *, *} a2b2 = {*, *, *} a2b3 = {*, *, *}
a3 a3b1 = {*, *, *} a3b2 = {*, *, *} a3b3 = {*, *, *}

 具体的にいえば、手順1とは{a1b1 vs a2b1, a1b1 vsa3b1, a2b1 vs a3b1}といったように列方向で多重比較していくということである。一方で手順2とは{a1b1 vs a1b2, a1b1 vs a1b3, a1b2 vs a1b3}といったように行方向で多重比較していくということである。

 しかし、全ての水準の組み合わせについて検定する必要などなく、主張する仮説によって比較する組み合わせを考えて(選んで)比較すればよい。場合によっては解釈がややこしくなるので、多重比較を行わない方がよい場合もありえる。

<< このページのトップへ



■ 一般線形モデルにおける多重比較

 ボンフェローニ法でもホルムの方法でも、1要因3水準のデータセットにおいて、{A vs B, A vs C, B vs C}という3つの帰無仮説について検定を行った結果、A vs Cのみが有意であることが認められた。

 分散分析をすることによって得られる結果(分散分析表)は次のようになる。

> A
[1] 12 18 18
> B
[1] 20 19 26
> C
[1] 26 32 31

# 群を表す変数: カテゴリカル型
> group <- c(1, 1, 1, 2, 2, 2, 3, 3, 3)
> group <- as.factor(group)

# 各群のデータを1つのベクトルにまとめる
> y <- c(A, B, C)

# 分散分析モデルの解析
> result <- lm(y ~ group)

# 分散分析表の出力
> summary(aov(result))
            Df  Sum Sq Mean Sq F value   Pr(>F)   
group        2 282.889 141.444  11.573 0.008725 **
Residuals    6  73.333  12.222                    

 これによれば、groupに対するp値はp = 0.008725となっており、これが示す事実は「全体として群間差が認められた」ということである。つまり、観測値yの変動は群を表す変数groupによってうまく説明することができるということである。

 しかし、この結果は全体として差は認められたものの、例えばA vs Bでは差が認められるのかという問いには答えていないことになる。それゆえに多重比較をする必要があるというわけである。

 ところが一般線形モデル、すなわち分散分析モデルとして解析した場合は、そもそも多重比較ということが必要でなくなる。なぜならば、一般線形モデルで扱った場合、分散分析表のみでなく係数表も得られるからである。

> summary(result)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   16.000      2.018   7.927 0.000214 ***
group2         5.667      2.854   1.985 0.094345 .  
group3        13.667      2.854   4.788 0.003039 **

 この係数は各群の効果量であるといえる。つまり、group2の効果量は5.667であり、group3の効果量は13.667である。これはgroup1を0(基準)として考えた場合の効果量であるから、group2はgroup1よりも5.667だけ観測値yに及ぼす影響の度合いが強いということができる(group3も同じ解釈ができる)。

 言い換えると、group2はgroup1よりも5.667だけ平均値が高いということである。つまり、各係数についてのp値が有意水準(例えば0.05)よりも小さければ、統計的に有意にgroup1との差が認められるということである。なお、この係数についてのp値は「推定された係数は0である」という帰無仮説を検定するためのものである。

 したがって、group2の係数についてのp値はp = 0.09なので有意とはいえない(つまり係数は0であるということ)。grou3の係数は有意である(p = 0.00)。これにより、group2の平均値は0であり、group1との差はないと主張できる。group3の平均値は16.000 + 13.667 = 29.667であり、group1との差があると主張できる。そしてこれは先ほどボンフェローニやホルムの方法によって多重比較した結果と一致している。

<< このページのトップへ


参考文献