# Exercise 3.9 in Davis, 2002

hour <- rep(c(1,2,3,6),2)
drug <- c(rep(c("A"),4),rep(c("B"),4))
m <- matrix(scan("c:\\data.ext\\davis\\ch2.dat"),ncol=9,byrow=T)
subject <- m[,1]
y <- m[,2:9]

# code for graphics, using nlme and lattice libraries
library(nlme)
drugs <- balancedGrouped( y ~ hour|subject, matrix(m[,2:9],nrow=5,ncol=8,dimnames=list(subject,hour)),labels=list(y="Antibiotic serum level"))
drug <- rep(drug,5)
# profile plot
plot(drugs, inner=~as.factor(drug), aspect=1)
# mean plot
drugs.means <- aggregate(drugs$y,list(drugs$hour,drug),mean,na.rm=TRUE)
names(drugs.means)[1]<-"hour"
drugs.means$hour <- as.numeric(levels(drugs.means$hour))[drugs.means$hour]
names(drugs.means)[2]<-"drug"
library(lattice)
xyplot(x~hour|drug,drugs.means,ylab="Mean Antibiotic serum level")

# (a) drug difference
diffdrug <- matrix(c(y[,1]-y[,5],y[,2]-y[,6],y[,3]-y[,7],y[,4]-y[,8]),nrow=5,ncol=4,byrow=FALSE)
summary(manova(diffdrug~1),intercept=TRUE)

# (b) time differences, both A and B
# Hotelling's T^2 cannot be used the dimension of the contrast matrix (2*3=6) exceeds n (number of subjects)

# (c) time differences A
difftimeA <- matrix(c(y[,1]-y[,2],y[,2]-y[,3],y[,3]-y[,4]),nrow=5,ncol=3,byrow=FALSE)
summary(manova(difftimeA~1),intercept=TRUE)

# (d) time differences B
difftimeB <- matrix(c(y[,5]-y[,6],y[,6]-y[,7],y[,7]-y[,8]),nrow=5,ncol=3,byrow=FALSE)
summary(manova(difftimeB~1),intercept=TRUE)

# (e) linear time effects, using divided differences
c <- matrix(c(-1,2,-1,0,rep(0,4),0,-1,1+1/3,-1/3,rep(0,4),rep(0,4),-1,2,-1,0,rep(0,4),0,-1,1+1/3,-1/3),nrow=4,ncol=8,byrow=TRUE)
divdiff <- y %*% t(c)
summary(manova(divdiff~1),intercept=TRUE)
divdiffA <- y %*% t(c[1:2,])
summary(manova(divdiffA~1),intercept=TRUE)
divdiffB <- y %*% t(c[3:4,])
summary(manova(divdiffB~1),intercept=TRUE)