# Additional example for Chapter 5 in Davis, 2002; Orthodont data in Pinheiro & Bates
# Continuation of Examples 3.3, 4.1 and 4.3

age <- c(8,10,12,14)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\pr.dat"),ncol=6,byrow=T)
sex <- m[,1]
subject <- m[,2]
y <- m[,3:6]

library(nlme)
data(Orthodont)
anova(lme(distance~as.factor(age)*Sex,random=~1|Subject,data=Orthodont))

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