# Example 3.2 (Table 3.1) in Davis, 2002

age <- c(8,8.5,9,9.5)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\eg.dat"),ncol=5,byrow=T)
boy <- m[,1]
y <- m[,2:5]

# code for graphics, using nlme and lattice libraries
library(nlme)
eg <- balancedGrouped( y ~ age|boy, matrix(m[,2:5],nrow=20,ncol=4,dimnames=list(boy,age)),
                           labels=list(y="Ramus bone heights (mm)"))
# profile plot
plot(eg, aspect=1)
# mean plot
eg.means <- aggregate(eg$y,list(eg$age),mean,na.rm=TRUE)
names(eg.means)[1]<-"Age"
eg.means$Age <- as.numeric(levels(eg.means$Age))[eg.means$Age]
library(lattice)
xyplot(x~Age,eg.means,ylab="Ramus bone heights (mm)")

# test for constant values over time
diff <- matrix(c(y[,1]-y[,2],y[,2]-y[,3],y[,3]-y[,4]),nrow=20,ncol=3,byrow=FALSE)
summary(manova(diff~1),intercept=TRUE)

# test for linearity using quadratic and cubic contrasts
nlincon <- matrix(c(y[,1]-y[,2]-y[,3]+y[,4],-y[,1]+3*y[,2]-3*y[,3]+y[,4]),nrow=20,ncol=2,byrow=FALSE)
summary(manova(nlincon~1),intercept=TRUE)

# test for linearity using divided differences
nlinddiff <- matrix(c(-y[,1]+2*y[,2]-y[,3],-y[,2]+2*y[,3]-y[,4]),nrow=20,ncol=2,byrow=FALSE)
summary(manova(nlinddiff~1),intercept=TRUE)