共分散分析

共分散分析（ANCOVA）を学習した記録です。
参考にさせていただいたサイト様.

• 統計を学ぶ化学系技術者の記録:共分散分析.
  模擬データも、こちらのサイト様のものを使わせていただきました。

回帰直線の差を共分散分析（ANCOVA）を使って検定する。

• 傾きが等しいか（平行性の検定）.
• パラメーターに意味があるか（回帰の有意性の検定）.
• 切片が等しいか（群間差の検定）.

https://y2pt.com/8214/

1 共変量xで３群比較.

1.1 ３群のデータ.

data1 <- data.frame(x=c(1:5),y=c(2,6,6,9,6),z=rep(1,5))
data2 <- data.frame(x=c(2:6),y=c(0,4,4,3,6),z=rep(2,5))
data3 <- data.frame(x=c(0:4),y=c(3,5,7,9,7),z=rep(3,5))
data <- rbind(data1,data2,data3)
apply(data,2,mean)

##        x        y        z 
## 3.000000 5.133333 2.000000

mx <- mean(data$x)
my <- mean(data$y)
x1 <- data1$x
y1 <- data1$y
x2 <- data2$x
y2 <- data2$y
x3 <- data3$x
y3 <- data3$y

#par(mfrow=c(1,2))
plot(y ~ x, data, pch=16,
     xlab="x", ylab="y",col=c(1:3),ylim=c(-1,10))
legend("topleft", legend=c("1", "2", "3"), pch=16, col=1:3, cex=0.8) #凡例の作成

1.2 傾き共通モデル.

３群をまとめて単回帰分析をすると、

　　　全変動＝群内変動＋群間変動.

により群間変動の影響を受け共通の回帰係数が正しく計算されない.
そこで、
各群の重心を(0,0)に移動し、群間変動をゼロにしたうえで分散・共分散を計算すると、それは各群の分散・共分散の合計になる.
回帰係数は、
　　　\(\begin{align*} \hat{b} = \frac{S_{xy1}+S_{xy2}+S_{xy3}}{S_{xx1}+S_{xx2}+S_{xx3}}=\frac{11+11+12}{10+10+10}=1.133... \end{align*}\).

によって算出される.
（分母は群内変動、自由度で割ればプールされた分散であり、分子はプールされた共分散である。２次元混合分布の平方和分解とも言える）.

詳しくは、
　共分散分析（ANCOVA）：平行性の検定①傾き共通モデルを作る.
に説明があります.

1.3 pooled covariance.

#各群の偏差積和
cov(data1$x,data1$y)*4

## [1] 11

cov(data2$x,data2$y)*4

## [1] 11

cov(data3$x,data3$y)*4

## [1] 12

# pooled covariance（自由度12）.  
pool_cov <- cov(data1$x,data1$y)*4+cov(data2$x,data2$y)*4+cov(data3$x,data3$y)*4;pool_cov

## [1] 34

#プールしない共分散（自由度12）
covsum <- cov(data$x,data$y)*14;covsum

## [1] 20

# pooled variance（自由度12）
pool_varx <- var(data1$x)*4+var(data2$x)*4+var(data3$x)*4;pool_varx

## [1] 30

#プールしない分散（自由度14）
varx <- var(data$x)*14;varx

## [1] 40

• プール分散 ＜ 全分散.
  
• プール共分散 ＞ 全共分散.
  （共分散はマイナスがありうるためと考えられる）.

1.4 共通の傾き.

プールされた共分散・分散から算出する.

#回帰係数（共通傾き）. 
b_poor <- pool_cov/pool_varx;b_poor

## [1] 1.133333

#全データの傾き
b <- covsum/varx;b

## [1] 0.5

群間変動の影響で、両者の値は異なっている.

1.5 回帰.

plot(y ~ x, data, pch=16,
     xlab="x", ylab="y",col=c(1:3),ylim=c(-1,10),asp=1)
Regression0 <- lm(y ~ x,data)
Regression1 <- lm(y ~ x, subset(data, z==1))
Regression2 <- lm(y ~ x, subset(data,z==2))
Regression3 <- lm(y ~ x, subset(data,z==3))
Regression0$ coefficients

## (Intercept)           x 
##    3.633333    0.500000

Regression1$ coefficients

## (Intercept)           x 
##         2.5         1.1

Regression2$ coefficients

## (Intercept)           x 
##        -1.0         1.1

Regression3$ coefficients

## (Intercept)           x 
##         3.8         1.2

abline(Regression0, col="grey",lwd=1) #全体の回帰
abline(Regression1, col=1) #種1の回帰線
abline(Regression2, col=2) #種2の回帰線
abline(Regression3, col=3) #種3の回帰線
abline(my-mx*pool_cov/pool_varx,pool_cov/pool_varx, col="lightblue",lwd=3) 

