【統計】共分散分析（ANCOVA） - こちにぃるの日記

共分散分析についてめも。

【目次】

帰無仮説、対立仮設
計算式等
計算例
プログラムコード
交互作用項での平行の評価
参考

共分散分析は、分散分析と重回帰分析を組み合わせた手法と言われる。
分散分析では説明変数に集団などカテゴリー変数（要因）のみを投入するが、共分散分析ではカテゴリー変数（要因）と連続変数（共変量）の両方を投入するのが特徴。

分散分析を ANOVA (ANalysis Of VAriance) と言うのに対し、共分散分析は ANCOVA (ANalysis of COVAriance) と呼称される。

帰無仮説、対立仮設

帰無仮説 $H_{0}:$ 要因の水準間の主効果が同じ（平均に差がない）
対立仮設 $H_{1}:$ 要因の水準間の主効果に差がある（平均に差がある）

計算式等

ANOVAでは、水準間の変動が残差の変動（ランダムな変動）より有意に大きいかどうかで、水準間の差を評価する。
ANCOVAでは、目的変数と相関（線形関係）を持つ共変量を投入し、その影響を調整（共通の傾きを持つ回帰直線を作成）した上で、水準間の変動が残差の変動より有意に大きいかを評価する。

以下、回帰式を記載。重回帰分析の記事と同じ。
※カテゴリー変数１つ（２水準）、連続変数（共変量）１つの例。

■ 水準毎の回帰直線（水準毎の条件付き期待値）
$\quad \widehat{ys}_{ki} = \beta_{0_{k}} + \beta_{1_{k}}c_{ki}$
　※ k：要因の水準、c：共変量
　※ 例えば、要因の水準がA薬・B薬の２つなら薬剤毎に計２つの回帰直線を作成。

■ 共通の回帰直線（共通の条件付き期待値）
$\quad \widehat{yc}_{i} = \gamma_{0} + \gamma_{1}c_{i} + \gamma_{2} x_{i}$
　※ c：共変量、x：要因
　※ 共変量 c の傾きを一定として、要因 x による違いを表現した直線。
　※ 例えば、要因の水準がA薬・B薬の２つなら、x=0、x=1の2値で $\gamma_{2}$ 分が期待値の差になる。

■ 共通の回帰直線の重心（総平均）
$\quad \delta_{0} = \sum \{ \left( \gamma_{0} + \gamma_{2} x_{k} \right) N_{k} \} / N$
$\quad G = \delta_{0} + \gamma_{1} \overline{c}$
　※ $\gamma_{2} x_{k}$ ：水準kの回帰係数（xは０or１のダミー変数）、 $N_{k}$ =水準kのN数

■ 分散分析表

変動要因		平方和	自由度	平均平方	F値
B ：共通傾きでの回帰直線(水準毎) vs. 共通の回帰直線の重心		$S_{B}=\sum\left( \widehat{yc}_{i} - G \right)^{2}$	$df_{B}=a-1$	$V_{B}=S_{B}/df_{B}$	$F_{B}=V_{B}/V_{W}$
W：観測値 vs. 共通傾きでの回帰直線(水準毎)		$S_{W}=\sum\left( y_{i} - \widehat{yc}_{i} \right)^{2}$	$df_{W}=df_{T}-df_{B}$	$V_{W}=S_{W}/df_{W}$	$F_{B}=V_{B}/V_{W}$
	R：回帰直線(水準毎) vs. 共通傾きでの回帰直線(水準毎)	$S_{R}=\sum\sum\left( \widehat{ys}_{ik} - \widehat{yc}_{i} \right)^{2}$	$df_{R}=a-1$	$V_{R}=S_{R}/df_{R}$	$F_{R}=V_{R}/V_{e}$
	e：観測値 vs. 回帰直線(水準毎)	$S_{e}=\sum\sum\left( y_{i} - \widehat{ys}_{ik} \right)^{2}$	$df_{e}=df_{T}-df_{R}$	$V_{e}=S_{e}/df_{e}$	$F_{R}=V_{R}/V_{e}$
T ：観測値 vs. 共通の回帰直線の重心		$S_{T}=\sum\left( y_{i} - G \right)^{2}$	$df_{T}=n-1$	$V_{T}=S_{T}/df_{T}$	---

　※k=水準k、a=要因の水準数＋共変量の数。

（各平方和の図解）

