理学療法士がまとめたノート

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

共分散分析 ANCOVA

サンプル(勝手に作ったサンプルです。実存しません。)
治療Aと治療Bの降圧効果に関する検証

治療 <- c("B","A","B","A","A","B","B","A","A","B","A","B","B","B","A","B",
    "B","B","A","A")
治療前BP <- c(160,135,177,141,142,155,175,145,149,155,135,156,170,
    150,138,160,150,180,155,140)
治療後BP <- c(148,133,160,138,139,145,157,142,142,145,132,145,157,
    142,134,150,140,160,149,136)
前後差 <- 治療前BP - 治療後BP
dat1 <- data.frame(治療, 治療前BP, 治療後BP, 前後差)
dat1$治療 <- as.factor(dat1$治療)
head(dat1)

  治療 治療前BP 治療後BP 前後差
1    B      160      148     12
2    A      135      133      2
3    B      177      160     17
4    A      141      138      3
5    A      142      139      3
6    B      155      145     10

f:id:yoshida931:20210529081839p:plain:w500
治療Bが治療Aよりも降圧効果があるようです!平均値の差の検定をやってみましょう。

t.test(前後差~治療, data=dat1)
t = -6.7187, df = 13.82, p-value = 1.047e-05
95%信頼区間: -11.543307  -5.951643
A群とB群の平均値
3.888889       12.636364 

差がありました。95%信頼区間から6~11程度の差があるようです。しかし、差が大きいのは治療前BPが高い人では・・・という疑問が残ります。

治療前BPと前後差の散布図と回帰直線

fitAll <- lm(前後差 ~  治療前BP,  data = dat1)
anova(fitAll)
fitAllhat <- fitAll$coef[1] + fitAll$coef[2] *dat1$治療前BP
plot(dat1$治療前BP,dat1$前後差, cex=1.5 ,xlab="治療前BP", ylab="前後差") 
lines(range(治療前BP), fitAll$coef[1] + fitAll$coef[2] * range(治療前BP))#全体の回帰

f:id:yoshida931:20210530104037p:plain

やはり、想定したように治療前の血圧が高い人は治療効果も高くなるようです。この散布図をA群・B群に色分けします。

fig1 <- function() 
    {
    pchAB <- ifelse(dat1$治療 == "A", 19, 21)
    plot(dat1$治療前BP,dat1$前後差, pch=pchAB, cex=1.5 , 
    xlab="治療前BP",  ylab="前後差")
    lines(range(治療前BP),fitAll$coef[1]+fitAll$coef[2]*range(治療前BP),
    col=1)#全体の回帰
    legend("topleft", legend = c("治療A", "治療B"), pch = c(19, 21))
    }
fig1()

f:id:yoshida931:20210530103932p:plain
治療B群にはもともと血圧の高い人が集まっていたようで、治療効果も大きかったようです。治療前BPが交絡変数として考えられます。この治療前BP因子を考慮したうえでA群とB群の差を考えなければいけません。それぞれの群の回帰直線を挿入してみます。

setA <- subset(dat1, dat1$治療 == "A")
setB <- subset(dat1, dat1$治療 == "B")
fitA <- lm(前後差~治療前BP, data=setA)
fitB <- lm(前後差~治療前BP, data=setB)
fig2 <- function() 
    {
    fig1()
    lines(range(治療前BP),fitA$coef[1]+fitA$coef[2]*range(治療前BP),col="blue")
    lines(range(治療前BP),fitB$coef[1]+fitB$coef[2]*range(治療前BP),col="green")
    }
fig2()

f:id:yoshida931:20210530104336p:plain

単純に平均値だけを考えても答えは出ないようです。治療前BPが影響していることは明らかなので、治療前BPを共変量として解析します。(今回は、例として治療前BPを共変量としましたが、実際の解析では多くの変数の中から選択することになります)

統計学入門−第8章 より引用
このように他のデータの影響を考慮して目的のデータを総合的に比較する手法を共分散分析(ANCOVA:analysis of covariance、アンコバ)といい、影響を考慮する変数のことを共変数(covariable)といいます。 この手法は比較したいデータを目的変数にし、影響を考慮するデータを説明変数にした回帰分析と、群を要因にした一元配置分散分析を組み合わせた手法に相当します。 共分散とは回帰分析で重要な役目をする値であり、この手法が共分散分析と呼ばれるのは共分散を利用した分散分析だからです。

単回帰分析の傾きの推定値(最小二乗法)
\frac{\sum{(x_i-\bar{x})(y_i-\bar{y})}}{\sum{(x_i-\bar{x})}^2}

Rを使って解析やってみます。方法は至って簡単で、重回帰分析の説明変数の1つがカテゴリーになっているというだけです。

fit1 <- lm(前後差 ~  治療前BP + 治療, data = dat1)
summary(fit1)

             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -39.79522    4.71524  -8.440 1.75e-07 ***
治療前BP      0.30715    0.03301   9.304 4.41e-08 ***
治療B         2.50511    0.89016   2.814   0.0119 *

共通の傾きは0.30715、2群の切片の差は2.50511。つまり、治療Bの前後差平均値は、治療Bより平均して2.50511高いと推定できます。平均値の差の検定では6~11の差を見積もっていたのですが、誤算であったことが理解できました。

fit1$coef[2]#共通の傾き
fit1$coef[3]+fit1$coef[1]#Bの切片
fit1$coef[1]#Aの切片

