myscratchpad

dat <- data.frame(
ps=paste("p", 1:27, sep=""),
trt=gl(3, 9, labels=c("a1", "a2", "a3")), # "notreat", "spaced", "massed"
grp=rep(gl(3, 3, labels=c("b1", "b2", "b3")),3), #"washers", "checkers", "seekers"
value=c(5,4,5,4,5,3,3,5,4,5,6,3,4,6,3,12,10,13,10,12,16,11,10,12,13,12,10)
)

agdat <- aggregate(dat[4], list(dat[,2], dat[,3]), mean)

lmres <- lm(value~trt*grp, dat)
library(car)
Anovares <- Anova(lmres)
Anovares
smry <- summary(lmres)
mse <- smry$sigma^2 # 残差の平均平方和
by(dat, dat$trt, function(x) Anova(lm(value~grp, data=x), error=lmres))
by(dat, dat$grp, function(x) Anova(lm(value~trt, data=x), error=lmres))

vnm <- paste(agdat[,1], agdat[,2], sep="")
for (i in 1:9) assign(vnm[i], agdat[,3][i])

cv <- c(1,-1,1,-1) # 対比係数
n <- 3 # 各群の人数

nm <- (sum(cv*c(a1b1, a2b1, a1b3, a2b3)))^2 # (a1b1-a2b1+a1b3-a2b3)　b1水準のa1, a2の差とb3水準のa1, a2の差の差をとって二乗しただけ
dn <- (sum((cv^2)*mse))/n
(fv <- nm/dn)
(tv <- sqrt(fv))
1-pf(fv, df1=1, df2=18) # 分子の自由度は1。分母の自由度は分散分析の誤差自由度
# F(1, 18) = 18.03, p < .001. b1群とb3群はa1 -> a2への反応に違いがあるということ (原文: The significant F indicates that Washers (b1)  and Seekers (b3)  respond differently at the two levels of Factor A (a1 and a2) .
(1-pt(tv, 18))*2


参考:
http://www.psychology.emory.edu/clinical/mcdowell/PSYCH560/factorw.htm
http://core.ecu.edu/psyc/wuenschk/docs30/Factorial-Computations.doc
http://www.sprweb.org/articles/Keselman98.pdf