f:id:cochineal19:20210222182649p:plain

■ 調整平均（LS mean：Least Square mean）

水準毎の調整平均　※要因 x を水準AとBで考えた場合。
$\quad \overline{yc}_{A} = \gamma_{0} + \gamma_{1}\overline{c} + \gamma_{2} x_{A}$
$\quad \overline{yc}_{B} = \gamma_{0} + \gamma_{1}\overline{c} + \gamma_{2} x_{B}$
調整平均の差とその信頼区間
$\quad d_{y_{\verb|adj|}} = \gamma_{2}$
$\quad CL=d_{y_{\verb|adj|}} \pm \verb|T.INV|\left(α,\ n-3\right) \cdot se[\gamma_{2}]$
　※ $se[\gamma_{2}]$ ：共通の回帰直線における、共変量の回帰係数 $\gamma_{2}$ の標準誤差（Standard Error）

（補足）通常の平均の差とその信頼区間
$\quad d_{y}=\overline{y_{B}} - \overline{y_{A}}$
$\quad CL=d_{y} \pm \verb|T.INV|\left(α,\ n-2\right) \cdot \sqrt{\dfrac{\sigma ^{2}}{n_{A}}+\dfrac{\sigma ^{2}}{n_{B}}}$

計算例

今回も架空データで計算手順を確認。
　BASE = {21, 15, 18, 16, 26, 25, 22, 21, 16, 17, 18}
　TRT01AN = {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1}
　CHG = {-7, -2, -5, -4, -12, -15, -12, -12, -4, -7, -7}
　AVAL = {14, 13, 13, 12, 14, 10, 10, 9, 10, 10, 11}

特定時点における検査Xのベースラインからの変化量 (CHG) を目的変数として、A薬（TRT01AN=0）とB薬（TRT01AN=1）による効果の違いを評価するシナリオとする。この時、共変量として検査Xのベースライン値 (BASE) を投入。

■ 分析モデル：
　ベースラインからの変化量＝投与群（要因）＋ベースライン値（共変量）

回帰の同質性の検定

先ず、水準毎の回帰直線（CHG=BASE）を描画する。
図のとおり、若干傾きが異なり、A薬とB薬の差をどこで見るかで結果が変わってしまうことが分かる。
$\quad \hat{ys_{Ai}}=10.325 - 0.850 c_{Ai}$
$\quad \hat{ys_{Bi}}=11.011 - 1.051 c_{Bi}$

f:id:cochineal19:20210222231956p:plain

ここに、共通の傾きの回帰直線（CHG=BASE+TRT01AN）を重ねる。
この赤・青の実線であれば、共変量の値によらずA薬とB薬の差を一定に保つことができる。
$\quad \hat{yc_{i}} = 11.972 - 0.936 c_{i} - 3.240 x_{i}$

f:id:cochineal19:20210223115640p:plain

共分散分析では、この共通の傾きの回帰直線（赤・青の実線）を用いるが、そのためには、水準毎の回帰直線（赤・青の点線）が平行であるという前提が必要。
この前提の検証として、下図のとおり、 $S_{W}$ （共通の傾きの回帰直線からの残差平方和）を $S_{R}$ （水準毎の回帰直線からの残差平方和）と $S_{e}$ （２本の回帰直線の偏差平方和）に分解して F検定を行うことで、水準毎の回帰直線が平行と解釈できるかを評価する（回帰の同質性の検定、帰無仮説：水準毎の回帰直線が平行）。

f:id:cochineal19:20210222233701p:plain

先ほどの回帰式を個々の値に当てはめて、期待値と偏差平方を出したものが下表。
※Y共通傾き列：共通の傾きを用いた各水準の回帰直線での期待値、Y単体列：各水準単体の回帰直線での期待値、R列：Y共通傾きとY単体の偏差平方、e列：YCHGとY単体の偏差平方
＜再掲＞
$\quad \hat{ys_{Ai}}=10.325 - 0.850 c_{Ai}$
$\quad \hat{ys_{Bi}}=11.011 - 1.051 c_{Bi}$
$\quad \hat{yc_{i}} = 11.972 - 0.936 c_{i} - 3.240 x_{i}$

