Rで共分散分析

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

• 統計を学ぶ化学系技術者の記録:共分散分析.
  
• アセット･マネジメント･コンサルティング:Rで統計的仮説検定:共分散分析.
  
• 橋本洸哉:ANCOVA.
  
• こちにいるの日記.

1 サンプルデータ.

共分散分析では第１段階として、帰無仮説「各水準の回帰直線の傾きが等しい」を検定します。
アセット･マネジメント･コンサルティング様のコードをもとに作成しました。

傾きは同一で切片は異なるサンプルデータを作成する。

seed <- 20230915
set.seed(seed = seed)
n <- 20
x1 <- rnorm(n = n);x2 <- rnorm(n = n);x3 <- rnorm(n = n)
b1 <- 2;b2 <- 2;b3 <- 2
a1 <- 1;a2 <- 2;a3 <- 3
y1 <- b1 * x1 + a1 + rnorm(n)
y2 <- b2 * x2 + a2 + rnorm(n)
y3 <- b3 * x3 + a3 + rnorm(n)

1.1 データの内容.

第１水準の回帰直線について内容を確認します。
残差平方和・回帰係数・予測値を抽出し、回帰直線を表示します。

hist(x1)

fit1 <- lm(y1 ~ x1)
#str(fit1)
summary(fit1)

## 
## Call:
## lm(formula = y1 ~ x1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7091 -0.5450 -0.1864  0.3974  2.3536 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.0675     0.2226   4.795 0.000145 ***
## x1            1.7542     0.2510   6.989 1.59e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9699 on 18 degrees of freedom
## Multiple R-squared:  0.7307, Adjusted R-squared:  0.7158 
## F-statistic: 48.85 on 1 and 18 DF,  p-value: 1.586e-06

sum((fit1 $ residuals)^2)

## [1] 16.93436

fit1 $ coefficients[2]

##       x1 
## 1.754228

fit1 $ fitted.values

##             1             2             3             4             5 
## -4.239705e-01  4.370099e+00  1.083428e+00  1.072216e+00  1.361904e+00 
##             6             7             8             9            10 
##  3.696202e+00  1.335986e+00  2.607911e+00 -1.218334e+00  1.393849e+00 
##            11            12            13            14            15 
## -5.920356e-05  2.924964e+00  2.936905e-01  2.090119e+00  8.446076e-02 
##            16            17            18            19            20 
##  4.199896e-01  2.316133e-02  4.469376e+00  1.287782e+00  1.493201e+00

mean(x1)

## [1] 0.1999866

mean(y1)

## [1] 1.418299

plot(x1,y1)
abline(fit1)

#plot(fit1)

1.2 共通の傾き.

共通の傾きを求めるため、データの平均が０になるよう平行移動します。scale=F は標準偏差で割らないオプションです。

X1 <- scale(x1,center=T,scale=F)
Y1 <- scale(y1,center=T,scale=F)
X2 <- scale(x2,center=T,scale=F)
Y2 <- scale(y2,center=T,scale=F)
X3 <- scale(x3,center=T,scale=F)
Y3 <- scale(y3,center=T,scale=F)

Fit1 <- lm(Y1 ~ X1)
summary(Fit1)

## 
## Call:
## lm(formula = Y1 ~ X1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7091 -0.5450 -0.1864  0.3974  2.3536 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.681e-17  2.169e-01   0.000        1    
## X1           1.754e+00  2.510e-01   6.989 1.59e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9699 on 18 degrees of freedom
## Multiple R-squared:  0.7307, Adjusted R-squared:  0.7158 
## F-statistic: 48.85 on 1 and 18 DF,  p-value: 1.586e-06

plot(X1,Y1)
abline(Fit1)

第１水準の回帰直線ですが、原点を通るように平行移動されていることが分かります。

1.3 ３つの水準を合わせる.

３水準のデータを合併したデータフレームを作ります。水準はgroupeで、A,B,C で識別します。

df_sample <- data.frame(group = rep(c("A", "B",
    "C"), each = n), y = c(Y1, Y2, Y3),
    x = c(X1,X2,X3))
df_sample$group <- df_sample$group %>%
    factor(levels = c("A", "B", "C"))
is.data.frame(df_sample)

## [1] TRUE

glimpse(df_sample)

## Rows: 60
## Columns: 3
## $ group <fct> A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, B, B…
## $ y     <dbl> -2.28845741, 2.62205303, 0.01835167, -1.41097461, 2.29716311, 2.…
## $ x     <dbl> -1.05018826, 1.68267833, -0.19089364, -0.19728512, -0.03214797, …

1.4 ３水準の回帰直線.

