myscratchpad

carパッケージのAnova関数で混合要因の分散分析
データはデータ解析テクニカルブック
SPFp.qデザイン　データ数が同じ(p.107)

dat <- data.frame(
s=factor(c("s1","s2","s3","s4","s5")),
a=factor(c(rep("a1",5), rep("a2",5))),
b1=c(3,3,1,3,5,3,5,2,4,6), b2=c(4,3,4,5,7,2,6,3,6,4), b3=c(6,6,6,4,8,3,2,3,6,5), b4=c(5,7,8,7,9,2,3,3,4,6)
)
dat
aggregate(dat[3:6], list(dat[,2]), sum)

# メインの分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(cbind(b1, b2, b3, b4) ~ a , contrasts = list(a = contr.sum), data = dat)

bfact <- as.factor(c("b1", "b2", "b3", "b4"))
idat <- data.frame(bfact)
library(car)
Anovares <- Anova(lmres, idata = idat, idesign = ~ bfact, type = 3)
summary(Anovares)

# aovと比較
library(reshape)
dat.aov <- melt(dat, id=c("s", "a"))
names(dat.aov) <- c("s", "a", "b", "result")
print(summary(aov(result ~ a * b + Error(s:a + s:a:b), dat.aov)))

# 作図のため、平均値を出して行列にする
agdat <- aggregate(dat[3:6], list(dat[,2]), mean)
mdat <- as.matrix(agdat[2:5])
# 折れ線で作図
matplot(mdat, type="b", pch=1:4, xaxt="n")
axis(1, at=1:2, labels=c("a1",  "a2"))
legend("topright", legend=c("b1", "b2", "b3", "b4"), lty=1:4, pch=1:4, col=1:4)
# 棒グラフで作図
rownames(mdat) <- c("a1", "a2")
palette(grey.colors(nrow(t(mdat)))) # a要因を横軸にしたいのでtで転置
barplot(t(mdat), beside=T)
legend("topright", rownames(t(mdat)), fill = 1:4)


## 事後検定　テクニカルブックとは違って水準別
## b要因 (被験者内要因) の水準ごと
lmab1 <- lm(b1 ~ a, contrasts = list(a = contr.sum), data = dat)
Anv.ab1 <- Anova(lmab1, type = 3)
lmab2 <- lm(b2 ~ a, contrasts = list(a = contr.sum), data = dat)
Anv.ab2 <- Anova(lmab2, type = 3)
lmab3 <- lm(b3 ~ a, contrasts = list(a = contr.sum), data = dat)
Anv.ab3 <- Anova(lmab3, type = 3)
lmab4 <- lm(b4 ~ a, contrasts = list(a = contr.sum), data = dat)
Anv.ab4 <- Anova(lmab4, type = 3)

Anv.ab1
Anv.ab2
Anv.ab3
Anv.ab4

## 事後検定
## a要因 (被験者間要因) の水準ごと、aの１
levels(dat$a)
data1 <- subset(dat, a=="a1")
data1
lmba1 <- lm(cbind(b1, b2, b3, b4) ~ 1, contrasts=contr.sum, data = data1)
Anv.a1 <- Anova(lmba1, idata = idat, idesign = ~ bfact2, type = 3)
summary(Anv.a1)

## a要因 (被験者間要因) の水準ごと、aの2
data2 <- subset(dat, a=="a2")
data2
lmba2 <- lm(cbind(b1, b2, b3, b4) ~ 1, contrasts=contr.sum, data = data2)
Anv.a2 <- Anova(lmba2, idata = idat, idesign = ~ bfact2, type = 3)
summary(Anv.a2)

##多重比較 単純効果のあったa1で
#holm法
datpa <- stack(data1[3:6])
datpa$s <- rep(data1[,2], 4)
pttres <- pairwise.t.test(datpa$values, datpa$ind, p.adjust.method="holm", paired=TRUE)
pttres

#holm法2
p12 <- t.test(data1$b1, data1$b2, paired = T)
p13 <- t.test(data1$b1, data1$b3, paired = T)
p14 <- t.test(data1$b1, data1$b4, paired = T)
p23 <- t.test(data1$b2, data1$b3, paired = T)
p24 <- t.test(data1$b2, data1$b4, paired = T)
p34 <- t.test(data1$b3, data1$b4, paired = T)
p.adjust(c(p12$p.value,p13$p.value,p14$p.value,p23$p.value,p24$p.value,p34$p.value), method = "holm")