# Exercise 5.1 in Pinheiro & Bates, 2000

library(nlme)
data(BodyWeight)

# (a)
fm1BW.lme <- lme(weight~Time*Diet, data=BodyWeight, random=~Time)
plot(fm1BW.lme, resid(.)~as.integer(Diet), abline=0)

# (b)
fm1BW.lme.het <- update(fm1BW.lme, weights=varIdent(form=~1|Diet))
summary(fm1BW.lme.het)
intervals(fm1BW.lme.het)
# note error in wording of exercise: one should definitely not use the helmert contrasts
# CI's show that diet 3 differs significantly from diet 1, but diet 2 does not

# (c)
anova(fm1BW.lme,fm1BW.lme.het)
# clear significance based on LR test for different variances at the diets
fm2BW.lme <- update(fm1BW.lme, weights=varPower())
anova(fm1BW.lme.het,fm2BW.lme,test=F)
# log likelihood is slightly lower for fm1BW.lme.het but model has one parameter more,
# therefore fm2BW.lme could be preferred; as an alternative we fit a two-parameter Diet model
BodyWeight$Diet2 <- as.numeric(BodyWeight$Diet==3)
fm1BW.lme.het2 <- update(fm1BW.lme, weights=varIdent(form=~1|Diet2))
summary(fm1BW.lme.het2)
anova(fm1BW.lme.het,fm1BW.lme.het2,fm2BW.lme)
# the two-parameter Diet model has slightly higher logL than the varPower model
# but collapsing Diets 1 and 2 may be biologically questionable

# (d)
fm3BW.lme <- update(fm2BW.lme, corr=corExp(form=~Time))
# same as fm3BW.lme <- update(fm2BW.lme, corr=corCAR1(form=~Time))
fm1BW.gls <- gls(weight~Time*Diet, data=BodyWeight, corr=corCAR1(form=~Time|Rat), weights=varPower())
anova(fm3BW.lme, fm1BW.gls)
# the models ARE nested because fm3BW.lme has the random effects in addition to fm1BW.gls