The R Project SVN
R-packages
(root)/trunk/nlme/inst/scripts/ch04.R – Rev 8034
Rev 7855 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
#-*- R -*-# initializationlibrary(nlme)library(lattice)options(width = 65,## reduce platform dependence in printed output when testingdigits = 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 lmListfm1Orth.lm <- lm(distance ~ age, Orthodont)fm1Orth.lmpar(mfrow=c(2,2))plot(fm1Orth.lm) # Figure 4.1fm2Orth.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.2fm1Orth.lis <- lmList(distance ~ age | Subject, Orthodont)getGroupsFormula(Orthodont)fm1Orth.lis <- lmList(distance ~ age, Orthodont)formula(Orthodont)fm1Orth.lis <- lmList(Orthodont)fm1Orth.lissummary(fm1Orth.lis)pairs(fm1Orth.lis, id = 0.01, adj = -0.5) # Figure 4.3fm2Orth.lis <- update(fm1Orth.lis, distance ~ I(age-11))intervals(fm2Orth.lis)plot(intervals(fm2Orth.lis)) # Figure 4.5IGFplot(IGF) # Figure 4.6fm1IGF.lis <- lmList(IGF)coef(fm1IGF.lis)plot(intervals(fm1IGF.lis)) # Figure 4.7fm1IGF.lm <- lm(conc ~ age, data = IGF)summary(fm1IGF.lm)# 4.2 Fitting Linear Mixed-Effects Models with lmefm1Orth.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.lmefm2Orth.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))compOrthplot(compOrth, mark = fixef(fm1Orth.lme)) # Figure 4.8## Figure 4.9plot(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)pd1formula(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)pd2formula(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 contrastsoptions(contrasts = c("contr.treatment", "contr.poly"))fm1Assay <- lme(logDens ~ sample * dilut, Assay,random = pdBlocked(list(pdIdent(~ 1), pdIdent(~ sample - 1),pdIdent(~ dilut - 1))))fm1Assayanova(fm1Assay)formula(Oxide)fm1Oxide <- lme(Thickness ~ 1, Oxide)fm1Oxideintervals(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_BEGINsummary(fm1Wafer)## IGNORE_RDIFF_ENDfitted(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 Modelplot(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.lmeplot(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_BEGINintervals(fm2Wafer)## IGNORE_RDIFF_ENDqqnorm(fm2Wafer)qqnorm(fm2Orth.lme, ~ranef(.), id = 0.10, cex = 0.7)pairs(fm2Orth.lme, ~ranef(.) | Sex,id = ~ Subject == "M13", adj = -0.3)fm2IGF.lmec(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 fm1Machinefm1Machine <-lme(score ~ Machine, data = Machines, random = ~ 1 | Worker)(fm3Machine <- update(fm1Machine, random = ~Machine-1|Worker))# cleanupsummary(warnings())