<< 目次へ
最終更新日 ： 2008.10.7
2要因の分散分析モデル
2要因の分散分析モデルは1要因の分散分析モデルに2つの項が加わる。
Y = X1 + X2 + X1*X2
1要因の分散分析モデルにおいては、1つの従属変数Yにたいして1つのカテゴリカルな独立変数X1が存在した。ところが2要因の分散分析モデルでは、これに加えてもう1つのカテゴリカルな独立変数X2とさらにもう1つのX1*X2という交互作用項が現われる。したがって、2要因以上の分散分析モデルを考える場合、この交互作用項の解釈が重要になってくる。また、3要因以上は結果の解釈がとても複雑になるため、通常は2要因、多くても3要因までに留めておくべきであろう。
ここでは2要因の分散分析モデルについての例題を通して、主に交互作用の解釈について詳しく学んでいくことにしよう。表1のデータはある新薬の投与後に、ラットの体重（Weight）がどう変化したかを記録したものである。実験では平均体重が140gであるラットを18匹用意し、投与する薬の種類（Medicine）と投与量（Dose）という2つの変数が用いられた。Medicineのもつ水準は従来薬（=1）、新薬（=2）の2水準である。Doseのもつ水準は5mg（=1）、10mg（=2）、20mg（=3）の3水準である。
| Medicine | Dose | Before | After | Difference（= Weight） |
| 1 | 1 | 139.1 | 130.7 | 8.4 |
| 1 | 1 | 140.5 | 130.0 | 10.5 |
| 1 | 1 | 140.4 | 131.1 | 9.3 |
| 1 | 2 | 138.8 | 124.9 | 14.0 |
| 1 | 2 | 140.5 | 128.4 | 12.1 |
| 1 | 2 | 143.2 | 128.2 | 15.0 |
| 1 | 3 | 140.5 | 126.8 | 13.7 |
| 1 | 3 | 140.7 | 126.3 | 14.3 |
| 1 | 3 | 140.2 | 127.5 | 12.7 |
| 2 | 1 | 139.1 | 127.0 | 12.1 |
| 2 | 1 | 139.1 | 124.1 | 15.1 |
| 2 | 1 | 139.1 | 125.1 | 14.1 |
| 2 | 2 | 143.5 | 116.0 | 27.6 |
| 2 | 2 | 140.9 | 118.4 | 22.5 |
| 2 | 2 | 140.9 | 114.9 | 26.0 |
| 2 | 3 | 139.4 | 114.9 | 24.5 |
| 2 | 3 | 136.6 | 115.8 | 20.8 |
| 2 | 3 | 140.1 | 113.3 | 26.7 |
さて、表1のデータについてだが、これは経時測定データ（繰り返し測定データ）と呼ばれるものである。経時測定とは同じ個体から実験条件や時間をおいて再びその個体からデータをとることをいう。つまり今回の例だと、18匹のラットそれぞれから2回のデータを測定したことになる。
本来はこういった経時測定データには、一般線形モデルにおける分散分析モデルを適用することはできない（もっと別の特殊なモデルで解析されるべきものである）。しかし、投与前と投与後のデータの差をとることによって、分散分析モデルとして解析することが可能になる。したがって、今回の解析でも従属変数としてDifferenceを用いることにする。なお、従属変数に使ったWeightには表1のDifferenceが代入されており、Difference = Before - Afterである。
まずはRでの解析例を示す前にちょっとした便利な小技を紹介しておく。変数medicineには2つの水準をもつカテゴリカル型のデータが格納されており、doseには3水準のカテゴリカルなデータが格納されている。データをカテゴリカル型として指定するためにはas.factor()を使うのであったが、これをlevels()という関数を利用することによって、各水準のラベルを自由に変更することができる（もちろん、これは分析者の好みによるので、必ずこのようにするべきであるというものではない）。
# 変数の中身を確認してみる
> medicine
[1] 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2
Levels: 1 2
> dose
[1] 1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 3 3 3
Levels: 1 2 3
# medicineをMedicineに代入し、
# levels()を使ってラベルを変更する
> Medicine <- medicine
> levels(Medicine) <- c("Conv.", "New")
# doseをDoseに代入し、
# levels()を使ってラベルを変更する
> Dose <- dose
> levels(Dose) <- c("5mg", "10mg", "20mg")
# 再度、変数の中身を確認してみる
> Medicine
[1] Conv. Conv. Conv. Conv. Conv. Conv. Conv. Conv. Conv. New New New New
[14] New New New New New
Levels: Conv. New
> Dose
[1] 5mg 5mg 5mg 10mg 10mg 10mg 20mg 20mg 20mg 5mg 5mg 5mg 10mg 10mg 10mg 20mg
[17] 20mg 20mg
Levels: 5mg 10mg 20mg
当然、変数Weightには次のようにデータが格納されている。
> Weight [1] 8.3736 10.4541 9.3349 13.9660 12.0735 14.9959 13.7245 14.3414 12.6956 12.1474 [11] 15.0568 14.0570 27.5874 22.5469 26.0365 24.5165 20.7512 26.7186
それでは実際に解析をしてみよう。
# 交互作用項を指定したモデルをlm()で解析する
# それを変数model.einsに代入する
> model.eins <- lm(Weight ~ Medicine + Dose + Medicine*Dose)
# 分散分析表の出力
> summary(aov(model.eins))
Df Sum Sq Mean Sq F value Pr(>F)
Medicine 1 350.76 350.76 95.9731 4.470e-07 ***
Dose 2 232.22 116.11 31.7696 1.607e-05 ***
Medicine:Dose 2 46.08 23.04 6.3043 0.01345 *
Residuals 12 43.86 3.65
# 係数表の出力
> summary(model.eins)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.388 1.104 8.505 2.00e-06 ***
MedicineNew 4.366 1.561 2.797 0.01613 *
Dose10mg 4.291 1.561 2.749 0.01764 *
Dose20mg 4.200 1.561 2.690 0.01965 *
MedicineNew:Dose10mg 7.346 2.207 3.328 0.00603 **
MedicineNew:Dose20mg 6.042 2.207 2.737 0.01803 *
# よりモデル式に忠実な形で係数を出力
> dummy.coef(model.eins)
Full coefficients are
(Intercept): 9.387533
Medicine: Conv. New
0.0000 4.3662
Dose: 5mg 10mg 20mg
0.000000 4.290933 4.199633
Medicine:Dose: Conv.:5mg New:5mg Conv.:10mg New:10mg Conv.:20mg New:20mg
0.000000 0.000000 0.000000 7.345600 0.000000 6.042067
さて、dummy.coef()によって得られた結果を分かりやすく数式に書いておくことにする。
投与後はいずれの条件における、いずれのマウスも体重が減少しているわけである。しかし、符号はプラスになっているので、実際には符号を逆転させて解釈する必要がある。例えば、Medicineの効果はConv.を基準としているので（0としているので）NewはConv.よりも4.36だけ効果量が高いと考えればよい。言い換えれば、新薬を使うと約4gの減量効果があるといえるわけである。
なお、今回はマイナスの符号が出てくると（初学者にとっては）解釈が分かりにくくなると判断したので、符号がプラスになるようにBefore - After = Weightとしたが、通常はAfter - Before = Weightのようにすることが多い。ただこれはどちらでも問題になることはなく、要は符号の解釈をどうするかというだけのことである。
交互作用図
交互作用がみられるかどうかは交互作用図を描くことによって把握することができる。そしてこの交互作用図には5つのタイプがある。まずは今回の実験結果を元に交互作用図を描いてみると、それは図1のようになる。もしかしたら、ある大学の先生から「2つの折れ線が交差している場合は交互作用が存在しているときです」と教わったことがあるかもしれない。しかし、それは誤りである。重要なのは2つの折れ線が交差していることではなく、2つの折れ線が平行ではないということが有意に確認できることである。だから今回の例でも交互作用項は有意であるが、2つの折れ線は交差していないことが確認できる。
図1. 交互作用が存在する場合
では逆に交互作用が存在しない場合にはどのような折れ線が描かれるのだろうか。繰り返しになるが、交互作用が存在する場合は2つの折れ線は平行ではないことが確認できることである。それに対して交互作用が存在しない場合、次の4つの場合が考えられる。
1のようにどちらの要因も体重に影響しない場合、2つの折れ線は水平線になり、かつ2つの折れ線が完全に重なった状態になる。つまり、視覚的には水平線がただ1つ描かれるだけである（図2）。またこれは全平均を表した水平線であることを意味している。したがって、この場合のモデル式はWeight = α（αは定数項、つまり全平均である）となる。
図2. 薬の種類も投与量も影響に差がない場合
2のように薬の種類は体重に影響を与えるが、投与量は影響しない場合、2つの折れ線は平行になる。しかし、投与量によって変化が生じないのでそれぞれの折れ線は水平線の形になる（図3）。つまりこれはMedicineの項だけが有意である場合の交互作用図である。そしてこの場合のモデル式はWeight = Medicineである（WeightはMedicineによってのみ説明される）。
図3. 薬の種類での影響は異なり、投与量では異ならない場合
3のように投与量だけが体重に影響する場合の交互作用図は2つの折れ線が完全に一致した図が描かれる。ただし、図2のとは異なり、投与量が体重に影響するので水平線ではなくなる（図4）。つまり、視覚的には1つの折れ線が図面上に描かれたようになる。そしてこの場合のモデル式はWeight = Dose（WegihtはDoseのみによって説明される）となり、これはDose項のみが優位である場合である。
図4. 投与量での影響は異なり、薬の種類では異ならない場合
最後に紹介するのが4のように薬の種類と投与量が加法的に体重へ影響する場合の交互作用図である。これは2つの主効果が有意であるが、交互作用効果は認められない場合である。したがって、モデル式はWeight = Medicine + Doseとなり、このとき交互作用図は2つの水平線ではない折れ線が平行になる、という形になる（図5）。
図5. 薬の種類と投与量が加法的に影響をもつ場合
交互作用が有意でない場合に注意すべきこと
最初に「特に」重要なことを述べておく必要があるだろう。それはまず最初に交互作用項を含めたモデル解析を行ったとして、その結果、交互作用項が有意であった場合はそれをそのまま述べればよいだろう。だが、もし交互作用項が有意でなかった場合はどうするべきであろうか。このようなときに、交互作用項を含めないモデルを再度あてはめて解析する分析者がいるかもしれない。しかし、それは決定的な誤りである。
なぜなら、このような状況において交互作用項を除いて解析することは、最初から第2種の過誤を犯すことになるからである。つまり、最初に交互作用項を含んだモデルを解析し、その結果として交互作用項が有意でなかったのは、たまたまその実験で検出できなかったに過ぎないかもしれないということである（真のモデルでは交互作用項を含んでいるにもかかわらず、推定されたモデルでは交互作用項が有意でない確率は常に存在している）。
したがって、モデル解析の際には必ず交互作用項を含んだ最も複雑なモデルを解析する必要がある。それでもし交互作用項が有意でなかったとしても、交互作用項を除いた単純なモデルを再適合させて解析してはいけない。交互作用項が有意でなかった場合は、交互作用項を含んだモデルから「有意である主効果の項のみを取り出せばよい」のである。
実際、交互作用が認められない場合は、単純なモデルを再度当てはめて解析しても、主効果の項（2つの独立変数の係数）は同じになるが、、、
# 交互作用項を含むモデル（係数部分を再褐）
> dummy.coef(model.eins)
Full coefficients are
(Intercept): 9.387533
Medicine: Conv. New
0.0000 4.3662
Dose: 5mg 10mg 20mg
0.000000 4.290933 4.199633
Medicine:Dose: Conv.:5mg New:5mg Conv.:10mg New:10mg Conv.:20mg New:20mg
0.000000 0.000000 0.000000 7.345600 0.000000 6.042067
--- --- ---
# 交互作用項を含めない単純なモデルを解析
> model.zwei <- lm(Weight ~ Medicine + Dose)
# 係数の表示
> dummy.coef(model.zwei)
Full coefficients are
(Intercept): 7.156256
Medicine: Conv. New
0.000000 8.828756
Dose: 5mg 10mg 20mg
0.000000 7.963733 7.220667
と、このように交互作用が認められるモデルにおいて、単純なモデルを再適合させると異なる結果になってしまう。そもそも、検定を繰り返し行ってはいけないということを常に肝に銘じておけば、このような間違いはしなくて済むわけである。
Conv.かつ5mgの条件に当てはまる実験データから計算される平均値は、
(8.4 + 10.5 + 9.3) / 3 = 9.4
であるが、交互作用項を含んだモデル（model.eins）は9.387…と正確に予測できているが、交互作用項を含まないモデル（model.zwei）では7.156…とうまく予測できていないことが確認できるだろう。
結果の表示
例えば、薬の種類（従来薬と新薬）ごとにデータをプロットする場合、どのような形で出力するべきであろうか。これについてはいくつか方法があるが、代表的な方法として標準誤差のエラーバーを描く方法と標準偏差のエラーバーを描く方法とがある。結論からいって、好ましい結果の表示の仕方は、図6のように各群のデータをプロットし、「平均値±1.96*標準偏差」で表すことである。
標準誤差は平均値の推定精度を表しているのだが、標準偏差はデータの変動を表している。したがって、図5の見方としては「エラーバーの範囲内に元のデータ（母集団のデータ）の95%が含まれているのだな」ということになる。
また、少し面白い方法として差の標準誤差を表すということもできる。エラーバーを差の標準誤差で描いた場合、2つのエラーバーが重ならなければ、それは2群の平均値の差が認められるということが主張できる。しかし、この方法は好ましくない。その理由の1つは、分析者の目を平均値の比較にばかりに向かせがちになってしまうからである。もう1つの理由は比較する群が多くなるにつれて、第1種の過誤（実は正しいのに仮説を棄却してしまう確率）を犯す確率が高くなるからである。これは多重性の問題とよばれている。
図6. データを平均値±1.96*標準偏差で表す
参考文献
