\SweaveOpts{engine=R,eps=FALSE,pdf=TRUE,width=10,height=6.5,strip.white=all}
\SweaveOpts{prefix=TRUE,prefix.string=figs/large,include=TRUE}
\setkeys{Gin}{width=\textwidth}
<<preliminaries,echo=FALSE,results=hide>>=
options(width=65)
library(lattice)
library(Matrix)
library(lme4)
data(ScotsSec, package = "mlmRev")
#lattice.options(default.theme = standard.theme(color = FALSE))
@ 

\begin{frame}[fragile]\frametitle{A smaller, non-nested unbalanced
    example}
  
  To examine the structure of non-nested, unbalanced observational
  data we need a smaller example than the 1.68 million observations in
  the 10 years of grade point scores.

<<Scotssec>>=  
str(ScotsSec)
stab <- xtabs(~ second + primary, ScotsSec, sparse = TRUE)
@ 
$\mbox{}$
<<stabim,fig=TRUE,echo=FALSE,height=2.4,results=hide>>=
print(image(stab, xlab = "Primary school", ylab = "Secondary",sub=NULL))
@ 
\end{frame}

<<mattpatt,echo=FALSE>>=
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 = "")
@ 

\begin{frame}[fragile]\frametitle{Mean attainment by school}
<<mattain>>=  
head(patt)
head(satt)
@ 
\end{frame}

\begin{frame}[fragile]\frametitle{Normal probability plot of mean attainment by school}
<<attainQQ,fig=TRUE,echo=FALSE, height = 7>>=  
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")))
@ 
\end{frame}

\begin{frame}[fragile]\frametitle{Probability plot of pretest by posttest and sex}
<<attainQQ2,fig=TRUE,echo=FALSE, height = 7>>=  
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")))
@ 
\end{frame}

\begin{frame}[fragile]\frametitle{An LMM for the secondary school data}
<<scotsLMM,echo=FALSE>>=  
print(Sm1 <- lmer(attain ~ verbal + sex + (1|primary) + (verbal|second), ScotsSec), corr = FALSE)
@ 
$n=3435$, $p=3$, $k=2$, $n_1=148$, $n_2=19$, $q_1=1$, $q_2=2$, $q=186$
\end{frame}

\begin{frame}[fragile]\frametitle{Reduced LMM for the secondary school data}
<<scotsLMM2,echo=FALSE>>=  
print(Sm2 <- lmer(attain ~ verbal + sex + (1|primary) + (1|second), ScotsSec), corr = FALSE)
@ 
$n=3435$, $p=3$, $k=2$, $n_1=148$, $n_2=19$, $q_1=1$, $q_2=1$, $q=167$
\end{frame}

\begin{frame}[fragile]
  \frametitle{Reordering unbalanced, non-nested data}
  \begin{center}
<<reorder,fig=TRUE,echo=FALSE>>=  
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))
@ 
  \end{center}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Secondary school random effects are poorly defined}
  \begin{center}
<<Sm2ranef,fig=TRUE,echo=FALSE,height=5>>=    
print(dotplot(ranef(Sm2, post = TRUE), strip = FALSE)$second)
@ 
  \end{center}
These indicate we should re-fit the model.
\end{frame}

\begin{frame}[fragile]
  \frametitle{Model without random effects for second}
<<Sm3>>=  
Sm3 <- update(Sm2, . ~ . - (1|second))
anova(Sm3, Sm2)
@ 
\end{frame}

\begin{frame}[fragile]
\frametitle{Random effects for primary are not that remarkable}
\begin{center}
<<Sm3ranef,fig=TRUE,echo=FALSE,height=7.2>>=  
print(dotplot(ranef(Sm3, post = TRUE), strip = FALSE)$primary)
@ 
\end{center}
\end{frame}