共分散分析, ANCOVA, アンコバ

これは僕が作った架空の例題です。結果についてはご意見があると思いますが、ANCOVAの理解のための例題になっておりますのでご了承ください。

例題） 歩行可能な脳卒中患者を対象にした後ろ向き観察研究です。歩行練習A（10名）と歩行練習B（10名）を受けた20名の10m歩行時間の短縮時間を調べました。

データの準備と要約
2群の差の検定
Rで共分散分析を実行
非並行性の検定（交互性の検定）
両群に共通の回帰直線
共通回帰の検定

データの準備と要約

データセットancovaをdatに格納します（ファイルの読み込み）

ancova データのダウンロード

dat <- read.csv("ancova.csv", header=T, fileEncoding = "UTF-8")

歩行練習はfactor型に置換しておきましょう

dat$歩行練習 <- as.factor(dat$歩行練習)

datの先頭部分を見てみましょう

head(dat)

要約

summary(dat)

歩行練習の各群の要約

色々な方法があると思いますが、ここではパッケージを使用しない方法で

summary(dat[dat$歩行練習 == "A", 2:4])
summary(dat[dat$歩行練習 == "B", 2:4])

2群の差の検定

箱ひげ図

boxplot(短縮時間 ~ 歩行練習, data=dat)

歩行練習Aと歩行練習Bでは随分差があるようです（歩行練習A＞歩行練習B）

t検定

2群の差の平均値を比較検定してみます

 t.test(短縮時間 ~ 歩行練習, data=dat)

歩行練習の違いにより歩行時間短縮が有意に異なることが分かりました

練習Aの方が10m歩行時間短縮には有効である（95%信頼区間は2.76～6.16）・・・と単純に結果を出してもよいでしょうか？

練習前の10m歩行時間が影響しており、もともと時間がかかる人は短縮時間も長くなるのでは・・・？

その疑問に答えるために共分散分析を行います

Rで共分散分析を実行

歩行練習AとBをグラフで確認

fig1 <- function() 
    {
    pchAB <- ifelse(
        dat$歩行練習== "A", 19, 21
        )
    plot(
        dat$練習前, dat$短縮時間,
        pch=pchAB, cex=1.5 ,
        xlab="練習前の10ｍ歩行時間",  ylab="短縮時間"
        )
    legend(
        "topleft",
        legend = c("練習A", "練習B"),
        pch = c(19, 21)
        )
    }
fig1()

練習Aは運動機能が低下している（歩行速度が遅い）患者に適応されていたようです

このような場合には、歩行練習前の10m歩行時間を考慮した解析が必要になります（これが共分散分析になります）

歩行練習前の10m歩行時間を共変量として（歩行時間を考慮して）回帰分析を実施します

Rで実施する場合は、重回帰分析と同じ方法です

fit <- lm(短縮時間 ~ 歩行練習 + 練習前,  data=dat)
summary(fit)
confint(fit)

結果の解釈

練習前の時間を共変量に入れた（調整した）場合、歩行練習AとBの差は2.32（95%信頼区間：0.63～4.02）となりました

練習前の時間を考慮することで、さきほどのｔ検定の結果（95%信頼区間：2.76～6.16）よりも差が小さいことが理解できます

2群の差の検定だけでは読み取れなかった違いを吟味することができました。結果だけでよければここまでです。次からはもう少し詳しく勉強してみます。

非並行性の検定（交互性の検定）

各群の回帰直線

#各グループのデータセットを作成
setA <- subset(dat, dat$歩行練習 == "A")
setB <- subset(dat, dat$歩行練習 == "B")

#各グループの回帰直線
fitA <- lm(短縮時間 ~ 練習前, data=setA)
fitB <- lm(短縮時間 ~ 練習前, data=setB)

#作図
fig3 <- function() 
    {
    fig1()
   
    lines(
        range(setA$練習前),
        fitA$coef[1]+fitA$coef[2]*range(setA$練習前),
        col="red"
     )

    lines(
        range(setB$練習前),
        fitB$coef[1]+fitB$coef[2]*range(setB$練習前),
        col="red"
    )
}

fig3()

非平行性の検定のp値が大きい場合には、両群の回帰直線の傾きが非並行ということになります

両群の共通回帰直線は無意味になり、それぞれの回帰式を吟味しなければなりません

非平行性の検定のp値が小さい場合には、両群の回帰直線の傾きが並行ということになります

この場合には共通する傾きを利用した回帰直線で検証することとなります

回帰分析と分散分析で並行性を検定してみます

