library(nlme)

# Section 5.2
# Dialyzer data
data(Dialyzer)
plot(Dialyzer)
options(contrasts=c("contr.helmert","contr.poly")) # to get same values as in PB
fm1Dial.lme <- lme(rate~(pressure+I(pressure^2)+I(pressure^3)+I(pressure^4))*QB, data=Dialyzer, random=~pressure+I(pressure^2))
fm1Dial.lme
plot(fm1Dial.lme,resid(.) ~ pressure, abline=0)
plot(fm1Dial.lme,resid(.,type="p") ~ pressure, abline=0)
fm2Dial.lme <- update(fm1Dial.lme, weights=varPower(form=~pressure))
fm2Dial.lme
anova(fm1Dial.lme,fm2Dial.lme)
plot(fm2Dial.lme,resid(.) ~ pressure, abline=0)
plot(fm2Dial.lme,resid(.,type="p") ~ pressure, abline=0)
plot(fm2Dial.lme,resid(.) ~ pressure|QB, abline=0)
plot(fm2Dial.lme,resid(.,type="p") ~ pressure|QB, abline=0)
fm3Dial.lme <- update(fm2Dial.lme, weights=varPower(form=~pressure|QB))
fm3Dial.lme
anova(fm2Dial.lme,fm3Dial.lme)
fm4Dial.lme <- update(fm2Dial.lme, weights=varConstPower(form=~pressure))
fm4Dial.lme
anova(fm2Dial.lme,fm4Dial.lme)
plot(augPred(fm2Dial.lme),grid=T)
anova(fm2Dial.lme)
anova(fm2Dial.lme,Terms=8:10)

# BodyWeight data
data(BodyWeight)
options(contrasts=c("contr.treatment","contr.poly"))
plot(BodyWeight)
fm1BW.lme <- lme(weight~Time*Diet, data=BodyWeight, random=~Time)
fm1BW.lme
plot(fm1BW.lme, resid(.,type="p")~fitted(.), abline=0)
fm2BW.lme <- update(fm1BW.lme, weights=varPower())
plot(fm2BW.lme, resid(.,type="p")~fitted(.), abline=0)
anova(fm1BW.lme,fm2BW.lme)
summary(fm2BW.lme)
anova(fm2BW.lme, L=c("Time:Diet2"=1,"Time:Diet3"=-1))

# Section 5.3
# Ovary data
data(Ovary)
plot(Ovary)
fm1Ovar.lme <- lme(follicles~I(sin(2*pi*Time)) + I(cos(2*pi*Time)), data=Ovary, random=pdDiag(~I(sin(2*pi*Time))))
fm1Ovar.lme
ACF(fm1Ovar.lme)
plot(ACF(fm1Ovar.lme,maxLag=10),alpha=0.01)
fm2Ovar.lme <- update(fm1Ovar.lme,correlation=corAR1())
fm2Ovar.lme
fm2Ovar.lme <- update(fm1Ovar.lme,correlation=corAR1(value=c(0.6))) # values in PB
fm2Ovar.lme
anova(fm1Ovar.lme,fm2Ovar.lme)
plot(ACF(fm2Ovar.lme,resType="n",maxLag=10),alpha=0.01)
fm3Ovar.lme <- update(fm1Ovar.lme,correlation=corARMA(q=2))
fm3Ovar.lme
anova(fm2Ovar.lme,fm3Ovar.lme,test=FALSE)
plot(ACF(fm3Ovar.lme,resType="n",maxLag=10),alpha=0.01)
fm4Ovar.lme <- update(fm1Ovar.lme,correlation=corCAR1(form=~Time))
fm4Ovar.lme
fm4Ovar.lme <- update(fm1Ovar.lme,correlation=corCAR1(value=c(0.002),form=~Time)) # values in PB
fm4Ovar.lme
anova(fm2Ovar.lme,fm4Ovar.lme)
plot(ACF(fm4Ovar.lme,resType="n",maxLag=10),alpha=0.01)
fm5Ovar.lme <- update(fm1Ovar.lme,correlation=corARMA(value=c(0.6,0),p=1,q=1))
fm5Ovar.lme
anova(fm2Ovar.lme,fm5Ovar.lme)
plot(ACF(fm5Ovar.lme,resType="n",maxLag=10),alpha=0.01)

