# Exercise 4.8 in Davis, 2002

week <- c(0:8)
m <- matrix(scan("d:\\data.ext\\davis\\woolson1.dat"),ncol=11,byrow=T)
subject <- m[,2]
group <- m[,1]
y <- m[,3:11]
# ybase <- m[,3]
# yred <- m[4:11]

# code for graphics, using nlme and lattice libraries
library(nlme)
BPRS <- balancedGrouped( y ~ week|subject, matrix(m[,3:11],nrow=40,ncol=9,dimnames=list(subject,week)),
                           labels=list(y="BPRS measurements in 8 weeks"))
Group <- rep(group,rep(9,40))
# profile plots
plot(BPRS, outer=~as.factor(Group), aspect=1)
# mean plot
BPRS.means <- aggregate(BPRS$y,list(BPRS$week,Group),mean,na.rm=TRUE)
names(BPRS.means)[1]<-"Week"
BPRS.means$Week <- as.numeric(levels(BPRS.means$Week))[BPRS.means$Week]
names(BPRS.means)[2]<-"Group"
library(lattice)
xyplot(x~Week|Group,BPRS.means,ylab="Mean BPRS measurements in 8 weeks")

# (a) test of linearity using divided differences
c <- matrix(0,nrow=7,ncol=9)
d <- rep(1,8)
for (i in 1:7) {
  c[i,i] <- -d[i]
  c[i,i+1] <- d[i]+d[i+1]
  c[i,i+2] <- -d[i+1]
  }
nlinddiff <- y %*% t(c)
summary(manova(nlinddiff~as.factor(group)-1),test="Wilk")
group1 <- 2-group
group2 <- group-1
summary(manova(nlinddiff~group1+group2-1),test="Wilk")

# (b) assessment of linearity using orthogonal polynomials
# orthogonal T
orth <- poly(week,8)
t <- matrix(c(rep(1/3,9),orth),nrow=9,ncol=9,byrow=TRUE)
z <- y %*% t(t)
summary.aov(manova(z~as.factor(group)-1))
summary(manova(z[,5:9]~as.factor(group)-1),test="Wilk") # quartic and higher order
summary(manova(z[,4:9]~as.factor(group)-1),test="Wilk") # cubic and higher order
summary(manova(z[,3:9]~as.factor(group)-1),test="Wilk") # quadratic and higher order
summary.aov(manova(z~group1+group2-1))
summary(manova(z[,5:9]~group1+group2-1-1),test="Wilk") # quartic and higher order
summary(manova(z[,4:9]~group1+group2-1-1),test="Wilk") # cubic and higher order
summary(manova(z[,3:9]~group1+group2-1-1),test="Wilk") # quadratic and higher order
# quadratic model appropriate despite weak significance for cubic term in group 1

# (c) quadratic growth model, T=orthonormal
orth <- poly(week,2)
t <- matrix(c(rep(1/3,9),orth),nrow=3,ncol=9,byrow=TRUE)
library(MASS)
x <- matrix(c(group1,group2),nrow=40,ncol=2,byrow=FALSE)
bhat <- ginv(t(x) %*% x) %*% crossprod(x,y) # t(x) %*% y
s <- crossprod(y - x %*% bhat)/(40-2)
s.inv <- ginv(s)
z <- y %*% s.inv %*% t(t) %*% ginv(t %*% s.inv %*% t(t))
summary.aov(manova(z~as.factor(group)-1))
summary(manova(z~as.factor(group)-1),test="Wilk")
summary.aov(manova(z~group))
summary(manova(z~group),test="Wilk")

# (d) quadratic growth model, T=natural
t <- matrix(c(rep(1,9),week,week^2),nrow=3,ncol=9,byrow=TRUE)
z <- y %*% s.inv %*% t(t) %*% ginv(t %*% s.inv %*% t(t))
summary.aov(manova(z~as.factor(group)-1))
summary(manova(z~as.factor(group)-1),test="Wilk")
summary.aov(manova(z~group))
summary(manova(z~group),test="Wilk")
coef(manova(z~as.factor(group)-1))

# (e)
# summary statistic: slope
slope <- c()
for (i in 1:40) 
  slope=c(slope,coef(lm(y[i,]~week))[2])
names(slope) <- NULL
aggregate(slope,list(group),mean,na.rm=TRUE)
aggregate(slope,list(group),sd,na.rm=TRUE)
aggregate(slope,list(group),median,na.rm=TRUE)
t.test(slope ~ group)
wilcox.test(slope ~ group)
# analysis summary statistic insufficient to show any group differences

# comparison of profiles using T^2 statistic
summary(manova(y~group))
diff <- matrix(0,nrow=40,ncol=8)
for (i in 1:8)
  diff[,i] <- y[,i]-y[,i+1]
summary(manova(diff~group))
# profile analysis insufficient to detect differences in profiles