gistfile1.r

## Maxwell & Delaney p. 574
## data from RBF-pq design
md574 <- matrix(c(
420, 420, 480, 480, 600, 780,
420, 480, 480, 360, 480, 600,
480, 480, 540, 660, 780, 780,
420, 540, 540, 480, 780, 900,
540, 660, 540, 480, 660, 720,
360, 420, 360, 360, 480, 540,
480, 480, 600, 540, 720, 840,
480, 600, 660, 540, 720, 900,
540, 600, 540, 480, 720, 780,
480, 420, 540, 540, 660, 780), ncol=6, byrow=TRUE)

dfMD574 <- data.frame(
  rt   =as.vector(md574),
  id   =factor(rep(1:10, 6)),
  deg  =factor(rep(rep(c(0,4,8), c(10, 10, 10)), 2)),
  noise=factor(rep(c("no.noise", "noise"), c(30, 30))))

## average over second within factor (noise)
## -> first group of SPF-p.q design
dfYoung <- aggregate(rt ~ id + deg, FUN=mean, data=dfMD574)
md593a  <- matrix(dfYoung$rt, ncol=3)

## Maxwell & Delaney p. 593
## -> second group of SPF-p.q design
md593b <- matrix(c(
420, 570, 690,
600, 720, 810,
450, 540, 690,
630, 660, 780,
420, 570, 780,
600, 780, 870,
630, 690, 870,
480, 570, 720,
690, 750, 900,
510, 690, 810), ncol=3, byrow=TRUE)

Nj <- 10
P  <- 2
Q  <- 3

dfMD <- data.frame(
  rt =c(as.vector(md593a), as.vector(md593b)),
  id =factor(c(rep(1:10, times=Q), rep(11:20, times=Q))),
  age=factor(rep(c("young", "old"), each=Q*Nj)),
  deg=factor(rep(c(0,4,8,0,4,8), each=Nj)))

summary(aov(rt ~ age*deg + Error(id/deg), data=dfMD))

## Anova()
ageW  <- factor(rep(c("young", "old"), each=10))
fitPQ <- lm(rbind(md593a, md593b) ~ ageW)
inPQ  <- expand.grid(deg=gl(Q, 1))

library(car)
summary(Anova(fitPQ, idata=inPQ, idesign=~deg), multivariate=FALSE)

## manually
Mjk   <- ave(dfMD$rt,          dfMD$age, dfMD$deg, FUN=mean)
Mi    <- ave(dfMD$rt, dfMD$id,                     FUN=mean)
Mj    <- ave(dfMD$rt,          dfMD$age,           FUN=mean)
dfMD$IDxIV <- dfMD$rt - Mi - Mjk + Mj
(SSE  <- sum(dfMD$IDxIV^2))
dfSSE <- (Nj*P - P) * (Q-1)                     # df interaction = df error term
(MSE  <- SSE / dfSSE)                           # mean square error

## epsilon corrections
errMat <- data.matrix(unstack(dfMD, IDxIV ~ deg))
Serr   <- cov(errMat)
(epsGG <- (1 / (Q-1)) * sum(diag(Serr))^2 / sum(Serr^2))

dfId   <- (Nj-1) * P
(epsHF <- ((dfId + 1) * (Q-1) * epsGG - 2) / ((Q-1) * (dfId - (Q-1)*epsGG)))

## contrast
Mj     <- tapply(dfMD$rt, dfMD$deg, FUN=mean)   # group means within factor
Nj     <- table(dfMD$deg)                       # cell sizes
cntr   <- rbind(c(1, -2,  1))                   # contrast vector
psiHat <- cntr   %*% Mj                         # estimate psi-hat
lenSq  <- cntr^2 %*% (1/Nj)                     # squared length contrast vector
tStat  <- psiHat / sqrt(lenSq*MSE)              # t statistic
pVal   <- 2*(1-pt(abs(tStat), dfSSE))           # p-value
data.frame(psiHat, tStat, pVal)