library(nlme)

# Section 1.1
data(Rail)
plot(Rail)
fm1Rail.lme <- lme(travel ~1, data=Rail, random= ~ 1|Rail)
summary(fm1Rail.lme)
fm1Rail.lmeML <- update(fm1Rail.lme, method="ML")
summary(fm1Rail.lmeML)
plot(fm1Rail.lme)
intervals(fm1Rail.lme)
anova(fm1Rail.lme)

# Section 1.2
data(ergoStool)
ergoStool
plot(ergoStool)
plot.design(ergoStool)
fm2Stool <- lme(effort ~ Type, data=ergoStool, random=~1|Subject)
summary(fm2Stool)
anova(fm2Stool)
fm3Stool <- update(fm2Stool, effort~Type-1)
summary(fm3Stool)
plot(fm2Stool)
plot(fm2Stool, form=resid(., type="p") ~ fitted(.)|Subject, abline=0)

# Section 1.3
data(Machines)
plot(Machines)
attach(Machines)
interaction.plot(Machine, Worker, score)
detach(Machines)
fm1Machine <- lme(score~Machine, data=Machines, random=~1|Worker)
summary(fm1Machine)
fm2Machine <- update(fm1Machine, random=~1|Worker/Machine)
summary(fm2Machine)
anova(fm1Machine, fm2Machine)
intervals(fm2Machine)
fm3Machine <- update(fm1Machine, random=~Machine-1|Worker)
summary(fm3Machine)
anova(fm1Machine, fm2Machine, fm3Machine)

# Section 1.4
data(Orthodont)
plot(Orthodont)
OrthoFem <- Orthodont[Orthodont$Sex == "Female",]
getGroupsFormula(OrthoFem)
fm1OrthF.lis <- lmList(distance ~ age, data=OrthoFem)
coef(fm1OrthF.lis)
plot(intervals (fm1OrthF.lis))
fm2OrthF.lis <- update(fm1OrthF.lis, distance ~ I(age-11))
plot(intervals (fm2OrthF.lis))
fm1OrthF <- lme(distance ~ age, data=OrthoFem, random=~1|Subject)
summary(fm1OrthF)
fm1OrthFM <- update(fm1OrthF, method="ML")
fm2OrthF <- update(fm1OrthF, random= ~ age|Subject)
anova(fm1OrthF, fm2OrthF)
random.effects(fm1OrthF)
plot(compareFits(coef(fm1OrthF), coef(fm1OrthFM)))
plot(augPred(fm1OrthF), aspect="xy", grid=T)

# Section 1.5
data(Pixel)
plot(Pixel)
fm1Pixel <- lme(pixel ~ day + I(day^2), data=Pixel, random=list(Dog=~day, Side=~1))
fm2Pixel <- update(fm1Pixel, random=list(Dog=~day))
anova(fm1Pixel, fm2Pixel)
plot(fm1Pixel, form=resid(., type="p") ~ fitted(.)|Dog, abline=0)

# Section 1.6
data(Oats)
plot(Oats)
fm1Oats <-lme(yield ~ ordered(nitro)*Variety, data=Oats, random=~1|Block/Variety)
anova(fm1Oats)
fm2Oats <- update(fm1Oats, yield~ ordered(nitro) + Variety)
summary(fm2Oats)
fm3Oats <- update(fm1Oats, yield~ ordered(nitro))
summary(fm3Oats)
fm4Oats <- update(fm1Oats, yield~ nitro)
plot(augPred(fm4Oats), aspect=2.5, layout=c(6,3), between=list(x=c(0,0,0.5)))