C_BASE	X_TRT01AN	Y_CHG	Y_共通傾き	Y_単体	R	e
21	0	-7	-7.68	-7.53	0.02	0.28
15	0	-2	-2.07	-2.43	0.13	0.18
18	0	-5	-4.88	-4.98	0.01	0.00
16	0	-4	-3.00	-3.28	0.08	0.52
26	0	-12	-12.37	-11.78	0.34	0.05
25	1	-15	-14.67	-15.26	0.35	0.07
22	1	-12	-11.86	-12.11	0.06	0.01
21	1	-12	-10.93	-11.06	0.02	0.88
16	1	-6	-6.25	-5.80	0.19	0.04
17	1	-7	-7.18	-6.86	0.11	0.02
18	1	-7	-8.12	-7.91	0.04	0.82

R列、e列をそれぞれ足し合わせ平方和を算出し、F値、p値を求める。

変動要因	平方和	自由度	平均平方	F値	p値
R：回帰直線(水準毎) vs. 共通傾きでの回帰直線(水準毎)	1.357	2	0.679	1.4139	0.3140
e：観測値 vs. 回帰直線(水準毎)	2.880	6	0.480	1.4139	0.3140

p > 0.05 で非有意であれば、水準毎の回帰直線は平行であると解釈して、以降、共通の傾きでの回帰直線を用いて共分散分析を行う。
今回の架空データでは p=0.3140で非有意のため、A薬・B薬の回帰直線は平行と解釈し、共分散分析に進む。

（※ 水準毎の回帰直線が平行であることの評価方法として、交互作用項を含めたモデルを作り、交互作用項が非有意なら平行と解釈する方法もある。）

共分散分析

先ず、共通の回帰直線における重心（総平均）を考える。
※今回、A薬はN=5, B薬はN=6の全体N=11。A薬を x=0、B薬を x=1 。
$\quad \delta_{0} = \sum \{ \left( \gamma_{0} + \gamma_{2} x_{k} \right) N_{k} \} / N=\{ 11.972 \times 5 + \left(11.972 - 3.240 \right) \times 6 \} / 11 ≒ 10.205$
$\quad G = \delta_{0} + \gamma_{1} \overline{c} = 10.205 - 0.936 \times 19.546 ≒ -8.091$

重心が算出できたら同質性の検定時と同じ要領で偏差平方を求める。
※T列：YCHGと重心との偏差平方、B列：Y単体と重心との偏差平方、W列：YCHGとY共通傾きの偏差平方

f:id:cochineal19:20210223002545p:plain

C_BASE	X_TRT	Y_CHG	AVAL	Y_共通傾き	Y_単体	T	B	W
21	0	-7	14	-7.68	-7.53	1.19	0.16	0.47
15	0	-2	13	-2.07	-2.43	37.10	36.27	0.00
18	0	-5	13	-4.88	-4.98	9.55	10.33	0.02
16	0	-4	12	-3.00	-3.28	16.74	25.87	0.99
26	0	-12	14	-12.37	-11.78	15.28	18.27	0.13
25	1	-15	10	-14.67	-15.26	47.74	43.28	0.11
22	1	-12	10	-11.86	-12.11	15.28	14.22	0.02
21	1	-12	9	-10.93	-11.06	15.28	8.03	1.15
16	1	-6	10	-6.25	-5.80	4.37	3.41	0.06
17	1	-7	10	-7.18	-6.86	1.19	0.83	0.03
18	1	-7	11	-8.12	-7.91	1.19	0.00	1.25

T列、B列、W列をそれぞれ足し合わせ平方和を算出し、F値、p値を求める。

変動要因		平方和	自由度	平均平方	F値	p値
B ：共通傾きでの回帰直線(水準毎) vs. 共通の回帰直線の重心		160.672	2	80.336	151.6721	0.0000
W：観測値 vs. 共通傾きでの回帰直線(水準毎)		4.237	8	0.530	151.6721	0.0000
	R：回帰直線(水準毎) vs. 共通傾きでの回帰直線(水準毎)	1.357	2	0.679	1.4139	0.3140
	e：観測値 vs. 回帰直線(水準毎)	2.880	6	0.480	1.4139	0.3140
T ：観測値 vs. 共通の回帰直線の重心		164.909	10	16.491	---	---

　※薄黄色は先ほどの同質性の検定の部分。

この表の $S_{B}$ （水準間の平方和）と $S_{W}$ （共通の傾きの回帰直線からの残差平方和）の平均平方を比較することで、水準間の変動がランダムな変動より有意に大きいかを評価する。今回の架空データでは p < 0.001 で水準間に有意な変動がある。