fig3 <- function() 
    {
    fig2()
    #共通回帰式
    lines(range(治療前BP), (fit1$coef[3]+fit1$coef[1]*2)/2 + 
    fit1$coef[2]*range(治療前BP), col="red")
    #Aの回帰式(傾きは共通回帰式)
    lines(range(治療前BP), fit1$coef[1]+ fit1$coef[2]*range(治療前BP), 
    lty = 2, col="red")
    #Bの回帰式(傾きは共通回帰式)
    lines(range(治療前BP), fit1$coef[3]+fit1$coef[1]+ 
    fit1$coef[2]*range(治療前BP), lty = 2, col="red")
    }
fig3()

f:id:yoshida931:20210530104724p:plain 黒実線:全体の共通回帰
青実践:Aの回帰直線
緑実線:Bの回帰直線
赤実践:全体の共通回帰
赤点線:各群の共通回帰

共通回帰の傾きの求め方

まず最初に全体の平均と各群のデータの平均値が重なるようにデータを移動します。
イメージ図
f:id:yoshida931:20210529221110p:plain

m治療前BP <- mean(dat1$治療前BP)#治療前BPの全体平均
m前後差 <- mean(dat1$前後差)#前後差の全体平均
ma治療前BP <- tapply(dat1$治療前BP, dat1$治療, mean)[1]#Aの治療前BP平均
mb治療前BP <- tapply(dat1$治療前BP, dat1$治療, mean)[2]#Bの治療前BP平均
ma前後差 <- tapply(dat1$前後差, dat1$治療, mean)[1]#Aの前後差平均
mb前後差 <- tapply(dat1$前後差, dat1$治療, mean)[2]#Bの前後差平均
#移動後の治療前BP、前後差
移動後の治療前BP <- ifelse(dat1$治療=="A", dat1$治療前BP + (m治療前BP - ma治療前BP), dat1$治療前BP + (m治療前BP - mb治療前BP))
移動後の前後差 <- ifelse(dat1$治療=="A", dat1$前後差 + (m前後差 - ma前後差), dat1$前後差 + (m前後差 - mb前後差))
#共通の傾き
lm(移動後の前後差 ~ 移動後の治療前BP)$coef[2]

#  傾き=0.3071539 
# 各群の共通回帰の予測式 y = β0 + 0.3*x にそれぞれの群の平均値と共通回帰の傾きを当てはめて切片β0 を算出する

移動前後の散布図の比較

par(mfrow=c(1,2)) 
plot(治療前BP, 前後差, xlim=c(130,180), ylim=c(0, 20))
lines(range(130:180), fitAll$coef[1]+ fitAll$coef[2]*range(治療前BP))
plot(移動後の治療前BP, 移動後の前後差, xlim=c(130,180), ylim=c(0, 20))
f <- lm(移動後の前後差~移動後の治療前BP)
lines(range(130:180), f$coef[1]+ f$coef[2]*range(130:180))
par(mfrow=c(1,1)) 

f:id:yoshida931:20210529224051p:plain

共通回帰の検定

共通回帰直線の傾きが各群のもつ傾きと比較して意味を持つかどうかを検定します。この結果から「治療前BPを共変量として、前後差に影響する因子として考慮しなければならない」という解析の方向性を示すことになります.
f:id:yoshida931:20210529105002p:plain:w300
各群の共通回帰から得られる推定値と各群の平均値との差の平均平方和を残差の平均平方和で除したF値で検定します。共通回帰のF値が大きければ共通回帰が意味を持つことになる。小さい場合には、共通回帰の傾きが0に近いことを意味します。
F値= (AB群の共通回帰の推定値の平均平方和ー交互作用の平均平方和)÷ 残差平方和

setA <- subset(dat1, dat1$治療 == "A")#A群のみのセット
setB <- subset(dat1, dat1$治療 == "B")#B群のみのセット
fitA <- lm(前後差~治療前BP, data=setA)#A群のみの単回帰
fitB <- lm(前後差~治療前BP, data=setB)#B群のみの単回帰
fitAB <- lm(前後差~治療前BP*治療, data=dat1)#全データの回帰、交互作用あり
S1 <- anova(fitA)$Mean[1] + anova(fitA)$Mean[1]#各群の平方和
S2 <- anova(fitAB)$Mean[3]#交互作用の平均平方和
S3 <- anova(fitAB)$Mean[4]#残差の平方和
Fvalue <- (S1 - S2) / S3
#F検定
#帰無仮説:(AB群の共通回帰の推定値の平均平方和ー交互作用の平均平方和)の分散と残差平方和の分散が等しい
pf(Fvalue, 1, 16, lower.tail=F) 
# 0.00198736

非並行性の検定(交互性の検定)

共通回帰のF値が大きく、非平行性のF値が大きい場合には、両群の回帰直線の傾きが非並行ということになり、両群の共通回帰直線が意味を持つことになります。 共通回帰のF値が小さく、非平行性のF値も小さい場合には、共変量の影響を考慮する必要はなく分散分析で解析します。 ​ f:id:yoshida931:20210530090718p:plain

#交互作用の平均平方和 ÷ 残差平均平方和
f <- S2 / S3
#F検定
#帰無仮説:交互作用の平均平方和の分散と残差平方和の分散が等しい
pf(f, 1, 16, lower.tail=F) 
# 0.05857521

P=0.06ですので、有意水準をどのように設定するかで、A群とB群の非平行性の検定結果は異なります。有意水準は、検定の前に設定しなければなりません。p値から、どのような解析手法にするのか吟味しなければなりません。