ANCOVA
1 高校生ユニット 重回帰分析から分散分析へ
1.1 ２変量回帰分析.
真の回帰式として、
\(z=1+0.2 x-0.3
y\) を考え、誤差変動として \(\varepsilon_i \sim N(0,~0.1^2) ~:~ i.i.d.\)
を加えた
\(z_i=1+0.2 x_i-0.3
y_i+\varepsilon_i\) を回帰モデルとします。
xy平面の点として、次の9つの点を固定します。
library(tidyr)
data <- crossing(x = 1:3, y = 1:3);#data
data <- rbind(data,data,data)#３回繰り返す
data[1:9,]| x | y |
|---|---|
| 1 | 1 |
| 1 | 2 |
| 1 | 3 |
| 2 | 1 |
| 2 | 2 |
| 2 | 3 |
| 3 | 1 |
| 3 | 2 |
| 3 | 3 |
\((x_i,y_i)\) に対し \(z_i=1+0.2 x_i-0.3
y_i+\varepsilon_i\) により \(z_i\)
を３個ずつ出力させ、3Dプロットに表示します。
次に Rの線形回帰モデル lm（linear model）によって
回帰係数を求め回帰平面を表示します。
library(tidyr)
library(gapminder)
options(rgl.printRglwidget = TRUE)
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)
x <- data$x
y <- data$y
z <- 0.2*x - 0.3*y+ rnorm(27, sd = 0.1)#誤差分散 0.1^2
plot3d(x = x, y = y, z = z,type = "s", col = "blue", size = 1)
#open3d()
fit <- lm(z ~ x + y)
f <- fit$ fitted.values;#予測値
a <- fit$ coefficients;#回帰係数
e <- fit$ residuals;#e#残差
se <- sum(e^2);se#残差平方和## [1] 0.2095098s <- sqrt(se/(27-3));s#残差標準偏差## [1] 0.09343219ff<- (var(f)*26/2)/(se/24);ff#F値## [1] 152.345summary(fit)##
## Call:
## lm(formula = z ~ x + y)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.181329 -0.060407 -0.008652 0.061300 0.147369
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.005626 0.064832 0.087 0.932
## x 0.206615 0.022022 9.382 1.68e-09 ***
## y -0.324157 0.022022 -14.720 1.64e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09343 on 24 degrees of freedom
## Multiple R-squared: 0.927, Adjusted R-squared: 0.9209
## F-statistic: 152.3 on 2 and 24 DF, p-value: 2.297e-14plot3d(x, y, z, type = "s", col = "blue", size = 1)
coefs <- coef(fit)
a <- coefs["x"]
b <- coefs["y"]
c <- -1
d <- coefs["(Intercept)"]
#ax+by＋cz+d=0で表される平面
planes3d(a, b, c, d, alpha = 0.5)
1.2 aov分散分析
res <- aov(z ~ factor(x) + y)
summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## factor(x) 2 0.7705 0.3853 42.73 1.80e-08 ***
## y 1 1.8914 1.8914 209.77 4.76e-13 ***
## Residuals 23 0.2074 0.0090
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction.plot(x,y,z, type="b", pch=c(1,2))
2 市原p220
x <- c(54,50,39,44,35,41,32,46,38,30,43,39,49,42,35,31,36,34,27)
y <- c(rep(1,10),rep(2,9))
z <- c(17,15,15,14,13,12,11,10,10,9,22,21,19,18,16,15,13,11,9)
data <- data.frame(x=x,y=y,z=z)
res <- aov(z ~ x * factor(y))
summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 67.24 67.24 10.646 0.00524 **
## factor(y) 1 95.56 95.56 15.131 0.00145 **
## x:factor(y) 1 17.62 17.62 2.789 0.11562
## Residuals 15 94.74 6.32
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction.plot(x,y,z, type="b", pch=c(1,2))
par(las=1, family="sans", cex=1.2) #図の調整
plot(z ~ x, data, pch=16,
xlab="No. leaves", ylab="No. eggs",col=c(1:2))
simpleRegression1 <- lm(z ~ x, subset(data, y==1)) #単回帰
simpleRegression2 <- lm(z ~ x, subset(data, y==2)) #単回帰
abline(simpleRegression1, col=1)
abline(simpleRegression2, col=2)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## someANCOVA.m <- lm(z ~ x * factor(y), data) #ANCOVA用線形モデル
Anova(ANCOVA.m) | Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| x | 108.04386 | 1 | 17.106606 | 0.0008797 |
| factor(y) | 95.56324 | 1 | 15.130547 | 0.0014511 |
| x:factor(y) | 17.61743 | 1 | 2.789371 | 0.1156206 |
| Residuals | 94.73871 | 15 | NA | NA |
交互作用のF値は2.78となっており、市原p222の値 1.67 の二乗である。
x <- c(1:19)
y <- c(rep(1,10),rep(2,9))
z <- c(17,10,10,9,22,21,19,18,16,15,13,11,9,15,15,14,13,12,11)
#res <- aov(z ~ x + y)
#summary(res)
interaction.plot(x,y,z, type="b", pch=c(1,2))
連続変数の方を特に「共変量」
3 橋下
3.1 データ.
OvipPrefData1 <- read.csv("OvipPrefData12.csv")
str(OvipPrefData1)## 'data.frame': 21 obs. of 5 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ PlantID : chr "A_1" "A_2" "A_3" "A_4" ...
## $ Egg.No : int 26 18 20 31 12 17 15 19 28 24 ...
## $ Plant.Species: chr " A " " A " " A " " A " ...
## $ Leaf.No : int 46 31 40 59 16 26 29 29 53 40 ...head(OvipPrefData1)| X | 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 |
x <- OvipPrefData1$Leaf.No
y <- OvipPrefData1$Plant.Species
z <- OvipPrefData1$Egg.No
res <- aov(z ~ x + y)
summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 1447.8 1447.8 460.69 9.41e-14 ***
## y 2 117.4 58.7 18.68 5.12e-05 ***
## Residuals 17 53.4 3.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction.plot(x,y,z, type="b", pch=c(1,2))
x <- OvipPrefData1$Leaf.No
y <- OvipPrefData1$Plant.Species
z <- OvipPrefData1$Egg.No
res <- aov(z ~ x + y)
summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 1447.8 1447.8 460.69 9.41e-14 ***
## y 2 117.4 58.7 18.68 5.12e-05 ***
## Residuals 17 53.4 3.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction.plot(x,y,z, type="b", pch=c(1,2))
3.2 平均値の比較.
#SEを算出する関数"se"の定義
se <- function(x) sd(x)/sqrt(length(x))
#植物種ごとの平均値算出
Egg.mean <- with(OvipPrefData1, tapply(Egg.No, Plant.Species, mean))
Egg.mean## A B C
## 19.85714 26.42857 29.14286#植物種ごとのSE算出
Egg.se <- with(OvipPrefData1, tapply(Egg.No, Plant.Species, se))
Egg.se## A B C
## 2.482592 2.950844 4.008494#図の調整
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) 
3.3 one-way-anova.
simpleANOVA.m <- lm(Egg.No ~ Plant.Species, OvipPrefData1) #one-way ANOVA用線形モデルの作成
Anova(simpleANOVA.m) #one-way ANOVA| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| Plant.Species | 319.1429 | 2 | 2.210422 | 0.1385431 |
| Residuals | 1299.4286 | 18 | NA | NA |
3.4 共変量（葉数）と目的変数（卵数）の関係.
library(ggplot2)
library(tidyverse)
attach(OvipPrefData1)
OvipPrefData1$Egg.No## [1] 26 18 20 31 12 17 15 19 28 24 35 32 14 33 27 28 13 22 29 43 42OvipPrefData1$Leaf.No## [1] 46 31 40 59 16 26 29 29 53 40 56 50 18 52 35 40 14 28 38 69 71Plant.Species## [1] " A " " A " " A " " A " " A " " A "
## [7] " A " " B " " B " " B " " B " " B "
## [13] " B " " B " " C " " C " " C " " C "
## [19] " C " " C " " C "col=c(1,2,3);col## [1] 1 2 3par(las=1, family="sans", cex=1.2) #図の調整
plot(Egg.No ~ Leaf.No, OvipPrefData1, pch=16,
xlab="No. leaves", ylab="No. eggs",col=c(1:3))
plot(Egg.No ~ Leaf.No, OvipPrefData1, pch=16,
xlab="No. leaves", ylab="No. eggs",
col=c(1:3)) #散布図の作成
legend("topleft", legend=c("Species A", "Species B", "Species C"), pch=16, col=1:3, cex=0.8) #凡例の作成
head(OvipPrefData1)| X | 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 |
str(OvipPrefData1)## 'data.frame': 21 obs. of 5 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ PlantID : chr "A_1" "A_2" "A_3" "A_4" ...
## $ Egg.No : int 26 18 20 31 12 17 15 19 28 24 ...
## $ Plant.Species: chr " A " " A " " A " " A " ...
## $ Leaf.No : int 46 31 40 59 16 26 29 29 53 40 ...head(subset(OvipPrefData1, Plant.Species==" A "))| X | 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 |
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の回帰線
4 ANCOVAの実践
is.numeric(OvipPrefData1$Leaf.No) #Leaf.Noが数値ベクトルであることを確認## [1] TRUEANCOVA.m <- lm(Egg.No ~ Plant.Species*Leaf.No, OvipPrefData1) #ANCOVA用線形モデル
Anova(ANCOVA.m) #ANCOVAの実践| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| Plant.Species | 117.382559 | 2 | 17.7621480 | 0.0001108 |
| Leaf.No | 1246.004241 | 1 | 377.0868849 | 0.0000000 |
| Plant.Species:Leaf.No | 3.859987 | 2 | 0.5840873 | 0.5698060 |
| Residuals | 49.564343 | 15 | NA | NA |
5 平均値の比較ふたたび（共変量を考慮した場合）
library(lsmeans)## Loading required package: emmeans## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.(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#種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) #葉数の平均値における予測用データ作成rep(Leaf.mean, times=3)## [1] 40 40 40newdata_Lmean| Plant.Species | Leaf.No |
|---|---|
| A | 40 |
| B | 40 |
| C | 40 |
ANCOVA.m##
## Call:
## lm(formula = Egg.No ~ Plant.Species * Leaf.No, data = OvipPrefData1)
##
## Coefficients:
## (Intercept) Plant.Species B
## 3.96862 0.01509
## Plant.Species C Leaf.No
## 3.98327 0.45028
## Plant.Species B :Leaf.No Plant.Species C :Leaf.No
## 0.07695 0.05255sum.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) #エラーバーの作成