b1 <- Regression1$ coefficients[2]
b2 <- Regression2$ coefficients[2]
b3 <- Regression3$ coefficients[2]

プールされた分散・共分散の方がよく当てはまっている。
全データによる回帰係数の傾きは、この場合は意味をなさないことが分かる。

2 平行性の検定.

平行性の検定では、帰無仮説「すべての回帰直線の傾きが等しい」ことを検定数するため、

• 傾き共通の回帰式の残差平方和S1の残差.
• 個別の回帰式の残差平方和S2の残差.
  を比較する.

2.1 残差.

s1は残差の全変動に相当する。自由度は \(4\times 3 -1=11\)。
s2は残差の群内変動に相当する。１つの回帰の残差平方和の自由度は3. 自由度は合計 9 となる。

#残差の全変動
b_poor

## [1] 1.133333

s1 <- sum((b_poor*(x1-mean(x1))+mean(y1)-y1)^2)+sum((b_poor*(x2-mean(x2))+mean(y2)-y2)^2)+sum((b_poor*(x3-mean(x3))+mean(y3)-y3)^2);s1

## [1] 26.26667

s1/11

## [1] 2.387879

#残差の群内変動
s2 <- sum((Regression1$ residuals)^2)+sum((Regression2$ residuals)^2)+sum((Regression3$ residuals)^2);s2

## [1] 26.2

s2/9

## [1] 2.911111

#共通の傾きによる回帰変動
s3 <- ((b_poor-b1)^2*var(x1)+(b_poor-b2)^2*var(x2)+(b_poor-b3)^2*var(x3))*4;s3

##          x 
## 0.06666667

s3/2

##          x 
## 0.03333333

s1

## [1] 26.26667

s2+s3

##        x 
## 26.26667

変動＝群内変動＋群間変動　に相当する平方和分解は、

　　\(s1=s2+s3\)　　　であることが分かる.
数式での確認は次のようになると思う。

傾き共通の回帰式の残差平方和S1.

\(b=bpoor\)とおいて、

\(\displaystyle s1=\sum_i\left[b(x_{1i}-\bar{x_1})+\bar{y_1}-y_{1i}\right]^2+\sum_j\left[b(x_{2i}-\bar{x_2})+\bar{y_2}-y_{2i}\right]^2+\sum_k\left[b(x_{3i}-\bar{x_3})+\bar{y_3}-y_{3i}\right]^2\).
第１群の傾きを\(b_1\)とすると第一項は.

\(\displaystyle \sum_i\left[(b-b_1+b_1)(x_{1i}-\bar{x_1})+\bar{y_1}-y_{1i}\right]^2=\sum_i(b-b_1)^2(x_{1i}-\bar{x_1})^2+\sum_i\left[b_1(x_{1i}-\bar{x_1})+\bar{y_1}-y_{1i}\right]^2\).

第２項、第３項についても同様に平方和を分解し、辺々足すと.

\(\displaystyle s1=\sum_i(b-b_1)^2(x_{1i}-\bar{x_1})^2+\sum_i(b-b_2)^2(x_{2i}-\bar{x_2})^2+\sum_i(b-b_3)^2(x_{3i}-\bar{x_3})^2+s2\).

群間変動 \(s3\)は.

\(\displaystyle s3=(b-b_1)^2\sum_i(x_{1i}-\bar{x_1})^2+(b-b_2)^2\sum_i(x_{2i}-\bar{x_2})^2+(b-b_3)^2\sum_i(x_{3i}-\bar{x_3})^2\).

と考えられる。

一般に単回帰係数の分散は.