# BodyWeight data
data(BodyWeight)
fm2BW.lme <- lme(weight~Time*Diet, data=BodyWeight, random=~Time, weights=varPower())
Variogram(fm2BW.lme, form=~Time)
plot(Variogram(fm2BW.lme, form=~Time, maxDist=42))
fm3BW.lme <- update(fm2BW.lme, corr=corExp(form=~Time))
fm3BW.lme
anova(fm2BW.lme,fm3BW.lme)
plot(Variogram(fm3BW.lme,form=~Time,maxDist=42))
plot(Variogram(fm3BW.lme,form=~Time,maxDist=42,resType="n",robust=T))
# no convergence: fm4BW.lme <- update(fm2BW.lme, corr=corExp(form=~Time,nugget=TRUE))
fm5BW.lme <- update(fm2BW.lme, corr=corRatio(form=~Time))
fm6BW.lme <- update(fm2BW.lme, corr=corSpher(form=~Time))
fm7BW.lme <- update(fm2BW.lme, corr=corLin(form=~Time))
fm8BW.lme <- update(fm2BW.lme, corr=corGaus(form=~Time))
anova(fm3BW.lme,fm5BW.lme,fm6BW.lme,fm7BW.lme,fm8BW.lme)

# Section 5.4
# Orthodont data
data(Orthodont)
fm2Orth.lme <- lme(distance~Sex*I(age-11), random=~I(age-11)|Subject, data=Orthodont)
fm2Orth.lme
fm3Orth.lme <- update(fm2Orth.lme, weights=varIdent(form=~1|Sex))
fm3Orth.lme
fm1Orth.gls <- gls(distance~Sex*I(age-11), data=Orthodont, corr=corSymm(form=~1|Subject), weights=varIdent(form=~1|age))
fm1Orth.gls
intervals(fm1Orth.gls)
fm2Orth.gls <- update(fm1Orth.gls, corr=corCompSymm(form=~1|Subject))
fm2Orth.gls
anova(fm1Orth.gls,fm2Orth.gls)
fm3Orth.gls <- update(fm2Orth.gls, weights=NULL)
fm3Orth.gls
anova(fm2Orth.gls,fm3Orth.gls)
plot(fm3Orth.gls, resid(.,type="n")~age|Sex)
fm4Orth.gls <- update(fm2Orth.gls, weights=varIdent(form=~1|Sex))
fm4Orth.gls
anova(fm3Orth.gls,fm4Orth.gls)
# extra models
fm4Orth.lme <- update(fm3Orth.lme, random=~1|Subject)
fm4Orth.lme
boy <- as.numeric(Orthodont$Sex=="Male")
girl <- as.numeric(Orthodont$Sex=="Female")
fm5Orth.lme <- update(fm3Orth.lme, random=pdDiag(~boy+girl-1))
fm5Orth.lme
anova(fm4Orth.lme,fm5Orth.lme,fm4Orth.gls)

# Dialyzer data
data(Dialyzer)
options(contrasts=c("contr.helmert","contr.poly")) # to get same values as in PB
options(contrasts=c("contr.treatment","contr.poly"))
fm1Dial.gls <- gls(rate~(pressure+I(pressure^2)+I(pressure^3)+I(pressure^4))*QB, data=Dialyzer)
fm1Dial.gls
plot(fm1Dial.gls,resid(.) ~ pressure, abline=0)
fm2Dial.gls <- update(fm1Dial.gls, weights=varPower(form=~pressure))
fm2Dial.gls
anova(fm1Dial.gls,fm2Dial.gls)
plot(fm2Dial.gls,resid(.) ~ pressure, abline=0)
plot(fm2Dial.gls,resid(.,type="p") ~ pressure, abline=0)
ACF(fm2Dial.gls, form=~1|Subject)
plot(ACF(fm2Dial.gls, form=~1|Subject),alpha=0.01)
fm3Dial.gls <- update(fm2Dial.gls, weights=varPower(form=~pressure),corr=corAR1(0.771,form=~1|Subject))
fm3Dial.gls
intervals(fm3Dial.gls)
anova(fm2Dial.gls,fm3Dial.gls)
plot(ACF(fm3Dial.gls, form=~1|Subject,resType="n"),alpha=0.01)
anova(fm2Dial.gls,fm2Dial.lme,test=F)

# Wheat2 data
data(Wheat2)
plot(Wheat2)
fm1Wheat2 <- gls(yield~variety-1,Wheat2)
Variogram(fm1Wheat2,form=~latitude+longitude)
plot(Variogram(fm1Wheat2,form=~latitude+longitude,maxDist=32),xlim=c(0,32))
fm2Wheat2 <- update(fm1Wheat2, corr=corSpher(c(28,0.2),form=~latitude+longitude,nugget=TRUE))
fm2Wheat2
fm3Wheat2 <- update(fm1Wheat2, corr=corRatio(c(12.5,0.2),form=~latitude+longitude,nugget=TRUE))
fm3Wheat2
anova(fm1Wheat2,fm3Wheat2)
plot(Variogram(fm3Wheat2,resType="n"))
plot(fm3Wheat2, resid(.,type="n")~fitted(.),abline=0)
qqnorm(fm3Wheat2, ~ resid(.,type="n"))
fm4Wheat2 <- update(fm3Wheat2, model=yield~variety)
anova(fm4Wheat2)
anova(fm3Wheat2, L=c(-1,0,1))