（追記）
SASの Output の Type II または III を見ると F (1,1)=53.64, p<0.0001 で薬剤（TRT01AN）の主効果が有意だったことが分かる。Type X 平方和は、共分散分析モデルの要因・共変量（TRT01AN、BASE）を分解して、要因別の主効果の有無を評価したもの。

f:id:cochineal19:20210223215534p:plain

※ Type II, III 平方和の計算は省略。平方和の違いはいつかまとめたい。
※ Type I 平方和のTRT01ANは次のとおり。要否別で備忘録として。
$\quad \left( \overline{y_{\verb|A薬|}} - G \right)^{2} \cdot N_{\verb|A薬|} + \left( \overline{y_{\verb|B薬|}} - G \right)^{2} \cdot N_{\verb|B薬|} = \left( -6.000 - \left(-8.091 \right) \right)^{2} \cdot 5 + \left( -9.833 - \left(-8.091 \right) \right)^{2} \cdot 6 ≒ 40.0758$

調整平均（LS mean：Least Square mean）

共分散分析と一緒に調整平均の差とその信頼区間を示すこともあるため備忘録。
今回の架空データを ExcelのLINEST関数で実行した結果：

f:id:cochineal19:20210223104601p:plain
また、共変量（BASE）の平均は19.545だったため、調整平均は以下となる。
水準毎の調整平均
$\quad \overline{yc} = \gamma_{0} + \gamma_{1}\overline{c} + \gamma_{2} x=11.972-0.936\times19.545 + \begin{pmatrix} -3.240\times 0 \\ -3.240\times 1 \end{pmatrix}=\begin{pmatrix} -6.323 \\ -9.564\end{pmatrix}$
調整平均の差とその信頼区間
$\quad d_{y_{\verb|adj|}} = \gamma_{2}=-3.240$
$\quad CL=d_{y_{\verb|adj|}} \pm \verb|T.INV|\left(0.05, 8\right) \cdot se[\gamma_{2}]=-3.240 \pm 2.3060 \times 0.4424=[ -4.2608,\ -2.2202]$

これを通常の平均と比べると下表のとおり。

評価項目	A薬	B薬	差 (B-A)	95%信頼区間
Y_CHGの平均	-6.000	-9.833	-3.833	-8.9349	1.2682
Y_CHGの調整平均（LS mean）	-6.323	-9.564	-3.240	-4.2608	-2.2202

今回の架空データでは、通常の平均の差の信頼区間は０を挟むのに対し、調整平均では信頼区間の幅が狭まり、０を挟まなくなった（信頼区間下限でもB薬の方が効果を示している）。

Rでの実行：

library(tidyverse)
library(car)

#-- サンプルデータ
ADS <- data.frame(
  TRT01AN=c(0,0,0,0,0,1,1,1,1,1,1)
  ,BASE=c(21,15,18,16,26,25,22,21,16,17,18)
  ,AVAL=c(14,13,13,12,14,10,10,9,10,10,11))
ADS$CHG <- ADS$AVAL - ADS$BASE
ADS$TRT01AF <- relevel(factor(ifelse(ADS$TRT01AN==0,"A薬","B薬")), ref="A薬")

#-- 水準毎の回帰分析
ADS.A <- ADS[ADS$TRT01AF=="A薬",]; ANOVA1.A <- lm(CHG ~ BASE, data=ADS.A); summary(ANOVA1.A);
ADS.B <- ADS[ADS$TRT01AF=="B薬",]; ANOVA1.B <- lm(CHG ~ BASE, data=ADS.B); summary(ANOVA1.B)
g1 <- ggplot() +
  geom_point(data=ADS, mapping=aes(x=BASE, y=CHG, color=TRT01AF), size=5) +
  geom_smooth(data=ADS.A, mapping=aes(x=BASE, y=CHG), method="lm", se=F, color="red", lty=2) +
  geom_smooth(data=ADS.B, mapping=aes(x=BASE, y=CHG), method="lm", se=F, color="blue", lty=2) +
  labs(x="ベースラインの検査値（共変量）", y="ベースラインからの変化量"
       ,title="検査Xのベースラインからの変化量", color="投与群",caption="🄫Cochineal19") +
  theme(plot.title=element_text(size=30)
        ,axis.title=element_text(size=25)
        ,axis.text=element_text(size=20)
        ,legend.title=element_text(size=20)
        ,legend.text=element_text(size=20))
