# Example 5.2 (Table 5.4) in Davis, 2002

fraction <- c(0.9,0.7,0.5,0.3,0.25,0.1)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\crowder.dat"),ncol=8,byrow=T)
group <- m[,1]
subject <- m[,2]+m[,1]*10 # unique subject IDs
y <- m[,3:8]

# code for graphics, using nlme and lattice libraries
library(nlme)
dissolv <- balancedGrouped( y ~ fraction|subject, matrix(m[,3:8],nrow=17,ncol=6,dimnames=list(subject,fraction)),
                         labels=list(y="Logarithmic time to dissolve tablets"))
Group <- rep(group,rep(6,17))
plot(dissolv,outer=~as.factor(Group),aspect=1)

# mean plot
diss.means <- aggregate(dissolv$y,list(dissolv$fraction,Group),mean,na.rm=TRUE)
names(diss.means)[1]<-"Fraction"
diss.means$Fraction <- as.numeric(levels(diss.means$Fraction))[diss.means$Fraction]
names(diss.means)[2]<-"Group"
library(lattice)
xyplot(x~Fraction|Group,diss.means,ylab="Mean Logarithmic time to dissolve tablets")

# ordinary ANOVA but note slight difference in F for fraction
anova(lme(y~as.factor(fraction)*as.factor(Group),random=~1|subject,data=dissolv))

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