３水準の回帰直線を描き、３水準を合わせた回帰直線を太めに描きます。共通の傾きの回帰直線になります。

col=c(1,2,3)
par(las=1, family="sans", cex=1.2) #図の調整
plot(y ~ x, df_sample, pch=16,
     xlab="x", ylab="y",col=c(1:3))
legend("topleft", legend=c("A", "B", "C"), pch=16, col=1:3, cex=0.8) #凡例の作成
RegressionA <- lm(y ~ x, subset(df_sample, group=="A"))
RegressionB <- lm(y ~ x, subset(df_sample,group=="B"))
RegressionC <- lm(y ~ x, subset(df_sample,group=="C"))
Regression <- lm(y ~ x,df_sample)

abline(RegressionA, col=1) #種Aの回帰線
abline(RegressionB, col=2) #種Bの回帰線
abline(RegressionC, col=3) #種Cの回帰線
abline(Regression, col="lightblue",lwd=5) #全体の回帰

共通の傾きを抽出します。

summary(Regression)

## 
## Call:
## lm(formula = y ~ x, data = df_sample)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.56347 -0.58252 -0.06709  0.63856  2.35940 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.733e-17  1.341e-01    0.00        1    
## x           1.936e+00  1.345e-01   14.39   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.039 on 58 degrees of freedom
## Multiple R-squared:  0.7812, Adjusted R-squared:  0.7774 
## F-statistic:   207 on 1 and 58 DF,  p-value: < 2.2e-16

Regression $ coefficients[2]

##        x 
## 1.935922

2 平行性検定のF統計量.

https://chemstat.hatenablog.com/entry/2020/10/04/075824

2.1 個別回帰式の残差平方和.

s2は群内変動に相当する。１つの回帰の残差平方和の自由度は18.
自由度は合計 54 となる。

s2 <- sum((RegressionA $ residuals)^2)+sum((RegressionB $ residuals)^2)+sum((RegressionC $ residuals)^2)
MSE <- s2/54
MSE

## [1] 1.145471

個別回帰式から求めた残差分散である。

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

https://chemstat.hatenablog.com/entry/2020/09/12/180346 .
s1は全変動に相当する。
各水準の重心を通り共通の傾きをもつ回帰式による残差平方和を各水準で求め、３つの残差平方和を合計する。
自由度は\(19\times 3-1=56\)である。

\(s1>s2\) であるが、３つの回帰直線の傾きが共通傾きと同じなら、\(s1=s2\) となる。
「傾きが全て等しい」という帰無仮説を棄却するには、s1−s2 が大きくなければならない。 小さければ帰無仮説は棄却されず、「３つの回帰直線は平行」は維持されることになる。
s1−s2 の自由度は 56-54=2 である。

群間変動による残差分散
\(s1=(s1-s2)+s2\) であり、\((s1-s2)\)は群間変動に相当する。
群間変動による残差分散を求める。

s1 <- sum(((Regression $ coefficients[2])*(x1-mean(x1))+mean(y1)-y1)^2)+sum(((Regression $ coefficients[2])*(x2-mean(x2))+mean(y2)-y2)^2)+sum(((Regression $ coefficients[2])*(x3-mean(x3))+mean(y3)-y3)^2)
MSA <- (s1-s2)/2
MSA

## [1] 0.3493118

F値として、s1-s2 による残差分散MSA と s2による残差分散MSEの比を
　　　f=MSA/MSE とする。

f <- MSA/MSE
f

## [1] 0.3049503

値が小さいので傾きの差は有意とは言えないようである。これについては後半で述べる。

2.3 交互作用との関係.

共分散分析はgroup A,B,C と共変量xによる２元配置分散分析である。
交互作用が有意となるのは、groupによって傾きが異なる場合である。

df_sample3 <- data.frame(group = rep(c("A", "B",
    "C"), each = n), y = c(y1, y2, y3),
    x = c(x1,x2,x3))
df_sample3$group <- df_sample3$group %>%
    factor(levels = c("A", "B", "C"))
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

ANCOVA.m <- lm(y ~ x*group, df_sample3) #ANCOVA用線形モデル
Anova(ANCOVA.m)

	Sum Sq	Df	F value	Pr(>F)
x	223.2905510	1	194.9333451	0.0000000
group	38.8797847	2	16.9710864	0.0000019
x:group	0.6986236	2	0.3049503	0.7384211
Residuals	61.8554499	54	NA	NA

交互作用のF値は項目2.2で求めたf値と一致している。
交互作用は有意でなく帰無仮説は棄却されない。
「３つの回帰係数の傾きは共通傾きと同じである」は維持される。

3 切片の差の検定.

