# Exercise 6.2 in Davis, 2002
# Continuation of Examples 3.2 and 4.2

age <- c(8,8.5,9,9.5)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\eg.dat"),ncol=5,byrow=T)
boy <- m[,1]
y <- m[,2:5]

library(nlme)
eg <- balancedGrouped( y ~ age|boy, matrix(m[,2:5],nrow=20,ncol=4,dimnames=list(boy,age)),
                           labels=list(y="Ramus bone heights (mm)"))
# (a)
# model 1: unstructured, age as categorical
eg.unstr1 <- gls(y~as.factor(age),correlation=corSymm(form=~1|boy),weights=varIdent(form=~1|age),data=eg,method="ML")
summary(eg.unstr1)

# model 2: unstructured, age as linear
eg.unstr2 <- update(eg.unstr1, y~age)
summary(eg.unstr2)
# same estimates, but smaller SE's than multivariate approach
anova(eg.unstr1,eg.unstr2)

# (b)
# model 3: unstructured, homogeneous variances
eg.unstr.hom <- update(eg.unstr2,weight=NULL)
summary(eg.unstr.hom)

# model 4 (for 4 time points, equivalent to Toeplitz or banded) - does not work
#eg.toep <- update(eg.unstr.hom,correlation=corARMA(c(0.4,0.4,0.1),form=~1|boy,p=2,q=1))
#summary(eg.toep)
#corMatrix(Initialize(corARMA(c(0.0497636,0.9502312,0.9999999),form=~1|boy,p=2,q=1),data=egodont))[1]
# alternative ARMA(1,1)
eg.arma11 <- update(eg.unstr.hom,correlation=corARMA(c(0.9,0.1),form=~1|boy,p=1,q=1))
summary(eg.arma11)
corMatrix(Initialize(corARMA(c(0.9287301, 0.2334555),form=~1|boy,p=1,q=1),data=eg))[1]

# model 5: autoregressive (1)
eg.ar1 <- update(eg.unstr.hom,correlation=corAR1(form=~1|boy))
summary(eg.ar1)

# model 6: compound symmetry
eg.compsym <- update(eg.unstr.hom,correlation=corCompSymm(form=~1|boy))
summary(eg.compsym)

# model 7: independence
eg.ind <- update(eg.unstr.hom,correlation=corIdent(form=~1))
summary(eg.ind)
anova(eg.unstr2,eg.unstr.hom,eg.arma11,eg.ar1,eg.compsym,eg.ind)

# (c)
# model 8: random effects model with random slopes for age
eg.ranslp <- lme(y~age,random=~age|boy,data=eg,method="ML")
summary(eg.ranslp)

# model 9: random slopes with added ar(1) errors
eg.ranslp.ar1 <- update(eg.ranslp,random=~1|boy,correlation=corAR1(form=~1|boy))
summary(eg.ranslp.ar1)

# model 10: random intercepts with added ar(1) errors
eg.ranint.ar1 <- update(eg.ranslp,correlation=corAR1(form=~1|boy))
summary(eg.ranslp.ar1)
anova(eg.ranslp.ar1,eg.ranslp,eg.ranint.ar1,eg.compsym)

# best model seems to be ar(1) without any need for random slopes or intercepts.