# Exercise 5.9 in Davis, 2002
# Continuation of Exercises 2.4 and 4.3

week <- c(0:4)
m <- matrix(scan("d:\\data.ext\\davis\\box.dat"),ncol=7,byrow=T)
group <- m[,1]
rat <- m[,2]+m[,1]*10 # unique rat IDs
y <- m[,3:7]

library(nlme)
box <- balancedGrouped( y ~ week|rat, matrix(m[,3:7],nrow=27,ncol=5,dimnames=list(rat,week)),
                           labels=list(y="Weight of rats in 4 weeks"))
Group <- rep(group,rep(5,27))

# ordinary ANOVA, with a minor deviation in the F for week
anova(lme(y~as.factor(week)*as.factor(Group),random=~1|rat,data=box))

library(MASS)
tr <- function(M) sum(diag(M))
x <- matrix(c(rep(1,10),rep(0,27),rep(1,7),rep(0,27),rep(1,10)),nrow=27,ncol=3,byrow=FALSE)
bhat <- ginv(t(x) %*% x) %*% crossprod(x,y) # t(x) %*% y
n <- 27 # no. of subjects
p <- 5-1 # no. of time points -1
rX <- 3 # rank of X
s <- crossprod(y - x %*% bhat)/(n-rX)
c <- t(poly(week,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 condition not supported by the data but epsilon-correction
# does not affect conclusions of significance tests (they are significant anyway).