s3は全回帰の残差平方和で、自由度は58である。
s1は共通傾きによる３群の残差平方和の合計であり、自由度は56である。
切片に差がなければ、s1,s3 は一致する。
s3-s1 による残差分散が s1 による残差分散より大きければ、帰無仮説は棄却され、切片は異なることになる。

Regression3 <- lm(y ~ x,df_sample3)
#summary(Regression3)
s3 <- sum((Regression3 $ residuals)^2)
MSA3 <- (s3-s1)/2
MSA3

## [1] 19.43989

f_value <- MSA3/(s1/56)
f_value

## [1] 17.40309

rの関数を使う場合は、交互作用が有意でなかったので、交互作用のない２元配置分散分析を行う。

library(car)
ANOVA.m <- lm(y ~ x+group, df_sample3) 
Anova(ANOVA.m)

	Sum Sq	Df	F value	Pr(>F)
x	223.29055	1	199.89539	0.0e+00
group	38.87978	2	17.40309	1.3e-06
Residuals	62.55407	56	NA	NA

groupのF値は切片の差のF値と一致している。

3.1 まとめ.

３水準を合わせた全変動を \(S_T\) とすると、
　 \(S_T=(S_T-s3)+s3\)　　　（回帰変動＋残差変動）
残差変動をさらに分解し、
　\(S_T=(S_T-s3)+(s3-s1)+s1\)

　\(S_T=(S_T-s3)+(s3-s1)+(s1-s2)+s2\)

• \((s1-s2)\)により傾きが等しいかどうか検定し、有意でなければ.
  
• \((s3-s1)\)により切片が等しいかどうか検定する。

一般に２元配置分散分析では、

• はじめに交互作用が有意かどうか調べ、
• 交互作用が有意でなければ、交互作用のない検定をあらためて行う

傾きの検定は、交互作用の検定に現れる。

橋本洸哉先生のサイト参照.
ANCOVA.

２元配置分散分析（交互作用）で、傾きの検定と同じ内容が出力される。

• 傾き共通を調べる検定量　f <- \(MSA/MSE \sim F(2,54)\)　

• 切片の差を調べる統計量　f_value<- \(MSA3/(s1/56) \sim F(2,56)\)

これらは本当にF分布に従っているであろうか？
体感的に理解するには、シミュレーションが有効である。

4 傾きのシミュレーション.

帰無仮説のもとで傾きの検定量を100回出力し、F分布になるかどうか確認しよう。

X=c()
for (i in 1:100){
  n <- 20
  x1 <- rnorm(n = n);x2 <- rnorm(n = n);x3 <- rnorm(n = n)
  b1 <- 2;b2 <- 2;b3 <- 2
  a1 <- 1;a2 <- 2;a3 <- 3
  y1 <- b1 * x1 + a1 + rnorm(n)
  y2 <- b2 * x2 + a2 + rnorm(n)
  y3 <- b3 * x3 + a3 + rnorm(n)
  X1 <- scale(x1,center=T,scale=F)
  Y1 <- scale(y1,center=T,scale=F)
  X2 <- scale(x2,center=T,scale=F)
  Y2 <- scale(y2,center=T,scale=F)
  X3 <- scale(x3,center=T,scale=F)
  Y3 <- scale(y3,center=T,scale=F)
  df_sample <- data.frame(group = rep(c("A", "B",
    "C"), each = n), y = c(Y1, Y2, Y3),
    x = c(X1,X2,X3))
  df_sample$group <- df_sample$group %>%
    factor(levels = c("A", "B", "C"))
  RegressionA <- lm(y ~ x, subset(df_sample, group=="A"))
  RegressionB <- lm(y ~ x, subset(df_sample,group=="B"))
  RegressionC <- lm(y ~ x, subset(df_sample,group=="C"))
  Regression <- lm(y ~ x,df_sample)
  s2 <- sum((RegressionA $ residuals)^2)+sum((RegressionB $ residuals)^2)+sum((RegressionC $ residuals)^2)
  MSE <- s2/54
  s1 <- sum(((Regression $ coefficients[2])*(x1-mean(x1))+mean(y1)-y1)^2)+sum(((Regression $ coefficients[2])*(x2-mean(x2))+mean(y2)-y2)^2)+sum(((Regression $ coefficients[2])*(x3-mean(x3))+mean(y3)-y3)^2)
  MSA <- (s1-s2)/2
  f <- MSA/MSE
  X <- c(X,f)
}
hist(X,prob=T)
curve(df(x,2,54),add=T)

理論上のF分布のグラフと一致していることが確認できる。
F=0.3049 は有意でないことも確認できる。

5 市原.

傾きが等しいかどうか検定する。

