###################################################
### chunk number 1: preliminaries
###################################################
options(width = 70, show.signif.stars = FALSE)
library(lme4)
library(coda)
set.seed(123454321)

###################################################
### Slide 16
###################################################
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)"))
                    
###################################################
### Slide 20
###################################################
print(plot(confint(lmList(Reaction ~ Days | Subject, sleepstudy),
                   pooled = TRUE), order = 1))
                   
###################################################
### Slide 26
###################################################
(sm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))

###################################################
### Slide 28
###################################################
fixef(sm1)
VarCorr(sm1)$Subject
VarCorr(sm1)$Subject@factors$correlation

###################################################
### Slide 32
###################################################
ss1 <- mcmcsamp(sm1, 50000, trans = FALSE)
print(densityplot(ss1, plot.points = FALSE))

###################################################
### slide 35
###################################################
print(qqmath(window(ss1,thin = 10), pch = '.',
             cex = 2.5, type = c("g", "p"),
             xlab = "Quantiles of a standard normal distribution"))
             
###################################################
### Slide 38
###################################################
set.seed(123454321)
ss1a <- mcmcsamp(sm1, 50000)
print(densityplot(ss1a, plot.points = FALSE))

###################################################
### Slide 40
###################################################
print(qqmath(window(ss1a,thin = 10), pch = '.',
             cex = 2.5, type = c("g", "p"),
             xlab = "Quantiles of a standard normal distribution"))
             
###################################################
### Slide 43
###################################################
HPDinterval(ss1)
HPDinterval(ss1a)

###################################################
### Slide 47
###################################################
(sm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy))

###################################################
### Slide 49
###################################################
anova(sm2, sm1)

###################################################
### Slide 53
###################################################
ss2 <- mcmcsamp(sm2, 10000)
head(HPDinterval(ss2), 2)
sc2 <- summary(sm2)@coefs
sc2[,1] + qt(0.975, 144) * outer(sc2[,2], c(-1,1))

###################################################
### Slide 56
###################################################
sm3 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
anova(sm3, sm2)

###################################################
### Slide 63
###################################################
(rr1 <- ranef(sm1))

###################################################
### chunk number 22: rr1plot
###################################################
print(plot(rr1, aspect = 1, type = c("g", "p")))

###################################################
### Slide 65
###################################################
df <- coef(lmList(Reaction ~ Days | Subject, sleepstudy))
cc1 <- as.data.frame(coef(sm1)$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]))
               )))
               
###################################################
### Slide 67
###################################################
print(xyplot(Reaction ~ Days | Subject, sleepstudy, aspect = "xy",
             layout = c(9,2), type = c("g", "p", "r"),
             coef.list = df[,3:4],
             panel = function(..., packet.number, 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 25: caterpillar
###################################################
print(Matrix:::qqmath(ranef(sm1,post=TRUE))$Subject)