# Exercise 6.7 in Davis, 2002
# Continuation of Exercise 4.8

week <- c(0:8)
m <- matrix(scan("c:\\data.ext\\davis\\woolson1.dat"),ncol=11,byrow=T)
group <- m[,1]
subject <- m[,2]+100*group
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 <- as.factor(rep(group,rep(9,40)))

# (a)
# model 1: unstructured, natural time scale parametrization of quadratic model
# (note that orthogonal polynomials leads to same tests, despite different coefficients)
BPRS.unstr1 <- gls(y~Group*week+Group*I(week^2),correlation=corSymm(form=~1|subject),weights=varIdent(form=~1|week),data=BPRS,method="ML")
summary(BPRS.unstr1) # note (almost) same -2log likelihood as in SAS 
anova(BPRS.unstr1) # note differences to SAS and the strange denominator DF
# reparametrization without intercept (to get estimates)
BPRS.unstr2 <- update(BPRS.unstr1,y~Group-1+Group:week+Group:I(week^2))
coef(BPRS.unstr3) # identical to SAS
# extra models for combined tests
BPRS.unstr3 <- update(BPRS.unstr1,y~Group+week+I(week^2))
BPRS.unstr4 <- update(BPRS.unstr1,y~week+I(week^2))
anova(BPRS.unstr1,BPRS.unstr3) # likelihood-ratio test, not Wald tests based on F-dist
anova(BPRS.unstr1,BPRS.unstr4)

# (b) model with common intercept but different slopes (lin+quad)
BPRS.unstr5 <- update(BPRS.unstr1,y~week+Group:week+I(week^2)+Group:I(week^2))

# (c)
# banded structure not readily available in R
BPRS.ar1 <- update(BPRS.unstr5,weight=NULL,correlation=corAR1(form=~1|subject))
BPRS.cs <- update(BPRS.unstr5,weight=NULL,correlation=corCompSymm(form=~1|subject))
BPRS.ind <- update(BPRS.unstr5,weight=NULL,correlation=corIdent(form=~1))
anova(BPRS.unstr5,BPRS.ar1,BPRS.cs,BPRS.ind)

# extra models
BPRS.unstr.hom <- update(BPRS.unstr5,weight=NULL)