#-*- R -*-

# initialization

library(nlme)
library(lattice)
options(width = 65, digits = 5)
options(contrasts = c(unordered = "contr.helmert", ordered = "contr.poly"))
postscript(file = 'ch04.ps')

# Chapter 4    Fitting Linear Mixed-Effects Models

# 4.1 Fitting Linear Models in S with lm and lmList

##data(Orthodont, IGF, Oats, Assay, Oxide, Wafer, Machines, package = "nlme")
fm1Orth.lm <- lm(distance ~ age, Orthodont)
fm1Orth.lm
par(mfrow=c(2,2))
plot(fm1Orth.lm)                               # Figure 4.1
fm2Orth.lm <- update(fm1Orth.lm, formula = distance ~ Sex*age)
summary(fm2Orth.lm)
fm3Orth.lm <- update(fm2Orth.lm, formula = . ~ . - Sex)
summary(fm3Orth.lm)
bwplot(getGroups(Orthodont)~resid(fm2Orth.lm)) # Figure 4.2
fm1Orth.lis <- lmList(distance ~ age | Subject, Orthodont)
getGroupsFormula(Orthodont)
fm1Orth.lis <- lmList(distance ~ age, Orthodont)
formula(Orthodont)
fm1Orth.lis <- lmList(Orthodont)
fm1Orth.lis
summary(fm1Orth.lis)
pairs(fm1Orth.lis, id = 0.01, adj = -0.5)      # Figure 4.3
fm2Orth.lis <- update(fm1Orth.lis, distance ~ I(age-11))
intervals(fm2Orth.lis)
plot(intervals(fm2Orth.lis))                   # Figure 4.5
IGF
plot(IGF)                                      # Figure 4.6
fm1IGF.lis <- lmList(IGF)
coef(fm1IGF.lis)
plot(intervals(fm1IGF.lis))                    # Figure 4.7
fm1IGF.lm <- lm(conc ~ age, data = IGF)
summary(fm1IGF.lm)

# 4.2 Fitting Linear Mixed-Effects Models with lme

fm1Orth.lme <- lme(distance ~ I(age-11), data = Orthodont,
                     random = ~ I(age-11) | Subject)
fm1Orth.lme <- lme(distance ~ I(age-11), data = Orthodont)
fm1Orth.lme <- lme(fm2Orth.lis)
fm1Orth.lme
fm2Orth.lme <- update(fm1Orth.lme, distance~Sex*I(age-11))
summary(fm2Orth.lme)
fitted(fm2Orth.lme, level = 0:1)
resid(fm2Orth.lme, level = 1)
resid(fm2Orth.lme, level = 1, type = "pearson")
newOrth <- data.frame(Subject = rep(c("M11","F03"), c(3, 3)),
                        Sex = rep(c("Male", "Female"), c(3, 3)),
                        age = rep(16:18, 2))
predict(fm2Orth.lme, newdata = newOrth)
predict(fm2Orth.lme, newdata = newOrth, level = 0:1)
fm2Orth.lmeM <- update(fm2Orth.lme, method = "ML")
summary(fm2Orth.lmeM)
compOrth <-
      compareFits(coef(fm2Orth.lis), coef(fm1Orth.lme))
compOrth
                                        # Figure 4.8
plot(compOrth, mark = fixef(fm1Orth.lme))
                                        # Figure 4.9
plot(comparePred(fm2Orth.lis, fm1Orth.lme, length.out = 2),
     layout = c(8,4), between = list(y = c(0, 0.5, 0)))
plot(compareFits(ranef(fm2Orth.lme), ranef(fm2Orth.lmeM)),
     mark = c(0, 0))
fm4Orth.lm <- lm(distance ~ Sex * I(age-11), Orthodont)
summary(fm4Orth.lm)
anova(fm2Orth.lme, fm4Orth.lm)
#fm1IGF.lme <- lme(fm1IGF.lis)
#fm1IGF.lme
#intervals(fm1IGF.lme)
#summary(fm1IGF.lme)
pd1 <- pdDiag(~ age)
pd1
formula(pd1)
#fm2IGF.lme <- update(fm1IGF.lme, random = pdDiag(~age))
(fm2IGF.lme <- lme(conc ~ age, IGF,
                   random = pdDiag(~age)))
#anova(fm1IGF.lme, fm2IGF.lme)
anova(fm2IGF.lme)
#update(fm1IGF.lme, random = list(Lot = pdDiag(~ age)))
pd2 <- pdDiag(value = diag(2), form = ~ age)
pd2
formula(pd2)
lme(conc ~ age, IGF, pdDiag(diag(2), ~age))
fm4OatsB <- lme(yield ~ nitro, data = Oats,
                 random =list(Block = pdCompSymm(~ Variety - 1)))
summary(fm4OatsB)
corMatrix(fm4OatsB$modelStruct$reStruct$Block)[1,2]
fm4OatsC <- lme(yield ~ nitro, data = Oats,
        random=list(Block=pdBlocked(list(pdIdent(~ 1),
                                         pdIdent(~ Variety-1)))))
