# Exercise 3.6 in Davis, 2002

power <- rep(c(6,18,36,60),2)
eye <- c(rep(c("L"),4),rep(c("R"),4))
m <- matrix(scan("c:\\data.ext\\davis\\ch6.dat"),ncol=9,byrow=T)
subject <- m[,1]
y <- m[,2:9]

# code for graphics, using nlme and lattice libraries
library(nlme)
eyes <- balancedGrouped( y ~ power|subject, matrix(m[,2:9],nrow=7,ncol=8,dimnames=list(subject,power)),labels=list(y="Reaction time (ms)"))
Eye <- rep(eye,7)
# profile plot
plot(eyes, inner=~as.factor(Eye), aspect=1)
# mean plot
eyes.means <- aggregate(eyes$y,list(eyes$power,Eye),mean,na.rm=TRUE)
names(eyes.means)[1]<-"Power"
eyes.means$Power <- as.numeric(levels(eyes.means$Power))[eyes.means$Power]
names(eyes.means)[2]<-"Eye"
library(lattice)
xyplot(x~Power|Eye,eyes.means,ylab="Mean Reaction time (ms)")

# (a)
# Hotelling's T^2 cannot be used because t (number of time points) exceeds n (number of subjects)

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

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