Nothing
#-*- R -*-
# initialization
library(nlme)
library(lattice)
options(width = 65,
## reduce platform dependence in printed output when testing
digits = if(nzchar(Sys.getenv("R_TESTS"))) 3 else 5)
options(contrasts = c(unordered = "contr.helmert", ordered = "contr.poly"))
pdf(file = 'ch04.pdf')
# Chapter 4 Fitting Linear Mixed-Effects Models
# 4.1 Fitting Linear Models in S with lm and lmList
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
plot(compOrth, mark = fixef(fm1Orth.lme)) # Figure 4.8
## 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))))
## IGNORE_RDIFF_BEGIN
summary(fm1Wafer)
## IGNORE_RDIFF_END
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)
})
with(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), Wafer,
start = list(b3 = -0.0519, b4 = 0.1304, w = 4.19))
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)
## IGNORE_RDIFF_BEGIN
intervals(fm2Wafer)
## IGNORE_RDIFF_END
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
summary(warnings())
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.