###################################################
### 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(dotplot(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 <- dotplot(ranef(fm2, postVar = TRUE), strip = FALSE)
print(qrr2[[1]], pos = c(0,0,1,0.75), more = TRUE)
print(qrr2[[2]], pos = c(0,0.65,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 <- dotplot(ranef(fm3, postVar = TRUE), strip = FALSE)
print(qrr3[[1]], pos = c(0,0,1,0.75), more = TRUE)
print(qrr3[[2]], pos = c(0,0.65,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))