ANCOVA
（カテゴリ変数と連続変数Aが連続変数Bに与える影響が知りたい）
Embedded Files
ANCOVA（共分散分析）では、カテゴリ変数と連続変数を同時に説明変数に取って分析を行う。
ANCOVAにおける2つの説明変数のうち、連続変数の方を特に「共変量」と呼ぶことがある。
共変量を投入する主な目的は、「カテゴリ変数が目的変数に与える影響を精度高く検証する」ことである。
（なお、ANCOVAもANOVAや単回帰と同様、誤差に正規性と等分散性を仮定しているのだが、本稿ではとりあえずそれは満たされているものとする。）
※事前にTwo-way ANOVAと単回帰の項を読んでおくことをオススメします。
要点
要点
ANCOVAの効能：1. 誤差を小さくする
ANCOVAの効能：2. 見かけの差に騙されない
ANCOVAの効能：3. 真の差をあぶり出す
ANCOVAの効能その1：誤差を小さくする
ANCOVAの効能その1：誤差を小さくする
野外観察１
野外観察１
ある植食性昆虫の産卵選好性を調べるため、植物1株あたり産卵数を植物種ごとに数えた。その際、同時に株あたり葉数も数えた。その結果、以下のようなデータが得られた。
> OvipPrefData1
PlantID Egg.No Plant.Species Leaf.No
1 A_1 26 A 46
2 A_2 18 A 31
3 A_3 20 A 40
4 A_4 31 A 59
5 A_5 12 A 16
6 A_6 17 A 26
7 A_7 15 A 29
8 B_1 19 B 29
9 B_2 28 B 53
10 B_3 24 B 40
11 B_4 35 B 56
12 B_5 32 B 50
13 B_6 14 B 18
14 B_7 33 B 52
15 C_1 27 C 35
16 C_2 28 C 40
17 C_3 13 C 14
18 C_4 22 C 28
19 C_5 29 C 38
20 C_6 43 C 69
21 C_7 42 C 71
平均値の比較（共変量を考慮しない場合）
平均値の比較（共変量を考慮しない場合）
まず、葉の数のデータはなかったものとして、種ごとの卵数の平均値を比較してみる。
図1. 植物種ごとの卵数の平均値±SE。
se <- function(x) sd(x)/sqrt(length(x)) #SEを算出する関数"se"の定義
Egg.mean <- with(OvipPrefData1, tapply(Egg.No, Plant.Species, mean)) #植物種ごとの平均値算出
Egg.se <- with(OvipPrefData1, tapply(Egg.No, Plant.Species, se)) #植物種ごとのSE算出
par(las=1, family="sans", cex=1.2, lwd=2) #図の調整
x.at <- barplot(Egg.mean, xlab="Plant species", ylab="No. eggs", ylim=c(0, 35)) #棒グラフの作成
arrows(x.at, Egg.mean-Egg.se, x.at, Egg.mean+Egg.se, code=3, length=0.1, angle=90) #エラーバーの作成
右のスクリプトにより、左のような棒グラフを描くことができる。
植物種A, B, Cの順に卵数の平均値は高くなる傾向にあるようだが、SEはかなり広めであり、平均値の推定精度が低いことがうかがえる。実際、以下のように植物種を要因としたone-way ANOVAを行っても、植物種の効果は有意ではない。
> library(car)
> simpleANOVA.m <- lm(Egg.No ~ Plant.Species, OvipPrefData1) #one-way ANOVA用線形モデルの作成
> Anova(simpleANOVA.m) #one-way ANOVA
Anova Table (Type II tests)
Response: Egg.No
Sum Sq Df F value Pr(>F)
Plant.Species 319.14 2 2.2104 0.1385
Residuals 1299.43 18
このような場面で検討すべきなのが、葉数を共変量に入れたANCOVAである。
共変量（葉数）と目的変数（卵数）の関係
共変量（葉数）と目的変数（卵数）の関係
ANCOVAを行うためには、まず共変量となる連続変数と目的変数との関係を確認するのが有用である。
図2. 葉数と卵数との関係。
par(las=1, family="sans", cex=1.2) #図の調整
plot(Egg.No ~ Leaf.No, OvipPrefData1, pch=16,
xlab="No. leaves", ylab="No. eggs",
col=(1:3)[unclass(Plant.Species)]) #散布図の作成
legend("topleft", legend=c("Species A", "Species B", "Species C"),
pch=16, col=1:3, cex=0.8) #凡例の作成
simpleRegressionA <- lm(Egg.No ~ Leaf.No, subset(OvipPrefData1, Plant.Species=="A")) #種Aのデータで単回帰
simpleRegressionB <- lm(Egg.No ~ Leaf.No, subset(OvipPrefData1, Plant.Species=="B")) #種Bのデータで単回帰
simpleRegressionC <- lm(Egg.No ~ Leaf.No, subset(OvipPrefData1, Plant.Species=="C")) #種Cのデータで単回帰
abline(simpleRegressionA, col=1) #種Aの回帰線
abline(simpleRegressionB, col=2) #種Bの回帰線
abline(simpleRegressionC, col=3) #種Cの回帰線
右のスクリプトによって、上の散布図を描くことができる。
葉数と卵数との間には極めて明瞭な正の関係が存在している。植食性昆虫の多くは葉に産卵するから、葉が多ければ多いほど卵も多くなるのはいかにもありそうだ。しかし、今回の野外観察はどの植物種に多く産卵するかを知るのが目的であり、その目的のためには葉数の影響はなるべく排除したいところである。同じ葉数に揃えてから、植物種間で卵数を比較できれば良さそうである。
一方、種ごとの回帰線は、種A、B、Cごとに位置が高くなっていることがわかる。ANCOVAの目的は、この3つの回帰線の高さを種間で比較することである。回帰線の高さを比較するという事は、言い換えれば「葉数を種間で揃えてから、種間で卵数を比較する」とも解釈できる。これは、two-way ANOVAが「特定の変数の影響を、グループ毎に観察する」ことと似ている。
ANCOVAの実践
ANCOVAの実践
ANCOVAのやり方は、通常のtwo-way ANOVAとほとんど同じである。ただし、共変量となる連続変数（ここではLeaf.No）が数値ベクトルである必要がある。
is.numeric(OvipPrefData1$Leaf.No) #Leaf.Noが数値ベクトルであることを確認
ANCOVA.m <- lm(Egg.No ~ Plant.Species*Leaf.No, OvipPrefData1) #ANCOVA用線形モデル
Anova(ANCOVA.m) #ANCOVAの実践
> is.numeric(OvipPrefData1$Leaf.No) #Leaf.Noが数値ベクトルであることを確認
[1] TRUE
>
> ANCOVA.m <- lm(Egg.No ~ Plant.Species*Leaf.No, OvipPrefData1) #ANCOVA用線形モデル
> Anova(ANCOVA.m) #ANCOVAの実践
Anova Table (Type II tests)
Response: Egg.No
Sum Sq Df F value Pr(>F)
Plant.Species 117.38 2 17.7621 0.0001108 ***
Leaf.No 1246.00 1 377.0869 4.835e-12 ***
Plant.Species:Leaf.No 3.86 2 0.5841 0.5698060
Residuals 49.56 15
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
葉数という共変量を入れることにより、植物種が卵数に与える効果が有意になっている（F2,15 = 17.76, P < 0.001）。したがって、ANCOVAに基づけば、卵数は植物種間によって有意に異なると結論付けることができる。
なぜ葉数の投入によって植物種の影響が有意になるのかというと、「葉数が卵数に与える影響を統計的に消去することで、植物種が卵数に与える影響を精度高く検証できる」からである。葉数を投入したことにより、誤差の平方和が1299.43（one-way ANOVAのResiduals）から、49.56（ANCOVAのResiduals）と大幅に小さくなっている。これが、統計的な消去の意味するところである。F比の分母は誤差の平方和（を自由度で割ったもの）であり、これが小さいほどF比が大きくなるので、各項が有意になりやすくなる。なので、共変量の投入によって線形モデルの誤差が小さくなるのなら、各項の検証の精度が高くなるのである。
※注意！：カテゴリ変数と連続変数との交互作用の影響が大きそうな場合、ANCOVAによるカテゴリ変数の主効果の解釈には極めて注意を要する。このような困難さを避けるためか、伝統的なANCOVAではカテゴリ変数と連続変数間の交互作用は含めない。
平均値の比較ふたたび（共変量を考慮した場合）
平均値の比較ふたたび（共変量を考慮した場合）
ANCOVAの結果の図示は図2のような感じでも良いのだが、ここでは「共変量を考慮した平均値の比較」として、least square means (estimated marginal means)による棒グラフ的な結果の図示の方法を紹介する。Least square meansの活用は、ANCOVAがやってることを感覚的に掴むためにも役立つと思われる。
まず、least square meansを算出する。これは、パッケージ'lsmeans'内の関数lsmeansを用いて、以下のスクリプトによって算出できる。
library(lsmeans)
(lsm.ANCOVA.m <- lsmeans(ANCOVA.m, specs="Plant.Species")) #ANCOVAモデルからleast square meansを算出
> library(lsmeans)
> (lsm.ANCOVA.m <- lsmeans(ANCOVA.m, specs="Plant.Species")) #ANCOVAモデルからleast square meansを算出
NOTE: Results may be misleading due to involvement in interactions
Plant.Species lsmean SE df lower.CL upper.CL
A 22.0 0.730 15 20.4 23.5
B 25.1 0.700 15 23.6 26.6
C 28.1 0.691 15 26.6 29.5
Confidence level used: 0.95
※注意！："NOTE: Results may be misleading due to involvement in interactions"という警告メッセージが出ているが、これは一般的には無視できない。事前に交互作用の影響が小さいことを確認してから、least square meansの活用に移るべきである。今回は交互作用の影響は小さいと仮定して次のステップに進む。
散布図で理解するleast square means
散布図で理解するleast square means
ANCOVAのように、カテゴリ変数と連続変数が混在したモデルの場合、カテゴリ変数のleast square meansは連続変数の平均値におけるモデルの予測値で表される。言葉で書いてもわかりづらいが、以下の図のように表せば、その意味がわかるだろう。
図3. ANCOVAモデルの予測回帰線と葉数の平均値40における各種の予測値。
以上の散布図の白抜き〇が、lsmeansで算出された値そのものとなっているはずである。これは、least square meansは「葉数を種間で40（葉数の平均値）に揃えて、卵数を種間で比較している」ことを示す。このように、ANCOVAにおけるleast square meansとは、「連続変数を揃えて、カテゴリ変数の水準間で比較する」ためのものだと言うことができる。
#種Aにおける予測用データ作成
DataA <- subset(OvipPrefData1, Plant.Species=="A")
newdata_A <-
data.frame(Plant.Species="A",
Leaf.No=seq(min(DataA$Leaf.No), max(DataA$Leaf.No), by=1))
#種Bにおける予測用データ作成
DataB <- subset(OvipPrefData1, Plant.Species=="B")
newdata_B <-
data.frame(Plant.Species="B",
Leaf.No=seq(min(DataB$Leaf.No), max(DataB$Leaf.No), by=1))
#種Cにおける予測用データ作成
DataC <- subset(OvipPrefData1, Plant.Species=="C")
newdata_C <-
data.frame(Plant.Species="C",
Leaf.No=seq(min(DataC$Leaf.No), max(DataC$Leaf.No), by=1))
par(las=1, family="sans", cex=1.2) #図の調整
plot(Egg.No ~ Leaf.No, OvipPrefData1, pch=16, cex=0.8,
xlab="No. leaves", ylab="No. eggs",
col=(1:3)[unclass(Plant.Species)]) #散布図の作成
legend("topleft", legend=c("Species A", "Species B", "Species C"),
pch=16, col=1:3, cex=0.8) #凡例の作成
lines(newdata_A$Leaf.No, predict(ANCOVA.m, newdata_A), col=1) #種Aの回帰線
lines(newdata_B$Leaf.No, predict(ANCOVA.m, newdata_B), col=2) #種Bの回帰線
lines(newdata_C$Leaf.No, predict(ANCOVA.m, newdata_C), col=3) #種Cの回帰線
Leaf.mean <- mean(OvipPrefData1$Leaf.No) #葉数の平均値は40
newdata_Lmean <- data.frame(Plant.Species=c("A", "B", "C"),
Leaf.No=Leaf.mean) #葉数の平均値における予測用データ作成
points(rep(Leaf.mean, times=3), predict(ANCOVA.m, newdata_Lmean),
col=1:3, pch=1, lwd=2, cex=2) #予測値の上書き
abline(v=Leaf.mean, lty=2) #縦点線
segments(rep(10, times=3), predict(ANCOVA.m, newdata_Lmean),
rep(Leaf.mean, times=3), predict(ANCOVA.m, newdata_Lmean),
lty=2, col=1:3) #横点線
mtext(round(predict(ANCOVA.m, newdata_Lmean), digits=1),
side=2, las=1, line=-1.5,
at=predict(ANCOVA.m, newdata_Lmean)+1.1, col=1:3) #数値の書き込み
次に、算出したleast square meansを使って、最初に描いたのと同じような棒グラフを作成する。以下のようなスクリプトを書くことにより、least square means ± SEの棒グラフを描くことができる。
sum.lsm.ANCOVA.m <- summary(lsm.ANCOVA.m) #lsmeansのデータフレーム化
lsm <- sum.lsm.ANCOVA.m$lsmean #least square meansの格納
ui <- lsm + sum.lsm.ANCOVA.m$SE #least square means+SEの格納（エラーバーの上端）
li <- lsm - sum.lsm.ANCOVA.m$SE #least square means-SEの格納（エラーバーの下端）
par(las=1, family="sans", cex=1.2, lwd=2) #図の調整
x.at <- barplot(lsm, xlab="Plant species", ylab="No. eggs", ylim=c(0, 35),
names.arg=sum.lsm.ANCOVA.m[,1]) #棒グラフの作成
arrows(x.at, ui, x.at, li, code=3, length=0.1, angle=90) #エラーバーの作成
図4. Least square means ± SEで示した、植物種ごとの卵数。
「共変量を考慮した」この棒グラフを、図1と比較してみよう。
棒の高さも図1とは若干異なるが、それよりもはっきり違うのは、SEが図1に比べて大幅に狭くなっていることである。このことは、各植物種ごとの卵数の平均値の推定精度が大きく向上したことを意味する。これが、ANCOVAの効能の一つが「誤差を小さくする」ことを示す図なのである。誤差が小さくなっているので、平均値の差も有意になりやすくなっているのである。
Least square meansによる事後検定
Least square meansによる事後検定
Least square meansを用いれば、ANOVAのときと全く同じように、ANCOVAでも事後検定を行うことができる。
> pairs(lsm.ANCOVA.m, adjust="tukey") #事後検定
contrast estimate SE df t.ratio p.value
A - B -3.09 1.011 15 -3.059 0.0205
A - C -6.09 1.005 15 -6.054 0.0001
B - C -2.99 0.984 15 -3.042 0.0212
P value adjustment: tukey method for comparing a family of 3 estimates
以上のように、Tukey法を用いた多重比較により、全ての種間で卵数に有意な差が認められることがわかる。
ANCOVAの効能その２：見かけの差に騙されない
ANCOVAの効能その２：見かけの差に騙されない
野外観察２
野外観察２
ある植食性昆虫の産卵選好性を調べるため、植物1株あたり産卵数を植物種ごとに数えた。その際、同時に株あたり葉数も数えた。その結果、以下のようなデータが得られた。
> OvipPrefData2
PlantID Egg.No Plant.Species Leaf.No
1 A_1 3 A 6
2 A_2 12 A 17
3 A_3 10 A 4
4 A_4 5 A 10
5 A_5 12 A 19
6 A_6 7 A 7
7 A_7 9 A 8
8 B_1 26 B 47
9 B_2 33 B 49
10 B_3 28 B 51
11 B_4 32 B 56
12 B_5 29 B 46
13 B_6 24 B 45
14 B_7 24 B 49
15 C_1 39 C 61
16 C_2 31 C 69
17 C_3 38 C 65
18 C_4 30 C 53
19 C_5 42 C 72
20 C_6 38 C 68
21 C_7 48 C 77
平均値の比較（共変量を考慮しない場合）
平均値の比較（共変量を考慮しない場合）
観察２は、データの違いがもたらす結果を分かりやすく示すために、観察１と目的、デザイン共に全く同じにしている。
なので、今回もまずは同じように葉数を無視して単純に植物種ごとに平均値を比較してみる。
図5. 植物種ごとの卵数の平均値±SE。
se <- function(x) sd(x)/sqrt(length(x)) #事前に定義していればこの行は不要
Egg.mean <- with(OvipPrefData2, tapply(Egg.No, Plant.Species, mean)) #植物種ごとの平均値算出
Egg.se <- with(OvipPrefData2, tapply(Egg.No, Plant.Species, se)) #植物種ごとのSE算出
par(las=1, family="sans", cex=1.2, lwd=2) #図の調整
x.at <- barplot(Egg.mean, xlab="Plant species", ylab="No. eggs", ylim=c(0, 45)) #棒グラフの作成
arrows(x.at, Egg.mean-Egg.se, x.at, Egg.mean+Egg.se, code=3, length=0.1, angle=90) #エラーバーの作成
平均値を比較すると、卵数は植物種A, B, Cの順にはっきりと多くなっている。実際、植物種のみを説明変数にしてone-way ANOVAを行うと、植物種間で有意な差が検出される。
> library(car) #事前に呼び出していればここは不要
> simpleANOVA.m <- lm(Egg.No ~ Plant.Species, OvipPrefData2) #one-way ANOVA用線形モデルの作成
> Anova(simpleANOVA.m) #one-way ANOVA
Anova Table (Type II tests)
Response: Egg.No
Sum Sq Df F value Pr(>F)
Plant.Species 3200.4 2 75.913 1.688e-09 ***
Residuals 379.4 18
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
…さて、これで「卵数は植物種間で異なる→植物種に対する産卵選好性がある！」と結論付ける人は見かけの差に騙されている。結論を下すためには、種間で葉の枚数に差がないだろうか？また、もし種間で葉の枚数がいっしょなら？と考えなければならない。
共変量（葉数）と目的変数（卵数）の関係
共変量（葉数）と目的変数（卵数）の関係
そこで、葉数と卵数の関係を調べてみる。
図6. 葉数と卵数との関係。
par(las=1, family="sans", cex=1.2) #図の調整
plot(Egg.No ~ Leaf.No, OvipPrefData2, pch=16,
xlab="No. leaves", ylab="No. eggs",
col=(1:3)[unclass(Plant.Species)]) #散布図の作成
legend("topleft", legend=c("Species A", "Species B", "Species C"),
pch=16, col=1:3, cex=0.8) #凡例の作成
simpleRegressionA <- lm(Egg.No ~ Leaf.No, subset(OvipPrefData2, Plant.Species=="A")) #種Aのデータで単回帰
simpleRegressionB <- lm(Egg.No ~ Leaf.No, subset(OvipPrefData2, Plant.Species=="B")) #種Bのデータで単回帰
simpleRegressionC <- lm(Egg.No ~ Leaf.No, subset(OvipPrefData2, Plant.Species=="C")) #種Cのデータで単回帰
abline(simpleRegressionA, col=1) #種Aの回帰線
abline(simpleRegressionB, col=2) #種Bの回帰線
abline(simpleRegressionC, col=3) #種Cの回帰線
すると、植物種ごとの葉数はA, B, Cの順に多くなっていること、さらに卵の数は植物種に関わらず葉数が増えるほど多くなっており、その回帰線は種間でほとんど重なってしまっていることが分かる。このことは、もし葉の数が種間で同じになれば（A種の大きめの個体と、B,C種の小さめの個体で比べることをイメージすべし）、卵の数はほとんど種間で差がなくなるだろうことを意味している。
ANCOVAの実践
ANCOVAの実践
そこで、野外観察１と同じように、葉数をモデルに投入したANCOVAを実行すると、以下のような結果が得られる。
> ANCOVA.m <- lm(Egg.No ~ Plant.Species*Leaf.No, OvipPrefData2) #ANCOVA用線形モデル
> Anova(ANCOVA.m) #ANCOVAの実践
Anova Table (Type II tests)
Response: Egg.No
Sum Sq Df F value Pr(>F)
Plant.Species 2.901 2 0.1066 0.899536
Leaf.No 169.710 1 12.4769 0.003017 **
Plant.Species:Leaf.No 5.689 2 0.2091 0.813619
Residuals 204.029 15
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
葉数という共変量を入れることにより、植物種が卵数に与える効果が有意でなくなった（F2,15 = 0.11, P = 0.90）。したがって、ANCOVAに基づけば、植物種間で卵数に有意な差はないと結論付けることができる。これは要するに、「卵数の平均値が植物種間で異なったのは、単に葉の数が多いほど卵数が多いから」だろう、ということを意味している。
各項の平方和を見ると、one-way ANOVAの時と比べて、誤差項はほとんど減っていない一方、植物種の主効果の平方和が大幅に減っている（3200.4→2.901）。植物種のF比の分子は平方和（を自由度で割ったもの）なので、ここが減ったことは植物種の主効果が有意になりにくくなったことを示している。共変量の投入により、植物種の見かけの影響から、葉数の影響が統計的に消去されたのである。
交絡要因
交絡要因
野外観察２における葉数のような要因を「交絡要因」と呼ぶことがある。これは、注目している要因（ここでは植物種）と"相関"することで、注目する要因について見かけの影響をもたらす要因のことである。仮に、植物種Aに1, Bに2, Cに3という番号を振り、葉数との関係を散布図でプロットしてみる。
図7. 葉数と植物種との関係。
par(las=1, family="sans") #図の調整
plot(as.numeric(Plant.Species) ~ Leaf.No, OvipPrefData2, pch=1, cex=1.2, lwd=2,
xlab="No. leaves", ylab="Plant species", yaxt="n") #散布図の作成
axis(2, labels=c("1", "2", "3"), at=1:3) #y軸ラベルの作成
両者は明らかに相関していることがわかる。一般に、2つの説明変数が強く相関している場合、どちらの変数が目的変数に真に影響しているか結論を下すことは統計解析からは不可能である。これは、これら2つの説明変数が、目的変数への影響を共有してしまっているからである。言えるのは、共有部分を取り除いた後も個々の説明変数単独の影響が残るかどうかだけである。説明変数間に相関が存在することを「多重共線性」と呼ぶ。
さて、野外観察２におけるANCOVAの線形モデルはまさに多重共線性をはらんだモデルであるが、多重共線性を避けるため、葉数はモデルから除いた方が良いと提案したら、ANCOVAの説明をここまで読んできた人なら誰でもおかしいと思うだろう。これは、注目したい要因と相関する別の要因（共変量）を積極的にモデルに含めるべき場面も存在する例の一つである。
平均値の比較ふたたび（共変量を考慮した場合）
平均値の比較ふたたび（共変量を考慮した場合）
以上を踏まえ、least square meansによって植物種ごとの卵数の平均値を表すことを試みる。
図8. Least square means ± SEで示した、植物種ごとの卵数。
library(lsmeans) #事前に呼び出していればここは不要
lsm.ANCOVA.m <- lsmeans(ANCOVA.m, specs="Plant.Species") #ANCOVAモデルからleast square meansを算出
sum.lsm.ANCOVA.m <- summary(lsm.ANCOVA.m) #lsmeansのデータフレーム化
lsm <- sum.lsm.ANCOVA.m$lsmean #least square meansの格納
ui <- lsm + sum.lsm.ANCOVA.m$SE #least square means+SEの格納（エラーバーの上端）
li <- lsm - sum.lsm.ANCOVA.m$SE #least square means-SEの格納（エラーバーの下端）
par(las=1, family="sans", cex=1.2, lwd=2) #図の調整
x.at <- barplot(lsm, xlab="Plant species", ylab="No. eggs", ylim=c(0, 45),
names.arg=sum.lsm.ANCOVA.m[,1]) #棒グラフの作成
arrows(x.at, ui, x.at, li, code=3, length=0.1, angle=90) #エラーバーの作成
このように、産卵された卵数には植物種間でほとんど差がない、という図が描かれる。
注意すべきはエラーバー（SE）も図5に比べて大幅に広くなっていることである。これは、平均値の推定精度がかなり悪くなっていることを意味する。図6の散布図を見ればわかるように、このleast square meansは”外挿”であり、推定精度が低いのは当然である。もし、葉数30~50くらいにおける種A, B, Cそれぞれのデータを集めることができれば、もしかしたら卵数には種間で差があるかもしれない。したがって、このようなデータにおけるANCOVAは、真の差を導き出すものと言うよりは、「見かけの差から結論を下すのを防ぐ」ものと考えた方が良いかもしれない。
ANCOVAの効能その３：真の差をあぶり出す
ANCOVAの効能その３：真の差をあぶり出す
野外観察３
野外観察３
ある植食性昆虫の産卵選好性を調べるため、植物1株あたり産卵数を植物種ごとに数えた。その際、同時に株あたり葉数も数えた。その結果、以下のようなデータが得られた。
> OvipPrefData3
PlantID Egg.No Plant.Species Leaf.No
1 A_1 26 A 23
2 A_2 38 A 41
3 A_3 20 A 20
4 A_4 30 A 31
5 A_5 42 A 44
6 A_6 29 A 25
7 A_7 25 A 26
8 B_1 35 B 45
9 B_2 34 B 48
10 B_3 40 B 51
11 B_4 50 B 60
12 B_5 38 B 44
13 B_6 33 B 41
14 B_7 42 B 49
15 C_1 39 C 49
16 C_2 44 C 59
17 C_3 37 C 54
18 C_4 28 C 38
19 C_5 48 C 62
20 C_6 43 C 57
21 C_7 49 C 69
平均値の比較（共変量を考慮しない場合）
平均値の比較（共変量を考慮しない場合）
今回もまず、葉の数のデータはなかったものとして、種ごとの卵数の平均値を比較してみる。
図9. 植物種ごとの卵数の平均値±SE。
se <- function(x) sd(x)/sqrt(length(x)) #事前に定義していればこの行は不要
Egg.mean <- with(OvipPrefData3, tapply(Egg.No, Plant.Species, mean)) #植物種ごとの平均値算出
Egg.se <- with(OvipPrefData3, tapply(Egg.No, Plant.Species, se)) #植物種ごとのSE算出
par(las=1, family="sans", cex=1.2, lwd=2) #図の調整
x.at <- barplot(Egg.mean, xlab="Plant species", ylab="No. eggs", ylim=c(0, 45)) #棒グラフの作成
arrows(x.at, Egg.mean-Egg.se, x.at, Egg.mean+Egg.se, code=3, length=0.1, angle=90) #エラーバーの作成
SEがやや広めではあるが、種A, B, Cの順に卵数が多くなる傾向が見て取れる。このデータに対して、植物種を説明変数に取ったone-way ANOVAを行うと、やはり有意となる。
> library(car) #事前に呼び出していればここは不要
> simpleANOVA.m <- lm(Egg.No ~ Plant.Species, OvipPrefData3) #one-way ANOVA用線形モデルの作成
> Anova(simpleANOVA.m) #one-way ANOVA
Anova Table (Type II tests)
Response: Egg.No
Sum Sq Df F value Pr(>F)
Plant.Species 484.95 2 4.9954 0.01881 *
Residuals 873.71 18
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
…この結果をもって、「産卵選好性は植物種A, B, Cの順に高くなる」と結論付けるのはやはり早計である。ここでもやはり、葉数が卵数に与える影響について考えなければならない。
共変量（葉数）と目的変数（卵数）の関係
共変量（葉数）と目的変数（卵数）の関係
図10. 葉数と卵数との関係。
par(las=1, family="sans", cex=1.2) #図の調整
plot(Egg.No ~ Leaf.No, OvipPrefData3, pch=16,
xlab="No. leaves", ylab="No. eggs",
col=(1:3)[unclass(Plant.Species)]) #散布図の作成
legend("topleft", legend=c("Species A", "Species B", "Species C"),
pch=16, col=1:3, cex=0.8) #凡例の作成
simpleRegressionA <- lm(Egg.No ~ Leaf.No, subset(OvipPrefData2, Plant.Species=="A")) #種Aのデータで単回帰
simpleRegressionB <- lm(Egg.No ~ Leaf.No, subset(OvipPrefData2, Plant.Species=="B")) #種Bのデータで単回帰
simpleRegressionC <- lm(Egg.No ~ Leaf.No, subset(OvipPrefData2, Plant.Species=="C")) #種Cのデータで単回帰
abline(simpleRegressionA, col=1) #種Aの回帰線
abline(simpleRegressionB, col=2) #種Bの回帰線
abline(simpleRegressionC, col=3) #種Cの回帰線
葉数と卵数の間に明瞭な正の関係があること、葉数が植物種A, B, Cの順に多くなっていることは野外観察２と同じである。ここで注目すべきは、回帰線の高さが単純な平均値の大小関係とは逆に、植物種C, B, Aの順に高くなっていることである。このことは、図9の大小関係は葉数が交絡要因となって生じた見かけのものであること、大きめの植物種Aと小さめの植物種B, Cを比較すれば、植物種A上の卵数が3種中もっとも多いだろうことを示唆している。
ANCOVAの実践
ANCOVAの実践
葉数をモデルに投入したANCOVAを実行すると、以下のような結果が得られる。
> ANCOVA.m <- lm(Egg.No ~ Plant.Species*Leaf.No, OvipPrefData3) #ANCOVA用線形モデル
> Anova(ANCOVA.m) #ANCOVAの実践
Anova Table (Type II tests)
Response: Egg.No
Sum Sq Df F value Pr(>F)
Plant.Species 91.50 2 8.0764 0.004164 **
Leaf.No 783.36 1 138.2854 5.698e-09 ***
Plant.Species:Leaf.No 5.38 2 0.4746 0.631157
Residuals 84.97 15
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
植物種の主効果はF2,15 = 8.08, P < 0.01で有意であり、植物種間で卵数に有意差があると言える。しかし、ここでの問題は、「種間でどのような差があるのか？」である。それを端的に表すために、least square meansによる図示が極めて有効である。
平均値の比較ふたたび（共変量を考慮した場合）
平均値の比較ふたたび（共変量を考慮した場合）
卵数が種間でどのように異なっているか検証するため、上のANCOVAモデルを用いて、least square meansによって植物種ごとの卵数の平均値を表してみる。
図11. Least square means ± SEで示した、植物種ごとの卵数。
library(lsmeans) #事前に呼び出していればここは不要
lsm.ANCOVA.m <- lsmeans(ANCOVA.m, specs="Plant.Species") #ANCOVAモデルからleast square meansを算出
sum.lsm.ANCOVA.m <- summary(lsm.ANCOVA.m) #lsmeansのデータフレーム化
lsm <- sum.lsm.ANCOVA.m$lsmean #least square meansの格納
ui <- lsm + sum.lsm.ANCOVA.m$SE #least square means+SEの格納（エラーバーの上端）
li <- lsm - sum.lsm.ANCOVA.m$SE #least square means-SEの格納（エラーバーの下端）
par(las=1, family="sans", cex=1.2, lwd=2) #図の調整
x.at <- barplot(lsm, xlab="Plant species", ylab="No. eggs", ylim=c(0, 45),
names.arg=sum.lsm.ANCOVA.m[,1]) #棒グラフの作成
arrows(x.at, ui, x.at, li, code=3, length=0.1, angle=90) #エラーバーの作成
このように、植物種間における卵数の大小関係は図9と真逆になっている。したがって、ANCOVAに基づけば、「この植食性昆虫の産卵選好性は、植物種Aに対して最も高い」と結論付けることができる。これは、以下のように事後検定を行うことでも補強できる。
> pairs(lsm.ANCOVA.m, adjust="tukey") #事後検定、種Aの卵数は他の種に比べて有意に多い
contrast estimate SE df t.ratio p.value
A - B 6.08 2.08 15 2.922 0.0268
A - C 8.20 2.26 15 3.626 0.0066
B - C 2.12 1.76 15 1.207 0.4674
P value adjustment: tukey method for comparing a family of 3 estimates
以上のように、ANCOVAは、単に見かけの差に騙されないようにするだけでなく、時には見かけの差に惑わされずに真の差をあぶり出すことも可能である。なお、このようなケースは、注目するカテゴリ変数と共変量が相関しているからこそ生じていることにも留意が必要である。
'20 2/19 記
Google Sites
Report abuse