summary(fm4OatsC)
## establishing the desired parameterization for contrasts
options(contrasts = c("contr.treatment", "contr.poly"))
fm1Assay <- lme(logDens ~ sample * dilut, Assay,
   random = pdBlocked(list(pdIdent(~ 1), pdIdent(~ sample - 1),
                      pdIdent(~ dilut - 1))))
fm1Assay
anova(fm1Assay)
formula(Oxide)
fm1Oxide <- lme(Thickness ~ 1, Oxide)
fm1Oxide
intervals(fm1Oxide, which = "var-cov")
fm2Oxide <- update(fm1Oxide, random = ~ 1 | Lot)
anova(fm1Oxide, fm2Oxide)
coef(fm1Oxide, level = 1)
coef(fm1Oxide, level = 2)
ranef(fm1Oxide, level = 1:2)
fm1Wafer <- lme(current ~ voltage + I(voltage^2), data = Wafer,
              random = list(Wafer = pdDiag(~voltage + I(voltage^2)),
                            Site = pdDiag(~voltage + I(voltage^2))))
summary(fm1Wafer)
fitted(fm1Wafer, level = 0)
resid(fm1Wafer, level = 1:2)
newWafer <-
    data.frame(Wafer = rep(1, 4), voltage = c(1, 1.5, 3, 3.5))
predict(fm1Wafer, newWafer, level = 0:1)
newWafer2 <- data.frame(Wafer = rep(1, 4), Site = rep(3, 4),
                          voltage = c(1, 1.5, 3, 3.5))
predict(fm1Wafer, newWafer2, level = 0:2)

# 4.3 Examining a Fitted Model

plot(fm2Orth.lme, Subject~resid(.), abline = 0)
plot(fm2Orth.lme, resid(., type = "p") ~ fitted(.) | Sex,
      id = 0.05, adj = -0.3)
fm3Orth.lme <-
  update(fm2Orth.lme, weights = varIdent(form = ~ 1 | Sex))
fm3Orth.lme
plot(fm3Orth.lme, distance ~ fitted(.),
      id = 0.05, adj = -0.3)
anova(fm2Orth.lme, fm3Orth.lme)
qqnorm(fm3Orth.lme, ~resid(.) | Sex)
plot(fm2IGF.lme, resid(., type = "p") ~ fitted(.) | Lot,
      layout = c(5,2))
qqnorm(fm2IGF.lme, ~ resid(.),
        id = 0.05, adj = -0.75)
plot(fm1Oxide)
qqnorm(fm1Oxide)
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)
              })
attach(Wafer)
coef(lm(resid(fm1Wafer) ~ cos(4.19*voltage)+sin(4.19*voltage)-1))
nls(resid(fm1Wafer) ~ b3*cos(w*voltage) + b4*sin(w*voltage),
      start = list(b3 = -0.0519, b4 = 0.1304, w = 4.19))
detach()
fm2Wafer <- update(fm1Wafer,
      . ~ . + cos(4.5679*voltage) + sin(4.5679*voltage),
      random = list(Wafer=pdDiag(~voltage+I(voltage^2)),
             Site=pdDiag(~voltage+I(voltage^2))))
summary(fm2Wafer)
intervals(fm2Wafer)
qqnorm(fm2Wafer)
qqnorm(fm2Orth.lme, ~ranef(.),
         id = 0.10, cex = 0.7)
pairs(fm2Orth.lme, ~ranef(.) | Sex,
        id = ~ Subject == "M13", adj = -0.3)
fm2IGF.lme
c(0.00031074, 0.0053722)/abs(fixef(fm2IGF.lme))
fm3IGF.lme <- update(fm2IGF.lme, random = ~ age - 1)
anova(fm2IGF.lme, fm3IGF.lme)
qqnorm(fm1Oxide, ~ranef(., level = 1), id=0.10)
qqnorm(fm1Oxide, ~ranef(., level = 2), id=0.10)
#fm3Wafer <- update(fm2Wafer,
#              random = list(Wafer = ~voltage+I(voltage^2),
#                            Site = pdDiag(~voltage+I(voltage^2))),
#                   control = list(msVerbose = TRUE, msMaxIter = 200)
#                   )
#fm3Wafer
#anova(fm2Wafer, fm3Wafer)
#fm4Wafer <- update(fm2Wafer,
#                   random = list(Wafer = ~ voltage + I(voltage^2),
#                   Site = pdBlocked(list(~1,
#                   ~voltage+I(voltage^2) - 1))),
#                   control = list(msVerbose = TRUE,
#                   msMaxIter = 200))
#fm4Wafer
#anova(fm3Wafer, fm4Wafer)
#qqnorm(fm4Wafer, ~ranef(., level = 2), id = 0.05,
#        cex = 0.7, layout = c(3, 1))

# The next line is not in the book but is needed to get fm1Machine

fm1Machine <-
  lme(score ~ Machine, data = Machines, random = ~ 1 | Worker)

(fm3Machine <- update(fm1Machine, random = ~Machine-1|Worker))

# cleanup

proc.time()
q("no")
