# Exercise 4.3 in Davis, 2002
# Continuation of Exercise 2.4

week <- c(0:4)
#m <- matrix(scan("e:\\vhm\\vhm882\\data\\leppik.dat"),ncol=6,byrow=T)
m <- matrix(scan("d:\\data.ext\\davis\\box.dat"),ncol=7,byrow=T)
group <- m[,1]
rat <- m[,2]
y <- m[,3:7]

# (a) profile analysis
# test for parallel curves
diff <- matrix(c(y[,1]-y[,2],y[,2]-y[,3],y[,3]-y[,4],y[,4]-y[,5]),nrow=27,ncol=4,byrow=FALSE)
summary(manova(diff~as.factor(group)),test="Wilk")

# test for equality between groups, without assuming parallelism
summary(manova(y~as.factor(group)),test="Wilk")
summary.aov(manova(y~as.factor(group)))

# test for equality between weeks, without assuming parallelism
summary(manova(diff~as.factor(group)-1),test="Wilk")
control<- as.numeric(group==1)
thyroxin <- as.numeric(group==2)
thiouracil <- as.numeric(group==3)
summary(manova(diff~control+thyroxin+thiouracil-1))

# (b) growth curve analysis
# orthogonal T
orth <- poly(week,4)
t <- matrix(c(rep(sqrt(0.2),5),orth),nrow=5,ncol=5,byrow=TRUE)
z <- y %*% t(t)
summary.aov(manova(z~as.factor(group)-1))
summary(manova(z[,2:5]~as.factor(group)-1),test="Wilk") # parallel
summary(manova(z[,3:5]~as.factor(group)-1),test="Wilk") # nonlinear
summary(manova(z[,4:5]~as.factor(group)-1),test="Wilk") # high-order
# natural T
t <- matrix(c(rep(1,5),week,week^2,week^3,week^4),nrow=5,ncol=5,byrow=TRUE)
library(MASS)
z <- y %*% ginv(t)
coef(manova(z~as.factor(group)-1))