myscratchpad

#　対応のない1要因
# データ生成。テクニカルブックp. 87
dat <- data.frame(a = factor(c(rep("a1",8), rep("a2",8), rep("a3",8), rep("a4",8))), result = c(9,7,8,8,12,11,8,13, 6,5,6,3,6,7,10,9, 10,13,8,13,12,14,14,16, 9,11,13,14,16,12,15,14))

# 平均と合計
aggregate(dat[2], list(dat[,1]), mean)
aggregate(dat[2], list(dat[,1]), sum)

# lmで分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(result~a, dat)
# carのAnovaで分析
library(car)
Anovares <- Anova(lmres, type=3)
Anovares

# aovと比較
summary(aov(result~a, dat))


#　対応のある1要因
# データ生成。テクニカルブックp. 92
dat <- data.frame(s=factor(1:8), a1=c(9,7,8,8,12,11,8,13), a2=c(6,5,6,3,6,7,10,9), a3=c(10,13,8,13,12,14,14,16), a4=c(9,11,13,14,16,12,15,14))

# 平均と合計
sapply(dat[2:5], mean)
sapply(dat[2:5], sum)

# 分析
options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(cbind(a1, a2, a3, a4)~1, dat)
# carのAnovaで分析
library(car)
afact <- factor(c("a1", "a2", "a3", "a4"))
idat <- data.frame(afact)
Anovares <- Anova(lmres, idata=idat, idesign=~afact, type=3)
summary(Anovares)
(AnovaTable <- Anovastat(Anovares))

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