###################################################
### chunk number 1: preliminaries
###################################################
options(width=65)
library(lattice)
library(Matrix)
library(lme4)
data(ScotsSec, package = "mlmRev")
#lattice.options(default.theme = standard.theme(color = FALSE))


###################################################
### chunk number 2: Scotssec
###################################################
str(ScotsSec)
stab <- xtabs(~ second + primary, ScotsSec, sparse = TRUE)


###################################################
### chunk number 3: stabim
###################################################
print(image(stab, xlab = "Primary school", ylab = "Secondary",sub=NULL))


###################################################
### chunk number 4: mattpatt
###################################################
patt <- with(ScotsSec, data.frame(mattain = unname(tapply(attain, primary, mean)),
                                  n = as.vector(unname(table(primary))),
                                  type = "Primary"))
row.names(patt) <- paste(levels(ScotsSec$primary), "P", sep = "")
satt <- with(ScotsSec, data.frame(mattain = unname(tapply(attain, second, mean)),
                                  n = as.vector(unname(table(second))),
                                  type = "Secondary"))
row.names(satt) <- paste(levels(ScotsSec$second), "S", sep = "")


###################################################
### chunk number 5: mattain
###################################################
head(patt)
head(satt)


###################################################
### chunk number 6: attainQQ
###################################################
print(qqmath(~mattain, rbind(patt,satt), groups = type, aspect = 1, type = c("g", "p"),
             xlab = "Standard normal quantiles",
             ylab = "Mean attain score in school",
             auto.key = list(space = "right")))


###################################################
### chunk number 7: attainQQ2
###################################################
print(qqmath(~ verbal | sex, ScotsSec, groups = attain, aspect = 1.2, type = c("g", "p"),
             xlab = "Standard normal quantiles",
             ylab = "Verbal pretest score",
             auto.key = list(space = "right", title = "attain")))


###################################################
### chunk number 8: scotsLMM
###################################################
print(Sm1 <- lmer(attain ~ verbal + sex + (1|primary) + (verbal|second), ScotsSec), corr = FALSE)


###################################################
### chunk number 9: scotsLMM2
###################################################
print(Sm2 <- lmer(attain ~ verbal + sex + (1|primary) + (1|second), ScotsSec), corr = FALSE)


###################################################
### chunk number 10: reorder
###################################################
print(image(tcrossprod(Sm2@A), xlab = NULL, ylab = NULL, sub = "AA'"),
      position = c(0,0,0.5,1), more = TRUE)
print(image(Sm2@L, xlab = NULL, ylab = NULL, sub = "L", colorkey = FALSE),
      position = c(0.5,0,1,1))


###################################################
### chunk number 11: Sm2ranef
###################################################
print(dotplot(ranef(Sm2, post = TRUE), strip = FALSE)$second)


###################################################
### chunk number 12: Sm3
###################################################
Sm3 <- update(Sm2, . ~ . - (1|second))
anova(Sm3, Sm2)


###################################################
### chunk number 13: Sm3ranef
###################################################
print(dotplot(ranef(Sm3, post = TRUE), strip = FALSE)$primary)