Embedded Files
Two-way ANOVA（二元配置分散分析）とは、カテゴリ変数が2つのANOVAのことを指す。
ANOVAを含めて、線形モデルの強みの一つは、うまくモデリングすることで、2つ以上の変数が目的変数に与える影響を精度高く検証できることにある。
Two-way ANOVAは、このような線形モデルの長所が発揮される最も基本的な手法である。
要点
要点
Two-way ANOVAには交互作用が存在
要因Aと要因Bとの交互作用とは、「要因Aの効果が、要因Bの水準ごとに異なる」or「要因Bの効果が、要因Aの水準ごとに異なる」or「要因Aと要因Bの効果が相加的でない」こと（全部同値）
周辺平均を理解し、ANOVA表の各行が何をテストしているのか理解すべし
グループ間で繰り返し数が異なる実験では、car::Anova関数推奨
事後検定では、lsmeansまたはemmeansによって、どのグループ間で比較するかを自由に指定可能
例題1
例題1
実験
実験
ある植食性昆虫の幼虫の成長が、寄主植物の水不足によってどのように変化するか調べるため、幼虫を水を遣らない植物（渇水処理：Drought）と水を遣る植物（Control）を与えるグループに分けて飼育した。なお、実験は植物A、植物B、植物Cを使って行った。実験終了後、幼虫の体重を測定したところ、以下のようなデータが得られた（単位はmg）。
> larval.data
Weight Plant.Species Drought.Treatment
1 737.3546 PlantA Control
2 818.3643 PlantA Control
3 716.4371 PlantA Control
4 959.5281 PlantA Control
5 832.9508 PlantA Control
6 717.9532 PlantA Control
7 848.7429 PlantA Control
8 873.8325 PlantA Control
9 857.5781 PlantA Control
10 769.4612 PlantA Control
11 863.0415 PlantA Drought
12 687.6655 PlantA Drought
13 746.5206 PlantA Drought
14 693.5434 PlantA Drought
15 534.7529 PlantA Drought
16 650.2007 PlantA Drought
17 652.6852 PlantA Drought
18 692.8824 PlantA Drought
19 832.0030 PlantA Drought
20 791.5811 PlantA Drought
21 1502.3562 PlantB Control
22 1277.9686 PlantB Control
23 1075.7519 PlantB Control
24 757.0600 PlantB Control
25 1424.9862 PlantB Control
26 1191.0133 PlantB Control
27 1196.7619 PlantB Control
28 1388.7672 PlantB Control
29 1364.2442 PlantB Control
30 1318.7803 PlantB Control
31 486.8381 PlantB Drought
32 479.7311 PlantB Drought
33 555.7571 PlantB Drought
34 544.5331 PlantB Drought
35 444.8995 PlantB Drought
36 443.4004 PlantB Drought
37 529.1666 PlantB Drought
38 561.4826 PlantB Drought
39 491.0123 PlantB Drought
40 570.4886 PlantB Drought
41 782.7080 PlantC Control
42 770.3923 PlantC Control
43 706.7108 PlantC Control
44 520.9583 PlantC Control
45 755.7843 PlantC Control
46 694.9484 PlantC Control
47 685.9784 PlantC Control
48 567.6323 PlantC Control
49 656.9665 PlantC Control
50 737.6147 PlantC Control
51 315.9242 PlantC Drought
52 275.5189 PlantC Drought
53 313.6448 PlantC Drought
54 254.8255 PlantC Drought
55 357.3209 PlantC Drought
56 379.2160 PlantC Drought
57 285.3111 PlantC Drought
58 258.2346 PlantC Drought
59 322.7888 PlantC Drought
60 294.5978 PlantC Drought
要因実験（full factorial experiment）
要因実験（full factorial experiment）
上記の実験では、渇水処理の有無と3つの植物種の全ての組み合わせで、実験が10回繰り返されている。
このことは、以下のようなスクリプトによって確かめられる。
> table(larval.data$Drought.Treatment, larval.data$Plant.Species)
PlantA PlantB PlantC
Control 10 10 10
Drought 10 10 10
> #table関数は複数の要因間での全ての組み合わせ中のデータ数を返す
このように、複数のカテゴリ変数間での全ての組み合わせで処理を行う実験デザインを要因計画（full factorial design）と呼び、このようなデザインで行われる実験を要因実験（full factorial experiment）と呼ぶ。上記の実験では、渇水処理の有無という変数の水準数は2、植物種という変数の水準数は3である。このような実験は、2 x 3 要因計画と呼ばれる。
要因実験の場合、boxplot関数では以下のような記法を用いて、簡単に作図ができる。
boxplot(Weight ~ Drought.Treatment + Plant.Species, larval.data,
cex.axis=0.8, ylab="Larval weight (mg)")
この箱ひげ図を見ると、植物Aでは渇水処理による幼虫の体重の減少はわずかだが、植物BとCでは渇水処理によって幼虫の体重は大きく減少しており、特に植物Bではその傾向が顕著であることがわかる。要因実験の利点の一つは、このように「特定の変数（要因）の影響を、グループ毎に観察できる」ことである。
Two-way ANOVA（二元配置分散分析）は「特定の変数の影響を、グループ毎に観察」というのを統計的に行う手法と言える。
Two-way ANOVAの実践
Two-way ANOVAの実践
さて、one-way ANOVAと同様、two-way ANOVAも正規性・等分散性を仮定している。
ただ、これもone-way ANOVAと同様、特に正規性については最低限分布が一山型であればそれほど気にすることはない。
念のため、まずは見た目で等分散性のチェックのみ行っておく。
par(mfrow=c(1, 2), las=1, lwd=2, family="sans")
boxplot(Weight ~ Drought.Treatment+Plant.Species,
larval.data, ylab="Larval weight", main="raw")
boxplot(log(Weight) ~ Drought.Treatment+Plant.Species,
larval.data, ylab="Larval weight", main="ln")
生データ（raw）では値が小さいほど分散も小さい傾向にあるようだが、対数変換（ln）ではどの箱も同じ程度の分散にみえる。よって、今回は対数変換後のデータで解析を進めることにする。
Two-way ANOVAもone-way ANOVAと同じく、lm関数によって線形モデルを作成することで実行できる。One-way ANOVAと異なるのは、記号*によって検証したい2つの変数が繋げられていることである。
以下のスクリプトを実行すると、以下のような結果を得ることができる。
m <- lm(log(Weight) ~ Drought.Treatment*Plant.Species, larval.data)
anova(m) #Two-way ANOVA
> m <- lm(log(Weight) ~ Drought.Treatment*Plant.Species, larval.data)
> anova(m) #Two-way ANOVA
Analysis of Variance Table
Response: log(Weight)
Df Sum Sq Mean Sq F value Pr(>F)
Drought.Treatment 1 5.5696 5.5696 301.029 < 2.2e-16 ***
Plant.Species 2 3.7801 1.8900 102.154 < 2.2e-16 ***
Drought.Treatment:Plant.Species 2 1.7083 0.8541 46.166 2.044e-12 ***
Residuals 54 0.9991 0.0185
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
One-way ANOVAの出力と異なるのは、結果の行が2つの要因（Drought.TreatmentとPlant.Species）を含んでいることと、その2つの要因をコロン:で繋いだ”交互作用”（Drought.Treatment:Plant.Species）が存在することである。これらの行はそれぞれ、幼虫の体重に対する渇水処理の有無の影響、植物種の影響、渇水処理の有無と植物種との交互作用の影響を表している。
要因Aと要因Bとの交互作用とは、「要因Aの効果が、要因Bの水準ごとに異なる」ことを指す（ただし、これは「要因Aの効果が、要因Bの水準ごとに異なる」または「要因Aと要因Bの効果が相加的でない」ことと同値である）。交互作用に対して、要因AやB単独の効果を、要因AまたはBの主効果と呼ぶ。
今回の結果の場合、箱ひげ図からは「渇水処理の影響が植物の種によって異なっている」というパターンが読み取れる。Drought.Treatment:Plant.Speciesの行においてF2,54 = 46.17, P値が0.05未満であることは、上記のパターンが統計的に有意であるということを意味している。したがって、このtwo-way ANOVAの結果は、「渇水処理の影響は植物の種によって有意に異なる」ことを示している。
ところで、交互作用は普通、"Drought treatment × Plant species"のように表記される。交互作用を×（かける）ではなく:（コロン）で繋ぐのはR言語特有のものである。また、モデル式中での"Drought.Treatment*Plant.Species"という表記は、"Drought.Treatment + Plant.Species + Drought.Treatment:Plant.Species"を省略して書くためのものである。
Two-way ANOVAの場合のANOVA表は右のようになる。
この表を見れば、変動の約半分を渇水処理の有無だけで説明していること（5.57 / (5.57+3.78+1.71+1.00) = 0.46）、残差の分散が0.02とかなり小さいこと、そして、渇水処理の有無、植物の種、そして両者間の交互作用のすべてが有意に幼虫の体重に影響していることが読み取れる。
もちろん、one-way ANOVAの表2のように、sum of squaresやmean squaresを省略しても構わない。
表1. Two-way ANOVA表。Residuals: 残差、df: 自由度、Sum of squares: 平方和、Mean squares: 平均平方、F: F値、P: P値。
なお、同様の解析はパッケージ'car'内の関数Anovaでも行うことができる。
> library(car) #パッケージcarの読み込み
> Anova(m) #Two-way ANOVA
Anova Table (Type II tests)
Response: log(Weight)
Sum Sq Df F value Pr(>F)
Drought.Treatment 5.5696 1 301.029 < 2.2e-16 ***
Plant.Species 3.7801 2 102.154 < 2.2e-16 ***
Drought.Treatment:Plant.Species 1.7083 2 46.166 2.044e-12 ***
Residuals 0.9991 54
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
結果はanova関数を使った場合とまったく同じなのだが、下の補足で説明するように、各グループの繰り返し数がグループ間で異なる場合、anova関数には問題がある。なので、個人的には常にAnova関数を使った方が良いと思う。
Two-way ANOVAの解釈（周辺平均）
Two-way ANOVAの解釈（周辺平均）
ここでは、two-way ANOVAのANOVA表（表1）における各行が、実際には何をテストしているのかについて、少し詳細に解説する。
右の表（表2）は、渇水処理の有無と植物種の組み合わせごとに、対数変換後の幼虫の体重の平均値を計算したものである。
たとえば、無処理区・植物Aの幼虫体重の平均値は6.697である（上の箱ひげ図"ln"と比較してみよう）。
表2. 各要因の全ての組み合わせにおける（対数変換後の）幼虫の体重の平均値と、渇水処理の有無・植物種についての（対数変換後の）幼虫の体重の周辺平均値（太字）。
さて、ここで、ControlおよびDroughtの行の右端にある太字の数字6.779, 6.170に注目してみよう。
この数字は、
6.779 = (6.697+7.115+6.526) / 3
6.170 = (6.563+6.232+5.715) / 3
のように計算された数字で、無処理区および渇水処理区において、植物の種の違いを考慮せずに計算された、幼虫の体重の平均値である（各組み合わせにおけるN数はすべて10で同一なので、上のような計算式で計算することが可能）。このような平均値のことを、渇水処理の有無についての周辺平均と呼ぶ。PlantA, B, C列の下端にある太字の数字も同様に、植物種についての周辺平均である（渇水処理の有無についての周辺平均の計算式を参考に、植物種についての周辺平均を計算し、表2で答え合わせしてみるとよい）。
実は、ANOVA表（表1）において、Drought (D) 行は渇水処理の有無についての周辺平均6.779と6.170間の差を、Plant species (S) 行は植物種についての周辺平均6.630, 6.673, 6.121間の差をテストしている。
雑に言い換えれば、個々の要因に注目した場合の全体的な傾向についてテストしていると言ってもよい。実際に、植物Aでは渇水処理による幼虫の体重の減少はほとんど ないが、植物Bでも植物Cでも渇水処理によって幼虫の体重は大きく減少している。要するに渇水処理で幼虫の体重が減少するような植物の方が多数派であり、周辺平均において無処理よりも渇水処理の方が体重が低いのは、このような全体的な傾向を反映している。
そして、ANOVA表（表1）におけるD × S 行は、渇水処理の有無についての周辺平均に現れている「無処理の方が渇水処理よりも体重が重い（6.779 vs 6.170）」という全体のパターンと、植物A, B, Cそれぞれでの無処理と渇水処理との違い（6.697 vs 6.563, 7.115 vs 6.232, 6.526 vs 5.715）との間にズレが存在するかをテストしている（または、植物種についての周辺平均に現れている「植物Cでは他の植物に比べて体重が軽い（6.630 vs 6.673 vs 6.121）」という全体のパターンと、無処理または渇水処理それぞれでの植物種間の違い（6.697 vs 7.115 vs 6.526, 6.563 vs 6.232 vs 5.715）との間にズレが存在するかをテストしている）。
ここで少し考えなければならない問題は、「Drought (D) 行（主効果）が有意なので、渇水処理は無処理に比べて幼虫の体重が有意に低くなる」と結論付けて良いのかどうかである。
今回、交互作用D × S行が有意なので、渇水処理の影響は植物の種によって有意に異なるはずである。実際、植物Aでは無処理と渇水処理の差は0.13となっていて、他の種に比べてかなり小さく、ほとんど影響がないレベルと言ってよいかもしれない。これは、「渇水処理は無処理に比べて幼虫の体重が有意に低くなる」という結論と矛盾するように見える。
このような立場に立って、統計の教科書では「two-way ANOVAで交互作用が有意だった場合、個々の要因の効果（主効果と呼ばれる）の有意性は意味をなさないことが多い」と解説されることがある。しかし、個々の要因の効果が、その要因に関する全体的な傾向すなわち周辺平均の差をテストしていることを理解しておけば、交互作用が有意でも「個々の要因の効果は意味をなさない」とまでは言えないことが理解できるはずである。
したがって、交互作用が有意であっても、「全体の傾向としては、渇水処理は無処理に比べて幼虫の体重が有意に低くなる」という結論は間違いとは言えない。ただ、交互作用が有意なので、中には渇水処理の影響が小さかったり、場合によっては渇水処理の方が幼虫の体重が重くなっているような（今回は存在しないが）植物種も存在しうるので、それを確認する必要があるだけである。
結局のところ、two-way ANOVAの解釈は、実際のデータがどのようなものなのか、に依存している。
※以上の説明はすべて、繰り返し数が組み合わせ間で等しい場合にのみ当てはまることに注意！！ 繰り返し数が異なれば異なるほど、解釈が難しくなってくる。
回帰係数から見る交互作用
回帰係数から見る交互作用
上記の例で作った線形モデルについては、以下のようにsummary関数で詳細を確認することができる。
> m <- lm(log(Weight) ~ Drought.Treatment*Plant.Species, larval.data)
> summary(m)
Call:
lm(formula = log(Weight) ~ Drought.Treatment * Plant.Species,
data = larval.data)
Residuals:
Min 1Q Median 3Q Max
-0.48545 -0.06652 0.02291 0.09104 0.22313
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.69689 0.04301 155.692 < 2e-16 ***
Drought.TreatmentDrought -0.13390 0.06083 -2.201 0.03201 *
Plant.SpeciesPlantB 0.41799 0.06083 6.871 6.66e-09 ***
Plant.SpeciesPlantC -0.17080 0.06083 -2.808 0.00693 **
Drought.TreatmentDrought:Plant.SpeciesPlantB -0.74913 0.08603 -8.708 7.22e-12 ***
Drought.TreatmentDrought:Plant.SpeciesPlantC -0.67721 0.08603 -7.872 1.59e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.136 on 54 degrees of freedom
Multiple R-squared: 0.9171, Adjusted R-squared: 0.9095
F-statistic: 119.5 on 5 and 54 DF, p-value: < 2.2e-16
Coefficientsの中には6つの行があり、Estimate列にそれぞれの行に対応する回帰係数が並んでいる。
この数字について調べることで、交互作用がどのようなものか考えてみる。
右のスクリプトにより、対数を取った幼虫の体重の平均値を、右図のように各グループごとにプロットできる。
この図を用いて、回帰係数と各グループの平均値がどのような関係にあるのか調べてみよう。
larval.data.mean <- aggregate(log(Weight)~Drought.Treatment+Plant.Species, larval.data, mean)
par(las=1, cex=1.2, family="sans")
plot(c(1, 2, 3.5, 4.5, 6, 7), larval.data.mean[,3], xaxt="n", xlab="",
cex=2, pch=16, ylab="ln(larval weight)", xlim=c(0.5, 7.5), ylim=c(5.5, 7.5))
axis(1, at=c(1, 2, 3.5, 4.5, 6, 7), labels=c("PA,Co", "PA,Dr", "PB,Co", "PB,Dr", "PC,Co", "PC,Dr"))
text(x=c(1, 2, 3.5, 4.5, 6, 7), y=larval.data.mean[,3]+0.2,
labels=as.character(round(larval.data.mean[,3], digits=3)), cex=0.9)
まず、Interceptの係数がPlantAのコントロール区と一致していることが読み取れる（6.697）。
線形モデルでは、Intercept（切片）の値が基準点と想定される。
この例の場合、その基準点は植物A・無処理区となっている。
次に、Drought.TreatmentDroughtの係数は、植物A・無処理区と、植物A・渇水処理区との差と一致している（-0.134）。
したがって、Drought.TreatmentDroughtは、植物Aにおける「渇水処理の影響の大きさ」を表していると解釈できる。
さらに、Plant.SpeciesPlantBの係数は無処理区における植物Aと植物Bの差（0.418）、Plant.SpeciesPlantCの係数は無処理区における植物Aと植物Cの差（-0.171）に一致している。
したがって、Plant.SpeciesPlantBとPlant.SpeciesPlantCは、無処理区における「植物種という要因の影響の大きさ」を表していると解釈できる。
abline(h=larval.data.mean[1,3], lty=2, col="red") #水平線
co <- coef(m) #回帰係数の呼び出しと格納
arrows(c(2, 3.5, 6), co[1], c(2, 3.5, 6), c(co[1]+co[2], co[1]+co[3], co[1]+co[4]),
code=2, length=0.2, col="red", lwd=2) #矢印
次のステップとして、まず植物B、Cでも、渇水処理の影響の大きさは植物Aと同じだと仮定する。
これは、渇水処理と植物種の効果が「相加的」であることに等しい。
この仮定のもとに描かれる幼虫の体重は右図の灰色の点だが、これは実データとはきっかりDrought.TreatmentDrought:Plant.SpeciesPlantBの係数（-0.749）、Drought.TreatmentDrought:Plant.SpeciesPlantCの係数（-0.677）分異なっている。
これらの係数はともに交互作用項の係数であり、それぞれ植物Bにおける渇水処理区での逸脱分、植物Cにおける渇水処理区での逸脱分を表している。
言い換えれば、これらの係数は「植物B、Cでは、渇水処理の影響の大きさが植物Aと比べてどの程度異なるか」を表していると解釈できる。
このように、交互作用の本質は「変数間の相加的な効果から逸脱すること」だと考えられる。
Two-way ANOVAにおける渇水処理と植物種との交互作用の検定結果は、この図で示した逸脱をまとめてテストしていると解釈できる。
（この例を見ればわかるように、回帰係数の大きさ（及びその有意性）はどのグループを基準点(intercept)にするかで変わってしまう。説明変数がカテゴリ変数の場合、変数の有意性を検証するためには、回帰係数でなくANOVAのF検定を用いた方が良い理由はこの辺りにあると思う）
segments(c(3.5, 6), larval.data.mean[c(3, 5),3], c(4.5, 7), larval.data.mean[c(3, 5),3], lty=2) #水平線
points(c(4.5, 7), c(co[1]+co[3]+co[2], co[1]+co[4]+co[2]), pch=16, cex=2, col="grey") #予測値
arrows(c(4.5, 7), c(co[1]+co[3], co[1]+co[4]), c(4.5, 7), c(co[1]+co[3]+co[2], co[1]+co[4]+co[2]),
code=2, length=0.2, lwd=2) #矢印
arrows(c(4.5, 7), c(co[1]+co[3]+co[2], co[1]+co[4]+co[2]), c(4.5, 7), c(co[1]+co[3]+co[2]+co[5], co[1]+co[4]+co[2]+co[6]),
code=2, col="orange", length=0.2, lwd=3) #矢印
Unbalancedなデザインにおけるtwo-way ANOVA
Unbalancedなデザインにおけるtwo-way ANOVA
例題と同じ実験を行ったものの、各グループで10個の繰り返しを取ることができず、結果的に以下のようなデータになったとしよう。
> larval.data.unbalanced
Weight Plant.Species Drought.Treatment
1 737.3546 PlantA Control
2 818.3643 PlantA Control
19 832.0030 PlantA Drought
20 791.5811 PlantA Drought
22 1277.9686 PlantB Control
24 757.0600 PlantB Control
26 1191.0133 PlantB Control
27 1196.7619 PlantB Control
28 1388.7672 PlantB Control
29 1364.2442 PlantB Control
30 1318.7803 PlantB Control
37 529.1666 PlantB Drought
38 561.4826 PlantB Drought
39 491.0123 PlantB Drought
40 570.4886 PlantB Drought
47 685.9784 PlantC Control
50 737.6147 PlantC Control
51 315.9242 PlantC Drought
54 254.8255 PlantC Drought
56 379.2160 PlantC Drought
57 285.3111 PlantC Drought
58 258.2346 PlantC Drought
59 322.7888 PlantC Drought
60 294.5978 PlantC Drought
>
> with(larval.data.unbalanced, table(Drought.Treatment, Plant.Species)) #with関数内では、第1引数でデータを指定することで、変数名だけでベクトルを取り出すことができるので、"larval.data.unbalanced$"みたいなのを省略可能
Plant.Species
Drought.Treatment PlantA PlantB PlantC
Control 2 7 2
Drought 2 4 7
以上のように、2 × 3の組み合わせの各セルで繰り返し数が異なっているデザインを、unbalancedなデザインと呼ぶ。
（例題のように、組み合わせの各セルで繰り返し数が等しいデザインは、balancedなデザインと呼ばれる）
Unbalancedなデザインの場合、anova関数では変数の投入順序によって結果が異なる場合がある。
> m.1 <- lm(log(Weight) ~ Drought.Treatment*Plant.Species, larval.data.unbalanced)
> m.2 <- lm(log(Weight) ~ Plant.Species*Drought.Treatment, larval.data.unbalanced)
>
> anova(m.1) #Two-way ANOVA
Analysis of Variance Table
Response: log(Weight)
Df Sum Sq Mean Sq F value Pr(>F)
Drought.Treatment 1 4.5930 4.5930 206.317 2.656e-11 ***
Plant.Species 2 1.8847 0.9424 42.330 1.566e-07 ***
Drought.Treatment:Plant.Species 2 0.6101 0.3050 13.703 0.0002418 ***
Residuals 18 0.4007 0.0223
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> anova(m.2) #Two-way ANOVA
Analysis of Variance Table
Response: log(Weight)
Df Sum Sq Mean Sq F value Pr(>F)
Plant.Species 2 4.2958 2.14792 96.483 2.396e-10 ***
Drought.Treatment 1 2.1819 2.18191 98.010 1.042e-08 ***
Plant.Species:Drought.Treatment 2 0.6101 0.30505 13.703 0.0002418 ***
Residuals 18 0.4007 0.02226
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
今回の場合、有意性には変化はないが、m.1とm.2では、渇水処理と植物種の平方和(Sum Sq)、平均平方(Mean Sq)、F値が異なっている（1行目に来ている要因の平方和の方が大きくなる傾向にある）。
anova関数で使われている平方和の計算方法（Type I 平方和）では、グループ間で繰り返し数が異なる場合、このようなことが起こる。
これは実践的には結構問題である（変数の投入順序で有意か有意じゃないか変わってしまうかも）。
一方、car::Anovaでは、変数の投入順序で結果が変わるようなことは起こらない。
> Anova(m.1) #Two-way ANOVA
Anova Table (Type II tests)
Response: log(Weight)
Sum Sq Df F value Pr(>F)
Drought.Treatment 2.18191 1 98.010 1.042e-08 ***
Plant.Species 1.88471 2 42.330 1.566e-07 ***
Drought.Treatment:Plant.Species 0.61009 2 13.703 0.0002418 ***
Residuals 0.40072 18
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Anova(m.2) #Two-way ANOVA
Anova Table (Type II tests)
Response: log(Weight)
Sum Sq Df F value Pr(>F)
Plant.Species 1.88471 2 42.330 1.566e-07 ***
Drought.Treatment 2.18191 1 98.010 1.042e-08 ***
Plant.Species:Drought.Treatment 0.61009 2 13.703 0.0002418 ***
Residuals 0.40072 18
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
これは、Anova関数では平方和の計算法のデフォルトがType II 平方和というものになっているからである。
また、Anova関数ではType III平方和というもう一つの計算法も行うことができるが、これも順序には依存しない。
> m.1.1 <- lm(log(Weight) ~ Drought.Treatment*Plant.Species, larval.data.unbalanced,
+ contrasts=list(Drought.Treatment="contr.sum", Plant.Species="contr.sum")) #Type IIIの場合、線形モデル作成の際にcontrastsを"contr.sum"に指定する必要あり
> m.2.1 <- lm(log(Weight) ~ Plant.Species*Drought.Treatment, larval.data.unbalanced,
+ contrasts=list(Drought.Treatment="contr.sum", Plant.Species="contr.sum")) #Type IIIの場合、線形モデル作成の際にcontrastsを"contr.sum"に指定する必要あり
>
> Anova(m.1.1, type="III")
Anova Table (Type III tests)
Response: log(Weight)
Sum Sq Df F value Pr(>F)
(Intercept) 746.87 1 33549.135 < 2.2e-16 ***
Drought.Treatment 1.29 1 57.985 4.913e-07 ***
Plant.Species 1.30 2 29.295 2.187e-06 ***
Drought.Treatment:Plant.Species 0.61 2 13.703 0.0002418 ***
Residuals 0.40 18
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Anova(m.2.1, type="III")
Anova Table (Type III tests)
Response: log(Weight)
Sum Sq Df F value Pr(>F)
(Intercept) 746.87 1 33549.135 < 2.2e-16 ***
Plant.Species 1.30 2 29.295 2.187e-06 ***
Drought.Treatment 1.29 1 57.985 4.913e-07 ***
Plant.Species:Drought.Treatment 0.61 2 13.703 0.0002418 ***
Residuals 0.40 18
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
さて、Type IIとType IIIではともに各項の平方和は順序で変わらない。しかし、Type IIとIIIの間では、主効果の平方和が異なっている。で、結局どっちを使えば良いのか？については議論の分かれるところで、こうすればいつでもOK！とは言えない。が、おおざっぱに言えば以下のようにまとめられる（参考：Hector et al. 2010）。
交互作用が主効果に比べて重要でない（と考えられる）→ Type II
交互作用が主効果と同じかそれ以上に重要だ（と考えられる）→ Type III
事後検定 －どのグループ間で比較するか？－
事後検定 －どのグループ間で比較するか？－
渇水処理の影響が植物の種によって異なるなら、渇水処理はどの植物種で幼虫の体重に有意な影響を与えたと言えるのだろうか？
この問いに答えるための方法の一つは、やはり事後検定である。
Two-way ANOVAでも、one-way ANOVAと同じようにパッケージ'lsmeans'（またはemmeans）を用いることで簡単に事後検定ができる。One-way ANOVAと異なるのは、引数specsに2つの要因を両方指定することである。
m <- lm(log(Weight) ~ Drought.Treatment*Plant.Species, larval.data)
library(lsmeans)
(lsm.DP <- lsmeans(m, specs=c("Drought.Treatment", "Plant.Species")))
#格納しても、()で囲うと中身を出力してくれる
> m <- lm(log(Weight) ~ Drought.Treatment*Plant.Species, larval.data)
>
> library(lsmeans)
>
> (lsm.DP <- lsmeans(m, specs=c("Drought.Treatment", "Plant.Species")))
Drought.Treatment Plant.Species lsmean SE df lower.CL upper.CL
Control PlantA 6.70 0.043 54 6.61 6.78
Drought PlantA 6.56 0.043 54 6.48 6.65
Control PlantB 7.11 0.043 54 7.03 7.20
Drought PlantB 6.23 0.043 54 6.15 6.32
Control PlantC 6.53 0.043 54 6.44 6.61
Drought PlantC 5.71 0.043 54 5.63 5.80
Results are given on the log (not the response) scale.
Confidence level used: 0.95
> #格納しても、()で囲うと中身を出力してくれる
以上のように、specsで2つの要因を指定すると、両者の全ての組み合わせでleast square meanを推定してくれる（これは、実データのグループ毎の平均値と一致する）。あとは、これらのleast square meanをペアで比較してやればよい。
lsm.DPをpairsで囲めば、以下の結果が出力される。
> pairs(lsm.DP, adjust="tukey") #ここではTukey法を指定している
contrast estimate SE df t.ratio p.value
Control,PlantA - Drought,PlantA 0.1339 0.0608 54 2.201 0.2542
Control,PlantA - Control,PlantB -0.4180 0.0608 54 -6.871 <.0001
Control,PlantA - Drought,PlantB 0.4650 0.0608 54 7.645 <.0001
Control,PlantA - Control,PlantC 0.1708 0.0608 54 2.808 0.0715
Control,PlantA - Drought,PlantC 0.9819 0.0608 54 16.142 <.0001
Drought,PlantA - Control,PlantB -0.5519 0.0608 54 -9.073 <.0001
Drought,PlantA - Drought,PlantB 0.3311 0.0608 54 5.444 <.0001
Drought,PlantA - Control,PlantC 0.0369 0.0608 54 0.607 0.9901
Drought,PlantA - Drought,PlantC 0.8480 0.0608 54 13.941 <.0001
Control,PlantB - Drought,PlantB 0.8830 0.0608 54 14.516 <.0001
Control,PlantB - Control,PlantC 0.5888 0.0608 54 9.679 <.0001
Control,PlantB - Drought,PlantC 1.3999 0.0608 54 23.013 <.0001
Drought,PlantB - Control,PlantC -0.2942 0.0608 54 -4.837 0.0002
Drought,PlantB - Drought,PlantC 0.5169 0.0608 54 8.497 <.0001
Control,PlantC - Drought,PlantC 0.8111 0.0608 54 13.334 <.0001
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 6 estimates
このように、6C2 = 15通りのペアワイズの比較を、Tukey法による有意水準の調整をした上でやってくれる。この方法ではかなり詳しく結果が表示される一方で、比較が多すぎてわかりづらいという欠点もある。パッケージ'multcomp'のcld関数では、以下のようにもっと要約した結果も出力できる。
> library(multcomp)
> cld(lsm.DP, adjust="tukey", sort=FALSE, Letters=LETTERS)
Drought.Treatment Plant.Species lsmean SE df lower.CL upper.CL .group
Control PlantA 6.70 0.043 54 6.58 6.81 A
Drought PlantA 6.56 0.043 54 6.45 6.68 A
Control PlantB 7.11 0.043 54 7.00 7.23 B
Drought PlantB 6.23 0.043 54 6.11 6.35 C
Control PlantC 6.53 0.043 54 6.41 6.64 A
Drought PlantC 5.71 0.043 54 5.60 5.83 D
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for 6 estimates
P value adjustment: tukey method for comparing a family of 6 estimates
significance level used: alpha = 0.05
ということで、渇水処理がどの植物種で幼虫の体重に有意な影響を与えたのか？ という問いには、「植物BとCでは処理間に有意な差があるが、植物Aでは処理間に有意な差はない」と結論付けられる。
ただし、上の方法の場合、「植物Aの渇水処理と植物Bのコントロール」とか「植物Bの渇水処理と植物Cの渇水処理」とかといった、「渇水処理がどの植物種で幼虫の体重に有意な影響を与えたのか？」という問いとは関係のない比較も行われている。有意水準の調整を行っていることを考慮すると、このような比較は無駄に「有意になりにくさ」を増しているとも考えられる。渇水処理がどの植物種で幼虫の体重に有意な影響を与えたのか？」という問いに答えるには、植物種ごとにコントロールと渇水処理を比較すれば十分なはずである。
実は、lsmeans関数を上手く使うと、「どことどこの差を比較するか？」をかなり自由に指定可能である。
lsmeans関数には、byという引数があり、これによってleast square meanを以下のように「グループ分け」することが可能である。
> (lsm.DbyP <- lsmeans(m, specs="Drought.Treatment", by="Plant.Species"))
Plant.Species = PlantA:
Drought.Treatment lsmean SE df lower.CL upper.CL
Control 6.70 0.043 54 6.61 6.78
Drought 6.56 0.043 54 6.48 6.65
Plant.Species = PlantB:
Drought.Treatment lsmean SE df lower.CL upper.CL
Control 7.11 0.043 54 7.03 7.20
Drought 6.23 0.043 54 6.15 6.32
Plant.Species = PlantC:
Drought.Treatment lsmean SE df lower.CL upper.CL
Control 6.53 0.043 54 6.44 6.61
Drought 5.71 0.043 54 5.63 5.80
Results are given on the log (not the response) scale.
Confidence level used: 0.95
このような「グループ分けされたleast square mean」をpairs関数で囲むと、グループ毎にペアワイズの比較を行ってくれる。
> pairs(lsm.DbyP, adjust="tukey")
Plant.Species = PlantA:
contrast estimate SE df t.ratio p.value
Control - Drought 0.134 0.0608 54 2.201 0.0320
Plant.Species = PlantB:
contrast estimate SE df t.ratio p.value
Control - Drought 0.883 0.0608 54 14.516 <.0001
Plant.Species = PlantC:
contrast estimate SE df t.ratio p.value
Control - Drought 0.811 0.0608 54 13.334 <.0001
Results are given on the log (not the response) scale.
この方法の場合、有意水準の調整は「グループ内でのみ」行われる（今回の例の場合、グループ内では比較は1つだけなので、結局は調整しない場合と変わらない）。
もし、この3つの比較全体で有意水準の調整を行いたい場合は、各グループの結果を"rbindする"という技が使える。
> rbind(pairs(lsm.DbyP), adjust="tukey")
Plant.Species contrast estimate SE df t.ratio p.value
PlantA Control - Drought 0.134 0.0608 54 2.201 0.0800
PlantB Control - Drought 0.883 0.0608 54 14.516 <.0001
PlantC Control - Drought 0.811 0.0608 54 13.334 <.0001
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates
以上のような方法で、比較したいペアだけに注目するような事後検定を行うことができる。
注目すべきは、植物Aにおける渇水処理とコントロールとの差（のP値）である。すべての組み合わせでの比較を行った場合、植物Aにおける渇水処理とコントロールとの比較ではP = 0.25でまったく有意ではないが、植物種内でのみ比較を行った場合、P値は0.03もしくは0.08となっている。したがって、植物種内でのみ比較を行った場合、植物Aにおいて渇水処理は有意な影響またはmarginally significantな影響を与えた、と解釈することが可能である。
（ただ、比較の個数によってこんな簡単に結論がコロコロ変わっていいのかという疑問は残る。結局のところ、差の有る/無しという二者択一の基準に頼るのではなく、「植物Aでは渇水処理による幼虫の体重の減少は100mg程度だが、植物BとCでは渇水処理によって幼虫の体重は400mg~800mg程度減少した」などのように、定量的な記載を交互作用項の有意性とともに報告するのが望ましいだろう）
2つの要因があるときの、 barplot関数による棒グラフの作図
2つの要因があるときの、 barplot関数による棒グラフの作図
barplot関数で、2つの要因があるときの棒グラフを作成してみよう。
やり方はone-way ANOVAのときとほとんど変わらないが、データを行列の形で用意する点が異なる。
以下のスクリプトにより、各グループのleast square meanを2行3列の行列に変換できる。
lsm.DP.summary <- summary(lsm.DP) #least square meanをsummary関数でdata.frame型に変換
mean <- matrix(lsm.DP.summary$lsmean, ncol=3) #lsmean列を2行3列の行列化
error <- matrix(lsm.DP.summary$SE, ncol=3) #SE列を2行3列の行列化
mean
error
> lsm.DP.summary <- summary(lsm.DP) #least square meanをsummary関数でdata.frame型に変換
> mean <- matrix(lsm.DP.summary$lsmean, ncol=3) #lsmean列を2行3列の行列化
> error <- matrix(lsm.DP.summary$SE, ncol=3) #SE列を2行3列の行列化
> mean
[,1] [,2] [,3]
[1,] 6.696894 7.114888 6.526092
[2,] 6.562991 6.231859 5.714980
> error
[,1] [,2] [,3]
[1,] 0.04301367 0.04301367 0.04301367
[2,] 0.04301367 0.04301367 0.04301367
あとは、この行列を指数関数でback transformして実数スケールに戻し、barplot関数とarrows関数に渡せばよい。
par(las=1, cex=1.5, family="sans", lwd=2)
x.at <- barplot(exp(mean), beside=TRUE, #beside引数をTRUEにすることで、行列データを列ごとに横に並べた棒グラフにする
names.arg=c("Plant A", "Plant B", "Plant C"),
ylab="Larval weight (mg)", ylim=c(0, 1600),
legend.text=c("Control", "Drought")) #legend.textによって凡例をつけることができる
arrows(x.at, exp(mean-error), x.at, exp(mean+error),
angle=90, length=0.1, code=3)
text(c(mean(x.at[1:2]), mean(x.at[3:4]), mean(x.at[5:6])),
c(max(exp(mean+error)[,1]), max(exp(mean+error)[,2]), max(exp(mean+error)[,3]))+120,
labels=c("t54 = 2.2\nP = 0.08", "***", "***"), cex=0.8)
ポイントは引数besideで、これをTRUEにすることで、以下のような棒グラフを描くことができる。
ここでは、植物種内で比較を行い、3つの比較全体でTukey法による有意水準の調整を行った結果をtext関数で付記している。
このように、比較対象を片方の水準内に限定させたい場合には、このようなグループ化した棒グラフが視覚的に有効だろう。
'19 12/27
Google Sites
Report abuse