【統計】共分散分析(ANCOVA)

本記事では、共分散分析について記載します。

【目次】

 

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

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

 

帰無仮説、対立仮設


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

 

計算式等


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

以下、回帰式を記載しますが、重回帰分析の記事 と同じです。
※カテゴリー変数1つ(2水準)、連続変数(共変量)1つの例としています。

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

■ 共通の回帰直線(共通の条件付き期待値)
\quad \widehat{yc}_{i} = \gamma_{0} + \gamma_{1}c_{i} + \gamma_{2} x_{i}
 ※ c:共変量、x:要因
 ※ 共変量 c の傾きを一定として、要因 x による違いを表現した直線。
 ※ 例えば、要因の水準がA薬・B薬の2つなら、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は0or1のダミー変数)、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}  
   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}  
 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}(2本の回帰直線の偏差平方和)に分解して 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}  

CBASE XTRT01AN YCHG 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

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

CBASE XTRT YCHG 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
   R回帰直線(水準毎) vs. 共通傾きでの回帰直線(水準毎) 1.357 2 0.679 1.4139 0.3140
   e観測値 vs. 回帰直線(水準毎) 2.880 6 0.480
 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%信頼区間
YCHG の平均 -6.000 -9.833 -3.833 -8.9349 1.2682
YCHG の調整平均(LS mean) -6.323 -9.564 -3.240 -4.2608 -2.2202

今回の架空データでは、通常の平均の差の信頼区間は0を挟むのに対し、調整平均では信頼区間の幅が狭まり、0を挟まなくなったことが分かります(信頼区間下限でも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のコード

整備中

 

雑談


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

本記事の架空データでの例:
① CHG=BASE + TRT01AN + BASE*TRT01AN を実行する。
② BASE*TRT01AN が非有意なら、CHG=BASE + TRT01AN のモデルでANCOVAを実行する。

 

参考


統計学(出版:東京図書), 日本統計学会編

多変量解析実務講座テキスト, 実務教育研究所

https://www.sas.com/content/dam/SAS/ja_jp/doc/event/sas-user-groups/usergroups09-a-17h.pdf

http://nfunao.web.fc2.com/files/R-intro/R-stat-intro_06.pdf

 


サイトマップ

cochineal19.hatenablog.com

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。