統計学備忘録様 ANCOVA.
dat <- read.csv("ancova.csv", header = TRUE)
dat$歩行練習 <- as.factor(dat$歩行練習)
head(dat)
| 1 |
B |
22.5 |
5.6 |
| 2 |
A |
25.5 |
7.8 |
| 3 |
B |
19.3 |
3.3 |
| 4 |
A |
25.6 |
9.6 |
| 5 |
A |
31.2 |
10.8 |
| 6 |
B |
16.5 |
4.0 |
dim(dat)
## [1] 20 4
par(family = "HiraKakuProN-W3") #日本語フォントの指定
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()
fit <- lm(短縮時間 ~ 歩行練習 + 練習前, data=dat)
summary(fit)
##
## Call:
## lm(formula = 短縮時間 ~ 歩行練習 + 練習前, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0202 -0.7083 -0.0401 1.0303 1.7947
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.61018 2.42519 -0.252 0.80436
## 歩行練習B -2.32654 0.80219 -2.900 0.00996 **
## 練習前 0.32873 0.08474 3.879 0.00121 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.306 on 17 degrees of freedom
## Multiple R-squared: 0.8119, Adjusted R-squared: 0.7898
## F-statistic: 36.69 on 2 and 17 DF, p-value: 6.797e-07
confint(fit)
## 2.5 % 97.5 %
## (Intercept) -5.7268761 4.5065125
## 歩行練習B -4.0190108 -0.6340757
## 練習前 0.1499362 0.5075235
#各グループのデータセットを作成
setA <- subset(dat, dat$歩行練習 == "A")
setB <- subset(dat, dat$歩行練習 == "B")
#各グループの回帰直線
fitA <- lm(短縮時間 ~ 練習前, data=setA)
fitB <- lm(短縮時間 ~ 練習前, data=setB)
#作図
par(family = "HiraKakuProN-W3") #日本語フォントの指定
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()
fit2 <- lm(短縮時間 ~ 歩行練習*練習前, data=dat)
summary(fit2)
##
## Call:
## lm(formula = 短縮時間 ~ 歩行練習 * 練習前, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1653 -0.5919 0.2472 0.6316 2.0641
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.5320 2.8928 -1.221 0.239776
## 歩行練習B 4.4047 4.0947 1.076 0.298011
## 練習前 0.4323 0.1016 4.254 0.000606 ***
## 歩行練習B:練習前 -0.2791 0.1668 -1.673 0.113726
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.242 on 16 degrees of freedom
## Multiple R-squared: 0.8399, Adjusted R-squared: 0.8099
## F-statistic: 27.98 on 3 and 16 DF, p-value: 1.334e-06
anova(fit2)
| 歩行練習 |
1 |
99.458000 |
99.458000 |
64.50441 |
0.0000005 |
| 練習前 |
1 |
25.657365 |
25.657365 |
16.64032 |
0.0008738 |
| 歩行練習:練習前 |
1 |
4.316568 |
4.316568 |
2.79955 |
0.1137261 |
| Residuals |
16 |
24.670067 |
1.541879 |
NA |
NA |
par(family = "HiraKakuProN-W3") #日本語フォントの指定
#全体の回帰分析（群分け無し）
all <- lm(短縮時間 ~ 練習前, data=dat)
#全体の回帰分析の傾き
all$coef[2]
## 練習前
## 0.4972373
#全体の回帰分析の切片
all$coef[1]
## (Intercept)
## -5.978558
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()
