# Example 5.1 in Davis, 2002
# Continuation of Example 2.1

temp <- c(-10,25,37,50,65,80)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\deal.dat"),ncol=7,byrow=T)
subject <- m[,1]
y <- m[,2:7]

library(nlme)
deal <- balancedGrouped( y ~ temp|subject, matrix(m[,2:7],nrow=8,ncol=6,dimnames=list(subject,temp)),
                           labels=list(y="Ventilation volumes"))
anova(lme(y~as.factor(temp),random=~1|subject,data=deal)) # ordinary ANOVA

library(MASS)
tr <- function(M) sum(diag(M))
x <- matrix(rep(1,8),nrow=8,ncol=1,byrow=FALSE)
bhat <- ginv(t(x) %*% x) %*% crossprod(x,y) # t(x) %*% y
n <- 8 # no. of subjects
p <- 6-1 # no. of time points -1
rX <- 1 # rank of X
s <- crossprod(y - x %*% bhat)/(n-rX)
c <- t(poly(temp,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))