plot(g1)

#-- 交互作用項が非有意かどうかで平行性を評価。
ANCOVA0 <- lm(CHG ~ BASE + TRT01AF + BASE*TRT01AF, data=ADS);
summary(ANCOVA0)
car::Anova(ANCOVA0)

#-- ANCOVA
ANCOVA1 <- lm(CHG ~ BASE + TRT01AF, data=ADS)
(res <- summary(ANCOVA1))
car::Anova(ANCOVA1) #-- Type 2 平方和

#-- 共通の回帰直線の重心（総平均）
(b0.adj <- (res$coefficients[1] * nrow(ADS.A) + (res$coefficients[1]+res$coefficients[3]) * nrow(ADS.B)) / nrow(ADS))
(yint <- b0.adj + mean(ADS$BASE) * res$coefficients[2])

#-- 共通傾きでの回帰直線（水準毎）
ANCOVA1.A <- predict(ANCOVA1, newdata=ADS.A, se=F)
ANCOVA1.B <- predict(ANCOVA1, newdata=ADS.B, se=F)

#-- 共通の回帰直線と共通傾きの回帰直線（水準毎）をオーバーラップ
g2 <- g1 +
  geom_smooth(data=ADS, mapping=aes(x=BASE, y=ANCOVA1$fitted.values), stat="smooth", method="lm", se=F, color="darkgreen") +
  geom_smooth(data=ADS.A, mapping=aes(x=BASE, y=ANCOVA1.A), stat="smooth", method="lm", se=F, color="red") +
  geom_smooth(data=ADS.B, mapping=aes(x=BASE, y=ANCOVA1.B), stat="smooth", method="lm", se=F, color="blue") 
  geom_hline(yintercept=yint, lty=2, color="gray60") +
  geom_vline(xintercept=mean(ADS$BASE), lty=2, color="gray60")
plot(g2)

Rの ANCOVAテーブル

> #-- 交互作用項が非有意かどうかで平行性を評価。
> car::Anova(lm(CHG ~ BASE + TRT01AF + BASE*TRT01AF, data=ADS))
Anova Table (Type II tests)
Response: CHG
              Sum Sq Df F value    Pr(>F)    
BASE         120.596  1 293.113 5.692e-07 ***
TRT01AF       28.413  1  69.058 7.144e-05 ***
BASE:TRT01AF   1.357  1   3.299    0.1122    
Residuals      2.880  7                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> #-- ANCOVA
> car::Anova(ANCOVA1) #-- Type 2 平方和
Anova Table (Type II tests)
Response: CHG
           Sum Sq Df F value    Pr(>F)    
BASE      120.596  1 227.682 3.680e-07 ***
TRT01AF    28.413  1  53.642 8.196e-05 ***
Residuals   4.237  8                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

SASでの実行：

data ADS;
input BASE TRT01AN CHG AVAL 8. @@;
cards;
21 0 -7 14
15 0 -2 13 
18 0 -5 13
16 0 -4 12
26 0 -12 14
25 1 -15 10
22 1 -12 10
21 1 -12 9
16 1 -6 10
17 1 -7 10
18 1 -7 11
;run;
proc glm data=ADS;
 class TRT01AN; /* 要因を指定 */
 model CHG = TRT01AN BASE / ss1 ss2 ss3 e solution;
 lsmeans TRT01AN / cl pdiff=control('0');
run;

f:id:cochineal19:20210223215441p:plain
f:id:cochineal19:20210223112602p:plain

プログラムコード

■ Rのコード

#-- 交互作用項が非有意かどうかで平行性を評価。
ANCOVA.0 <- lm(Y ~ X1 + C1 + X1*C1, data=ADS)
summary(ANCOVA.0)
car::Anova(ANCOVA.0)

#-- ANCOVA
ANCOVA.1 <- lm(CHG ~ BASE + TRT01AF, data=ADS)
(res <- summary(ANCOVA.1))
car::Anova(ANCOVA.1) #-- Type 2 平方和

■ SASのコード

proc glm data=ADS;
 class X1; /* 要因を指定 */
 model Y = X1 C1;
 lsmeans X1 / cl pdiff=control('XXX'); /* 調整平均 controlでレファレンスを指定*/
 estimate "X1 XXX vs. YYY" X1 -1 1; /* 対比を用いる場合 */
run;

■ Pythonのコード