library(nlme)

# Section 4.1
data(Orthodont)
fm1Orth.lm <- lm(distance ~ age, Orthodont)
par(mfrow=c(2,2))
plot(fm1Orth.lm)

# Section 4.2.1
fm1Orth.lme <- lme(distance ~ I(age-11), Orthodont, random=~I(age-11)|Subject)
fm2Orth.lme <- update(fm1Orth.lme, fixed=distance~Sex*I(age-11))
summary(fm2Orth.lme)
fitted(fm2Orth.lme, level=0:1)
resid(fm2Orth.lme, level=1) # residuals() works as well
newOrth <- data.frame(Subject=rep(c("M11","F03"),c(3,3)),
Sex=rep(c("Male","Female"),c(3,3)), age=rep(16:18,2))
newOrth
predict (fm2Orth.lme, newdata=newOrth)

data(IGF)
formula(IGF)
summary(IGF)
fm1IGF.lis <- lmList(IGF)
fm1IGF.lme <- lme(fm1IGF.lis)
intervals(fm1IGF.lme)
fm2IGS.lme <- update(fm1IGF.lme, random=pdDiag(~age))
anova(fm1IGF.lme, fm2IGS.lme)

# Section 4.2.2
data(Oats)
fm4OatsB <- lme(yield~nitro, data=Oats, random=list(Block=pdCompSymm(~Variety - 1)))
summary(fm4OatsB)
fm4OatsC <- lme(yield~nitro, data=Oats, random=list(Block=pdBlocked(list(pdIdent(~1),pdIdent(~Variety-1)))))
summary(fm4OatsC)

data(Assay)
formula(Assay)
summary(Assay)
fm1Assay <- lme(logDens~ sample*dilut, data=Assay,
random=pdBlocked(list(pdIdent(~1),pdIdent(~sample-1),pdIdent(~dilut-1))))
fm1Assay
anova(fm1Assay)

# Section 4.2.3
data(Wafer)
formula(Wafer)
fm1Wafer <- lme(current~ voltage + I(voltage^2), Wafer,
random=list(Wafer=pdDiag(~voltage+I(voltage^2)), Site=pdDiag(~voltage+I(voltage^2))))
summary(fm1Wafer)

# Section 4.3.1
plot(fm2Orth.lme, Subject~resid(.), abline=0)
plot(fm2Orth.lme, resid(., type="p") ~ fitted(.)|Sex, id=0.05)
fm3Orth.lme <- update(fm2Orth.lme, weights=varIdent(form=~1|Sex))
fm3Orth.lme
plot(fm3Orth.lme, distance~fitted(.), id=0.05)
qqnorm(fm3Orth.lme, ~resid(.)|Sex)

data(Wafer)
fm1Wafer <- lme(current~ voltage + I(voltage^2), Wafer,
random=list(Wafer=pdDiag(~voltage+I(voltage^2)), Site=pdDiag(~voltage+I(voltage^2))))
plot(fm1Wafer, resid(.)~voltage|Wafer)
plot(fm1Wafer, resid(.)~voltage|Wafer, panel=function(x,y,...){
panel.grid()
panel.xyplot(x,y)
panel.loess(x,y,lty=2)
panel.abline(0,0)})
fm2Wafer <- update(fm1Wafer, fixed=.~.+cos(4.5679*voltage)+sin(4.5679*voltage),
random=list(Wafer=pdDiag(~voltage+I(voltage^2)),Site=pdDiag(~voltage+I(voltage^2))))
plot(fm2Wafer, resid(.)~voltage|Wafer)
plot(fm2Wafer, resid(.)~voltage|Wafer, panel=function(x,y,...){
panel.grid()
panel.xyplot(x,y)
panel.loess(x,y,lty=2)
panel.abline(0,0)})

# Section 4.3.2
qqnorm(fm2Orth.lme, ~ranef(.), id=0.1, cex=0.7)
pairs(fm2Orth.lme, ~ranef(.)|Sex, id=~Subject=="M13")
qqnorm(fm3Orth.lme, ~ranef(.), id=0.1, cex=0.7)
pairs(fm2Wafer, ~ranef(.))