################################################### ### chunk number 1: preliminaries ################################################### options(width=65) library(lattice) library(Matrix) library(lme4) lattice.options(default.theme = standard.theme(color = FALSE)) ################################################### ### chunk number 2: installlme4 eval=FALSE ################################################### ## install.packages("lme4") ################################################### ### chunk number 3: require eval=FALSE ################################################### ## require(lme4) ################################################### ### chunk number 4: attach eval=FALSE ################################################### ## library(lme4) ################################################### ### chunk number 5: Dyestuffstr ################################################### str(Dyestuff) summary(Dyestuff) ################################################### ### chunk number 6: Dyestuffplot ################################################### print(dotplot(reorder(Batch, Yield) ~ Yield, Dyestuff, ylab = "Batch", jitter.y = TRUE, xlab = "Yield of dyestuff (grams of standard color)", type = c("p", "a"))) ################################################### ### chunk number 7: fm1 ################################################### fm1 <- lmer(Yield ~ 1 + (1|Batch), Dyestuff) print(fm1) ################################################### ### chunk number 8: op ################################################### op <- options(digits=5) ################################################### ### chunk number 9: extractors ################################################### fixef(fm1) ranef(fm1, drop = TRUE) fitted(fm1) ################################################### ### chunk number 10: unop ################################################### options(op) ################################################### ### chunk number 11: fm1ranef ################################################### print(qqmath(ranef(fm1, postVar = TRUE), strip = FALSE)[[1]]) ################################################### ### chunk number 12: Batchmm ################################################### with(Dyestuff, as(Batch, "sparseMatrix")) ################################################### ### chunk number 13: fm1refit ################################################### invisible(update(fm1, verbose = TRUE)) ################################################### ### chunk number 14: op1 ################################################### op <- options(digits = 5) ################################################### ### chunk number 15: fm1TS ################################################### efm1 <- expand(fm1) efm1$S efm1$T ################################################### ### chunk number 16: reconstructing ################################################### (fm1S <- tcrossprod(efm1$T %*% efm1$S)) ################################################### ### chunk number 17: fm1Images ################################################### print(image(efm1[["T"]], xlab = NULL, ylab = NULL, sub = "T"), position = c(0,0,0.36,1), more = TRUE) print(image(efm1[["S"]], xlab = NULL, ylab = NULL, sub = "S"), position = c(0.32,0,0.68,1), more = TRUE) print(image(fm1S, xlab = NULL, ylab = NULL, sub = expression(Sigma)), position = c(0.64,0,1,1)) ################################################### ### chunk number 18: unop1 ################################################### options(op) ################################################### ### chunk number 19: update ################################################### (fm1M <- update(fm1, REML = FALSE)) ################################################### ### chunk number 20: fm1call ################################################### fm1@call ################################################### ### chunk number 21: Penicillinstr ################################################### str(Penicillin) xtabs(~ sample + plate, Penicillin) ################################################### ### chunk number 22: PenicillinPlot ################################################### print(dotplot(reorder(plate, diameter) ~ diameter, Penicillin, groups = sample, ylab = "Plate", xlab = "Diameter of growth inhibition zone (mm)", type = c("p", "a"), auto.key = list(columns = 6, lines = TRUE))) ################################################### ### chunk number 23: fm2 ################################################### (fm2 <- lmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin)) ################################################### ### chunk number 24: ################################################### op <- options(digits = 5) ################################################### ### chunk number 25: fixef2 ################################################### fixef(fm2) ranef(fm2, drop = TRUE) ################################################### ### chunk number 26: ################################################### options(op) ################################################### ### chunk number 27: fm2ranef ################################################### qrr2 <- qqmath(ranef(fm2, postVar = TRUE), strip = FALSE) print(qrr2[[1]], pos = c(0,0,0.5,1), more = TRUE) print(qrr2[[2]], pos = c(0.5,0,1,1)) ################################################### ### chunk number 28: fm2Z ################################################### print(image(fm2@Zt, xlab = NULL, ylab = NULL, sub = "Z'")) ################################################### ### chunk number 29: fm2images ################################################### print(image(expand(fm2)$S, xlab = NULL, ylab = NULL, sub = "S"), position = c(0,0,0.36,1), more = TRUE) print(image(tcrossprod(fm2@A), xlab = NULL, ylab = NULL, sub = "AA'"), position = c(0.32,0,0.68,1), more = TRUE) print(image(fm2@L, xlab = NULL, ylab = NULL, sub = "L", colorkey = FALSE), position = c(0.64,0,1,1)) ################################################### ### chunk number 30: fm2formula ################################################### fm2@call[["formula"]] ################################################### ### chunk number 31: Pastesstr ################################################### str(Pastes) xtabs(~ batch + sample, Pastes, sparse = TRUE) ################################################### ### chunk number 32: samplegen eval=FALSE ################################################### ## Pastes <- within(Pastes, sample <- (batch:cask)[drop = TRUE]) ################################################### ### chunk number 33: Pastesplot ################################################### Pastes <- within(Pastes, bb <- reorder(batch, strength)) Pastes <- within(Pastes, ss <- reorder(reorder(sample, strength), as.numeric(batch))) print(dotplot(ss ~ strength | bb, Pastes, strip = FALSE, strip.left = TRUE, layout = c(1, 10), scales = list(y = list(relation = "free")), ylab = "Sample within batch", type = c("p", "a"), xlab = "Paste strength", jitter.y = TRUE)) ################################################### ### chunk number 34: fm3 ################################################### (fm3 <- lmer(strength ~ 1 + (1|batch) + (1|sample), Pastes)) ################################################### ### chunk number 35: fm3ranef ################################################### qrr3 <- qqmath(ranef(fm3, postVar = TRUE), strip = FALSE) print(qrr3[[1]], pos = c(0,0,0.5,1), more = TRUE) print(qrr3[[2]], pos = c(0.5,0,1,1)) ################################################### ### chunk number 36: Z3fig ################################################### print(image(fm3@Zt, xlab = NULL, ylab = NULL, sub = NULL)) ################################################### ### chunk number 37: fm3images ################################################### print(image(tcrossprod(fm3@A), xlab = NULL, ylab = NULL, sub = "AA'"), position = c(0,0,0.36,1), more = TRUE) print(image(expand(fm3)$P, xlab = NULL, ylab = NULL, sub = "P"), position = c(0.32,0,0.68,1), more = TRUE) print(image(fm3@L, xlab = NULL, ylab = NULL, sub = "L", colorkey = FALSE), position = c(0.64,0,1,1)) ################################################### ### chunk number 38: fm3LRT ################################################### fm3M <- update(fm3, REML = FALSE) fm4M <- lmer(strength ~ 1 + (1|sample), Pastes, REML = FALSE) anova(fm4M, fm3M) ################################################### ### chunk number 39: fm4 ################################################### (fm4 <- update(fm4M, REML = TRUE))