# Example 6.1 in Davis, 2002; Orthodont data in Pinheiro & Bates
# Continuation of Examples 3.3, 4.1 and 4.3

library(nlme)
data(Orthodont)

# model 1: unstructured, age as categorical
Orth.unstr1 <- gls(distance~as.factor(age)*Sex,correlation=corSymm(form=~1|Subject),weights=varIdent(form=~1|age),data=Orthodont,method="ML")
summary(Orth.unstr1)
# same model, with subject random effects separated out
lme(distance~as.factor(age)*Sex,random=~1|Subject,correlation=corSymm(form=~1|Subject),weights=varIdent(form=~1|age),data=Orthodont,method="ML")

# model 2: unstructured, age (centered at 8 yrs) as linear effect
Orth.unstr2 <- update(Orth.unstr1,distance~I(age-8)*Sex)
summary(Orth.unstr2)

# model 3: unstructured, age as linear, parallel curves
Orth.unstr3 <- update(Orth.unstr2,distance~I(age-8)+Sex)
summary(Orth.unstr3)
anova(Orth.unstr1,Orth.unstr2,Orth.unstr3)

# model 4a: unstructured, homogeneous variances
Orth.unstr.hom <- update(Orth.unstr2,weight=NULL)
summary(Orth.unstr.hom)

# model 4 (for 4 time points, equivalent to Toeplitz or banded)
Orth.toep <- update(Orth.unstr.hom,correlation=corARMA(c(0.1,0.1,0.4),form=~1|Subject,p=2,q=1))
summary(Orth.toep)
corMatrix(Initialize(corARMA(c(0.1585499,0.5909483,0.2035497),form=~1|Subject,p=2,q=1),data=Orthodont))[1]

# model 4b: ARMA(1,1)
Orth.arma <- update(Orth.unstr.hom,correlation=corARMA(c(0.1,0.1),form=~1|Subject,p=1,q=1))
summary(Orth.arma)
corMatrix(Initialize(corARMA(c(0.9703622,-0.7162636),form=~1|Subject,p=1,q=1),data=Orthodont))[1]

# model 5: autoregressive (1)
Orth.ar1 <- update(Orth.unstr.hom,correlation=corAR1(form=~1|Subject))
summary(Orth.ar1)
# Orth.ar1b <- update(Orth.unstr.hom,correlation=corARMA(form=~1|Subject,p=1)) # same results
anova(Orth.unstr2,Orth.unstr.hom,Orth.toep,Orth.arma,Orth.ar1)

# model 6: random effects model with random slopes for age
Orth.ranslp <- lme(distance~I(age-1)*Sex,random=~age|Subject,,data=Orthodont,method="ML")
summary(Orth.ranslp)

# model 7: random intercept model
Orth.ran <- update(Orth.ranslp,random=~1|Subject)
summary(Orth.ran)
anova(Orth.ranslp,Orth.ran)

# model 7: compound symmetry
Orth.compsym <- update(Orth.unstr.hom,correlation=corCompSymm(form=~1|Subject))
summary(Orth.compsym)
anova(Orth.unstr2,Orth.unstr.hom,Orth.toep,Orth.arma,Orth.compsym)
anova(Orth.unstr2,Orth.unstr.hom,Orth.toep,Orth.compsym)

# model 8: independence
Orth.ind <- update(Orth.unstr.hom,correlation=corIdent(form=~1))
summary(Orth.ind)