fit2 <- lm(短縮時間 ~ 歩行練習*練習前,  data=dat)
summary(fit2)

p>0.05という結果になり交互作用は有意ではありません

この結果から歩行練習Aと歩行練習Bの回帰直線を求めることができます

歩行練習Aの切片$=-3.532$

歩行練習Bの切片$=-3.532+4.4047=0.8727$

歩行練習Aの傾き$=0.4323$

歩行練習Bの傾き$＝0.4323-0.2791=0.1532$

結果

歩行練習A$=-3.532+0.4323*$歩行練習前の10ｍ歩行時間

歩行練習B$=0.8727+0.1532*$歩行練習前の10ｍ歩行時間

注意）回帰直線は描けますが、外挿の予測はできません（結果を求めた標本の範囲内の予測のみ）

anova(fit2)

もちろん同じｐ値になりますが、p>0.05という結果になり交互作用は有意ではありません

（この検定結果で両群の回帰直線が非平行性を決める場合もあります）

交互作用項の平均平方和は各群の回帰式と共通する傾きをもつ回帰式の推定値から求めることができます（詳細は省略しますが、参考図のみ載せておきます）

両群に共通の回帰直線

両群の回帰直線が非並行性の検定で有意でなければ、共通した回帰直線が意味をもちます

全体の回帰直線（赤）
両群に共通した回帰直線（青）･･･fitで求めた直線
共通した傾きの歩行練習A群の回帰直線（緑）
共通した傾きの歩行練習B群の回帰直線（紫）

全て先ほどの回帰分析（fit）から求めることができます

比較のために群分け無しの全体の回帰直線（赤）も描いてみます

#全体の回帰分析（群分け無し）
all <- lm(短縮時間 ~ 練習前,  data=dat)
#全体の回帰分析の傾き
all$coef[2]
#全体の回帰分析の切片
all$coef[1]

fig2 <- function() 
{
    fig1()

#全体の回帰直線
    lines(
        range(dat$練習前), 
        all$coef[1] + all$coef[2]*range(dat$練習前) , 
        col="red"
    )

#共通回帰式
    lines(
        range(dat$練習前),
        (fit$coef[2] + fit$coef[1]*2)/2 + fit$coef[3]*range(dat$練習前),
        col="blue"
    )
    
#歩行練習A群の回帰式（傾きは共通回帰式）
    xa <- dat[dat$歩行練習=="A",  "練習前"]
    lines(
        range(xa),
        fit$coef[1] + fit$coef[3]*range(xa),
        col="green"
    )
   
#歩行練習B群の回帰式（傾きは共通回帰式）
    xb <- dat[dat$歩行練習=="B",  "練習前"]
    lines(
        range(xb),
        fit$coef[2] + fit$coef[1] + fit$coef[3]*range(xb),
        col="purple"
    )
}

fig2()

両群の回帰直線は外挿部分を除外してます

参考

両群に共通した回帰直線を求めるために以下のような作業を行います

まず歩行練習A群の平均を全体の平均までの移動させます（→）

次に、データを同じ距離を平行移動させます（→）

歩行練習B群も同じ操作を行い、全て移行させたデータから求めた回帰直線の傾きが共通した傾きとなります

共通回帰の検定

上のグラフの回帰直線（青）の検定

共通回帰直線の傾きが各群のもつ傾きと比較して意味を持つかどうかを検定します

summary(fit)

F値＝23.28, p<0.05でこのモデル全体が有意であることが分かりました

「練習前の10m歩行時間を歩行時間の短縮に影響する因子として考慮しなければならない」ということが言えます

F値の求め方（詳細は重回帰分析をご覧ください）

各群の共通回帰（緑と紫）から得られる推定値と各群の平均値との差の平均平方和を残差の平均平方和で除したF値で検定します

（平方和については分散分析、重回帰分析で解説してます）

共通回帰のF値が大きければ共通回帰が意味を持つことになります

小さい場合には、共通回帰の傾きが0に近いことを意味します

#共分散分析（fit）から得られた短縮時間の予測値
yhat <- predict(fit)

#回帰式から得た予測値と実際の値の平均値との差の二乗和
Sreg <- sum((yhat - mean(dat$短縮時間))^2)

#残差平方和
Sres <- sum((yhat - dat$短縮時間)^2)

#F値=(Sreg/自由度)÷(Sres/自由度)
F <- (Sreg/2)/(Sres/(20-2-1))
print(F)

#F検定
pf(F, 2, 17, lower.tail=FALSE)