x1 <- c(54,50,39,44,35,41,32,46,38,30)
x2 <- c(43,39,49,42,35,31,36,34,27)
y1 <- c(17,15,15,14,13,12,11,10,10,9)
y2 <- c(22,21,19,18,16,15,13,11,9)
X1 <- scale(x1,center=T,scale=F)
Y1 <- scale(y1,center=T,scale=F)
X2 <- scale(x2,center=T,scale=F)
Y2 <- scale(y2,center=T,scale=F)
df_sample <- data.frame(group = c(rep("A",10),rep("B",9)), y = c(Y1, Y2),x = c(X1,X2))
df_sample$group <- df_sample$group %>%
    factor(levels = c("A", "B"))
is.data.frame(df_sample)

## [1] TRUE

glimpse(df_sample)

## Rows: 19
## Columns: 3
## $ group <fct> A, A, A, A, A, A, A, A, A, A, B, B, B, B, B, B, B, B, B
## $ y     <dbl> 4.4, 2.4, 2.4, 1.4, 0.4, -0.6, -1.6, -2.6, -2.6, -3.6, 6.0, 5.0,…
## $ x     <dbl> 13.100000, 9.100000, -1.900000, 3.100000, -5.900000, 0.100000, -…

col=c(1,2)
par(las=1, family="sans", cex=1.2) #図の調整
plot(y ~ x, df_sample, pch=16,
     xlab="x", ylab="y",col=c(1:2))
legend("topleft", legend=c("A", "B"), pch=16, col=1:2, cex=0.8) #凡例の作成
RegressionA <- lm(y ~ x, subset(df_sample, group=="A"))
RegressionB <- lm(y ~ x, subset(df_sample,group=="B"))
Regression <- lm(y ~ x,df_sample)

abline(RegressionA, col=1) #種Aの回帰線
abline(RegressionB, col=2) #種Bの回帰線
abline(Regression, col="lightblue",lwd=4) #全体の回帰

Regression $ coefficients[2]

##         x 
## 0.3478553

summary(RegressionA)

## 
## Call:
## lm(formula = y ~ x, data = subset(df_sample, group == "A"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7880 -0.9515  0.3767  1.1808  2.8426 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 6.658e-16  6.459e-01   0.000   1.0000  
## x           2.329e-01  8.831e-02   2.638   0.0298 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.043 on 8 degrees of freedom
## Multiple R-squared:  0.4651, Adjusted R-squared:  0.3983 
## F-statistic: 6.957 on 1 and 8 DF,  p-value: 0.02982

summary(RegressionB)

## 
## Call:
## lm(formula = y ~ x, data = subset(df_sample, group == "B"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2682 -2.3073 -0.4246  2.2905  4.1341 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 1.174e-15  9.869e-01    0.00   1.0000  
## x           5.196e-01  1.565e-01    3.32   0.0128 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.961 on 7 degrees of freedom
## Multiple R-squared:  0.6116, Adjusted R-squared:  0.5561 
## F-statistic: 11.02 on 1 and 7 DF,  p-value: 0.01276

df_sample$y

##  [1]  4.4  2.4  2.4  1.4  0.4 -0.6 -1.6 -2.6 -2.6 -3.6  6.0  5.0  3.0  2.0  0.0
## [16] -1.0 -3.0 -5.0 -7.0

df_sample$group

##  [1] A A A A A A A A A A B B B B B B B B B
## Levels: A B

boxplot(y ~ group, df_sample)

t.test(y1, y2, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  y1 and y2
## t = -2.0551, df = 17, p-value = 0.05556
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -6.89044963  0.09044963
## sample estimates:
## mean of x mean of y 
##      12.6      16.0

data <- data.frame(group = c(rep("A",10),rep("B",9)), y = c(y1, y2),x = c(x1,x2))
ANOVA.m <- lm(y ~ group, data) 
summary(ANOVA.m)

## 
## Call:
## lm(formula = y ~ group, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -7.0   -2.6    0.0    2.4    6.0 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   12.600      1.139  11.066 3.44e-09 ***
## groupB         3.400      1.654   2.055   0.0556 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.601 on 17 degrees of freedom
## Multiple R-squared:  0.199,  Adjusted R-squared:  0.1519 
## F-statistic: 4.224 on 1 and 17 DF,  p-value: 0.05556

library(car)
ANCOVA.m <- lm(y ~ x*group, df_sample) #ANCOVA用線形モデル
Anova(ANCOVA.m)

	Sum Sq	Df	F value	Pr(>F)
x	108.04386	1	17.106606	0.0008797
group	0.00000	1	0.000000	1.0000000
x:group	17.61743	1	2.789371	0.1156206
Residuals	94.73871	15	NA	NA