サンプル（勝手に作ったサンプルです。実存しません。）
治療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
治療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))
やはり、想定したように治療前の血圧が高い人は治療効果も高くなるようです。この散布図を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()
治療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()
単純に平均値だけを考えても答えは出ないようです。治療前BPが影響していることは明らかなので、治療前BPを共変量として解析します。（今回は、例として治療前BPを共変量としましたが、実際の解析では多くの変数の中から選択することになります）
統計学入門−第8章 より引用
このように他のデータの影響を考慮して目的のデータを総合的に比較する手法を共分散分析(ANCOVA：analysis of covariance、アンコバ)といい、影響を考慮する変数のことを共変数(covariable)といいます。 この手法は比較したいデータを目的変数にし、影響を考慮するデータを説明変数にした回帰分析と、群を要因にした一元配置分散分析を組み合わせた手法に相当します。 共分散とは回帰分析で重要な役目をする値であり、この手法が共分散分析と呼ばれるのは共分散を利用した分散分析だからです。
単回帰分析の傾きの推定値（最小二乗法）
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.50511。つまり、治療Bの前後差平均値は、治療Bより平均して2.50511高いと推定できます。平均値の差の検定では6～11の差を見積もっていたのですが、誤算であったことが理解できました。
fit1$coef[2]
fit1$coef[3]+fit1$coef[1]
fit1$coef[1]
fig3 <- function()
{
fig2()
lines(range(治療前BP), (fit1$coef[3]+fit1$coef[1]*2)/2 +
fit1$coef[2]*range(治療前BP), col="red")
lines(range(治療前BP), fit1$coef[1]+ fit1$coef[2]*range(治療前BP),
lty = 2, col="red")
lines(range(治療前BP), fit1$coef[3]+fit1$coef[1]+
fit1$coef[2]*range(治療前BP), lty = 2, col="red")
}
fig3()
黒実線：全体の共通回帰
青実践：Aの回帰直線
緑実線：Bの回帰直線
赤実践：全体の共通回帰
赤点線：各群の共通回帰
共通回帰の傾きの求め方
まず最初に全体の平均と各群のデータの平均値が重なるようにデータを移動します。
イメージ図
m治療前BP <- mean(dat1$治療前BP)
m前後差 <- mean(dat1$前後差)
ma治療前BP <- tapply(dat1$治療前BP, dat1$治療, mean)[1]
mb治療前BP <- tapply(dat1$治療前BP, dat1$治療, mean)[2]
ma前後差 <- tapply(dat1$前後差, dat1$治療, mean)[1]
mb前後差 <- tapply(dat1$前後差, dat1$治療, mean)[2]
移動後の治療前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]
移動前後の散布図の比較
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))
共通回帰の検定
共通回帰直線の傾きが各群のもつ傾きと比較して意味を持つかどうかを検定します。この結果から「治療前BPを共変量として、前後差に影響する因子として考慮しなければならない」という解析の方向性を示すことになります．
各群の共通回帰から得られる推定値と各群の平均値との差の平均平方和を残差の平均平方和で除したF値で検定します。共通回帰のF値が大きければ共通回帰が意味を持つことになる。小さい場合には、共通回帰の傾きが0に近いことを意味します。
F値= （AB群の共通回帰の推定値の平均平方和ー交互作用の平均平方和）÷ 残差平方和
setA <- subset(dat1, dat1$治療 == "A")
setB <- subset(dat1, dat1$治療 == "B")
fitA <- lm(前後差~治療前BP, data=setA)
fitB <- lm(前後差~治療前BP, data=setB)
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
pf(Fvalue, 1, 16, lower.tail=F)
非並行性の検定（交互性の検定）
共通回帰のF値が大きく、非平行性のF値が大きい場合には、両群の回帰直線の傾きが非並行ということになり、両群の共通回帰直線が意味を持つことになります。 共通回帰のF値が小さく、非平行性のF値も小さい場合には、共変量の影響を考慮する必要はなく分散分析で解析します。
f <- S2 / S3
pf(f, 1, 16, lower.tail=F)
P=0.06ですので、有意水準をどのように設定するかで、A群とB群の非平行性の検定結果は異なります。有意水準は、検定の前に設定しなければなりません。ｐ値から、どのような解析手法にするのか吟味しなければなりません。
