# Exercise 5.4 in Davis, 2002
# Continuation of Exercise 3.6

power <- rep(c(6,18,36,60),2)
eye <- c(rep(c("L"),4),rep(c("R"),4))
m <- matrix(scan("c:\\data.ext\\davis\\ch6.dat"),ncol=9,byrow=T)
subject <- m[,1]
y <- m[,2:9]

library(nlme)
eyes <- balancedGrouped( y ~ power|subject, matrix(m[,2:9],nrow=7,ncol=8,dimnames=list(subject,power)),labels=list(y="Reaction time (ms)"))
Eye <- as.factor(rep(eye,7))

# hierarchical ANOVA, lacking subject*power term
anova(lme(y~Eye*as.factor(power),random=~1|subject/Eye,data=eyes))
# full ANOVA, with incorrect F-test denominators
anova(lm(y~Eye*as.factor(power)*as.factor(subject)-Eye:as.factor(power):as.factor(subject),data=eyes))

library(MASS)
tr <- function(M) sum(diag(M))
x <- matrix(rep(1,7),nrow=7,ncol=1)
bhat <- ginv(t(x) %*% x) %*% crossprod(x,y) # t(x) %*% y
n <- 7 # no. of subjects
rX <- 1 # rank of X

# sphericity on full set of points within one subject - does not work
p <- 8-1 # no. of time points -1
s <- crossprod(y - x %*% bhat)/(n-rX)
#c <- t(poly(power,p))
c <- t(poly(1:8,p))
r <- c %*% s %*% t(c)
w <- det(r)/(abs(tr(r)/p))^p
rho <- 1-(2*p*p+p+2)/(6*p*(n-rX)) # chi-square approx from SAS/SPSS
x2 <- -rho*(n-rX)*log(w)          # df=p*(p+1)/2-1, but maybe more complex formula
eps.GG <- tr(r)^2/(p*tr(r %*% r))
eps.HF <- (n*p*eps.GG-2)/(p*(n-rX-p*eps.GG))

# sphericity on powers only, same as in Crowder & Hand p.57-58
p <- 4-1 # no. of time points -1
s <- crossprod(y - x %*% bhat)/(n-rX)
c <- t(poly(power,p))
r <- c %*% s %*% t(c)
w <- det(r)/(abs(tr(r)/p))^p
rho <- 1-(2*p*p+p+2)/(6*p*(n-rX)) # chi-square approx from SAS/SPSS
x2 <- -rho*(n-rX)*log(w)          # df=p*(p+1)/2-1, but maybe more complex formula
1-pchisq(x2,p*(p+1)/2-1)
eps.GG <- tr(r)^2/(p*tr(r %*% r))
eps.HF <- (n*p*eps.GG-2)/(p*(n-rX-p*eps.GG))

# sphericity on power*eye interaction, same as in Crowder & Hand p.57-58
p <- 4-1 # no. of time points -1
s <- crossprod(y - x %*% bhat)/(n-rX)
c <- t(poly(power,p))*matrix(c(rep(1,12),rep(-1,12)),nrow=3,ncol=8,byrow=FALSE)
r <- c %*% s %*% t(c)
w <- det(r)/(abs(tr(r)/p))^p
rho <- 1-(2*p*p+p+2)/(6*p*(n-rX)) # chi-square approx from SAS/SPSS
x2 <- -rho*(n-rX)*log(w)          # df=p*(p+1)/2-1, but maybe more complex formula
1-pchisq(x2,p*(p+1)/2-1)
eps.GG <- tr(r)^2/(p*tr(r %*% r))
eps.HF <- (n*p*eps.GG-2)/(p*(n-rX-p*eps.GG))