　　　\(\sigma ^{2}\left( X'X\right) ^{-1}= \frac{\sigma^2}{\sum_{i=1}^{n}( x_{i}-\overline{x})^2} \begin{pmatrix} \overline {x^2} & -\overline x \\ -\overline x & 1 \end{pmatrix}\) であり、
\(\displaystyle V(b)= \frac{\sigma^2}{\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^{2}}\)　から.
　　\(\displaystyle V(b)\cdot {\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^2}= \sigma^2\)　　即ち　\(\displaystyle E\left[(b_1-b)^2\right]\cdot {\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^2}= \sigma^2\)

よって、
\(s3\)の各項の期待値は 残差分散 \(\sigma^2\) であり、
共通の傾き b を推定値として使っているので、自由度は１減じて２と考えられる。.

整理すると、

＜帰無仮説＞　\(b=b_1=b_2=b_3\) のもとで、

• \(E\left(\dfrac{s2}{9}\right)=\sigma^2\)　　　\(E\left(\dfrac{s3}{2}\right)=\sigma^2\).
  
• \(F=\dfrac{s3/2}{s2/9}\) は自由度\((2,9)\)のF分布に従う。

2.2 F検定.

MSA <- (s1-s2)/2
MSA

## [1] 0.03333333

#F値
 MSA/(s2/9)

## [1] 0.01145038

• F値は 0.0114504.
  
• P値は 0.9886293.

傾きが共通であるという帰無仮説は棄却されない。
これは、anova では交互作用が有意でないということに相当する.

そこで、単純主効果である切片の検定にうつる.

2.3 anova による出力.

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'.

#ANCOVA.m <- lm(Egg.No ~ Plant.Species*Leaf.No, OvipPrefData1) #ANCOVA用線形モデル
#Anova(ANCOVA.m) #ANCOVAの実践
#(lsm.ANCOVA.m <- lsmeans(ANCOVA.m, specs="Plant.Species"))

#data <- rbind(data1,data2,data3)
data

x	y	z
1	2	1
2	6	1
3	6	1
4	9	1
5	6	1
2	0	2
3	4	2
4	4	2
5	3	2
6	6	2
0	3	3
1	5	3
2	7	3
3	9	3
4	7	3

library(car)

## Loading required package: carData

x <- data$x
y <- data$y
z <- data$z
res <- aov(y ~ x + z)
interaction.plot(x,z,y, type="b",  pch=c(1,2))

summary(res)

##             Df Sum Sq Mean Sq F value Pr(>F)
## x            1  10.00  10.000   1.588  0.232
## z            1   2.16   2.160   0.343  0.569
## Residuals   12  75.57   6.298

#se <- function(x) sd(x)/sqrt(length(x))

simpleANOVA.m <- lm(y ~ z, data) #one-way ANOVA用線形モデルの作成
Anova(simpleANOVA.m) #one-way ANOVA

	Sum Sq	Df	F value	Pr(>F)
z	0.40000	1	0.059542	0.81103
Residuals	87.33333	13	NA	NA

library(ggplot2)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some()   masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

attach(data)

## The following objects are masked _by_ .GlobalEnv:
## 
##     x, y, z

par(las=1, family="sans", cex=1.2) #図の調整

plot(y ~ x, data, 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(data)

x	y	z
1	2	1
2	6	1
3	6	1
4	9	1
5	6	1
2	0	2

simpleRegressionA <- lm(y ~ x, subset(data, z=="1")) #種Aのデータで単回帰
simpleRegressionB <- lm(y ~ x, subset(data, z=="2")) #種Bのデータで単回帰
simpleRegressionC <- lm(y ~ x, subset(data, z=="3")) #種Cのデータで単回帰
plot(y ~ x, data, 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) #凡例の作成

abline(simpleRegressionA, col=1) #種Aの回帰線
abline(simpleRegressionB, col=2) #種Bの回帰線
abline(simpleRegressionC, col=3) #種Cの回帰線

ANCOVA.m <- lm(y ~ factor(z)*x, data) #ANCOVA用線形モデル
Anova(ANCOVA.m) #ANCOVAの実践

	Sum Sq	Df	F value	Pr(>F)
factor(z)	51.4666667	2	8.8396947	0.0075214
x	38.5333333	1	13.2366412	0.0054159
factor(z):x	0.0666667	2	0.0114504	0.9886293
Residuals	26.2000000	9	NA	NA

傾きの有意性は　factor(z):x　に現れる（交互作用）

重回帰分析の回帰平面が歪んでいるかの検定と同等である.
回帰平面がフラットならば、回帰直線は平行になる

3 切片差の検定.

s3 <- sum((Regression0$ residuals)^2);s3

## [1] 77.73333

s3/13

## [1] 5.979487

#残差の群内変動
s1 <- sum((b_poor*(x1-mean(x1))+mean(y1)-y1)^2)+sum((b_poor*(x2-mean(x2))+mean(y2)-y2)^2)+sum((b_poor*(x3-mean(x3))+mean(y3)-y3)^2);s1

## [1] 26.26667

s1/11

## [1] 2.387879

s3-s1

## [1] 51.46667

(s3-s1)/2

## [1] 25.73333

F <- ((s3-s1)/2)/(s2/11)
F

## [1] 10.80407

群z は因子(factor)としてanovaを適用する。

• 因子型だと、自由度は 3-1=2 であるが、
  
• 数値型だと、自由度は 1 として算出される.

x <- data$x
y <- data$y
z <- data$z
res <- aov(y ~ x + factor(z))
summary(res)

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## x            1  10.00  10.000   4.188 0.06538 . 
## factor(z)    2  51.47  25.733  10.777 0.00256 **
## Residuals   11  26.27   2.388                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rの関数 anova を使う場合は、交互作用が有意でなかったので、交互作用のない２元配置分散分析を行う。
factor(z)のF値は切片の差のF値と一致している。