###################################################
### chunk number 1: preliminaries
###################################################
options(width=65,show.signif.stars=FALSE)
library(lattice)
library(Matrix)
library(lme4)
#lattice.options(default.theme = standard.theme(color = FALSE))


###################################################
### chunk number 2: sleepxy
###################################################
print(xyplot(Reaction ~ Days | Subject, sleepstudy, aspect = "xy",
                    layout = c(9,2), type = c("g", "p", "r"),
                    index.cond = function(x,y) coef(lm(y ~ x))[1],
                    xlab = "Days of sleep deprivation",
                    ylab = "Average reaction time (ms)"))


###################################################
### chunk number 3: Sl1
###################################################
print(plot(confint(lmList(Reaction ~ Days | Subject, sleepstudy),
                   pooled = TRUE), order = 1))


###################################################
### chunk number 4: fm1
###################################################
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)


###################################################
### chunk number 5: sleepZ
###################################################
print(image(fm1@Zt,xlab=NULL,ylab=NULL,sub=NULL))


###################################################
### chunk number 6: sm1
###################################################
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))


###################################################
### chunk number 7: fm1Images
###################################################
efm1 <- expand(fm1)
print(image(efm1[["T"]], xlab = NULL, ylab = NULL, sub = "T"),
      position = c(0,0,0.27,1), more = TRUE)
print(image(tcrossprod(fm1@A), xlab = NULL, ylab = NULL, sub = "AA'"),
      position = c(0.24,0,0.51,1), more = TRUE)
print(image(efm1$P, xlab = NULL, ylab = NULL, sub = "P"),
      position = c(0.49,0,0.76,1), more = TRUE)
print(image(fm1@L, xlab = NULL, ylab = NULL, sub = "L", colorkey = FALSE),
      position = c(0.73,0,1,1))


###################################################
### chunk number 8: fm2
###################################################
(fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy))


###################################################
### chunk number 9: anovafm1fm2
###################################################
anova(fm2, fm1)


###################################################
### chunk number 10: fm3anova
###################################################
fm3 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
anova(fm3, fm2)


###################################################
### chunk number 11: rr1
###################################################
(rr2 <- ranef(fm2))


###################################################
### chunk number 12: rr2plot
###################################################
print(plot(rr2, aspect = 1, type = c("g", "p"))[[1]])


###################################################
### chunk number 13: shrinkage
###################################################
df <- coef(lmList(Reaction ~ Days | Subject, sleepstudy))
cc1 <- as.data.frame(coef(fm2)[["Subject"]])
names(cc1) <- c("A", "B")
df <- cbind(df, cc1)
with(df,
     print(xyplot(`(Intercept)` ~ Days, aspect = 1,
            x1 = B, y1 = A, 
            panel = function(x, y, x1, y1, subscripts, ...) {
                panel.grid(h = -1, v = -1)
                x1 <- x1[subscripts]
                y1 <- y1[subscripts]
                panel.arrows(x, y, x1, y1, type = "closed", length = 0.1,
                             angle = 15, ...)
                panel.points(x, y,
                             col = trellis.par.get("superpose.symbol")$col[2])
                panel.points(x1, y1)
            },
                   key = list(space = "top", columns = 2,
                   text = list(c("Mixed model", "Within-group")),
                   points = list(col = trellis.par.get("superpose.symbol")$col[1:2],
                   pch = trellis.par.get("superpose.symbol")$pch[1:2]))
               )))


###################################################
### chunk number 14: shrinkfit
###################################################
print(xyplot(Reaction ~ Days | Subject, sleepstudy, aspect = "xy",
             layout = c(9,2), type = c("g", "p", "r"),
             coef.list = df[,3:4],
             panel = function(..., coef.list) {
                 panel.xyplot(...)
                 panel.abline(as.numeric(coef.list[packet.number(),]),
                              col.line = trellis.par.get("superpose.line")$col[2],
                              lty = trellis.par.get("superpose.line")$lty[2]
                              )
             },
             index.cond = function(x,y) coef(lm(y ~ x))[1],
             xlab = "Days of sleep deprivation",
             ylab = "Average reaction time (ms)",
             key = list(space = "top", columns = 2,
             text = list(c("Mixed model", "Within-group")),
             lines = list(col = trellis.par.get("superpose.line")$col[2:1],
             lty = trellis.par.get("superpose.line")$lty[2:1]))))


###################################################
### chunk number 15: caterpillar
###################################################
print(dotplot(ranef(fm1,post=TRUE),
              scales = list(x = list(relation = 'free')))[["Subject"]])