% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.

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

\begin{frame}
  \frametitle{Random-effects variances estimated as zero?}
  \begin{itemize}
  \item The estimated variance of the random effects for a given term
    must be non-negative. It can be shown that they must be also be
    finite.  However, it is possible for the estimate to be exactly
    zero, indicating an overspecified model.  In these cases we should
    respecify the model to have fewer random effects.
  \item Because the estimates of these variances are printed
    explicitly, it is obvious to us when such a situation occurs
    (although you need to watch for estimated variances that are very
    close to but not exactly zero).
  \item Box and Tiao in ``Bayesian Inference in Statistical Analysis''
    (Addison-Wesley, 1973) provide simulated data similar to the
    \code{Dyestuff} data but with much lower batch-to-batch variability.
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Dyestuff2 data plot}
  \begin{center}
<<Dyestuff2plot,fig=TRUE,echo=FALSE>>=    
print(dotplot(reorder(Batch, Yield) ~ Yield, Dyestuff2,
              ylab = "Batch", jitter.y = TRUE,
              xlab = "Yield of dyestuff (grams of standard color)",
              type = c("p", "a")))
@ 
  \end{center}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Model fit for Dyestuff2}
<<Dyestuff2mod>>=
lmer(Yield ~ 1 + (1|Batch), Dyestuff2)
@ 
\end{frame}

\begin{frame}
  \frametitle{Singular variance-covariance matrices}
  \begin{itemize}
  \item With vector-valued random effects we can obtained a singular
    or near-singular estimate of a variance-covariance matrix that has
    safely non-zero variances for all components.  It is harder to
    spot such singularities.  If there are only two random effects per
    level of the grouping factor we can recognize this situation as
    generating highly correlated (sometimes perfectly correlated)
    random effects.  For dimensions greater than 2, even this signal
    may fail to be present.
  \item I will likely add a diagnostic method with a name like
    \code{rcond} (reciprocal of the condition number) to \code{lme4}
    to make this situation easier to spot.
  \end{itemize}
\end{frame}
\begin{frame}
  \frametitle{The early childhood cognitive study data}
  \begin{itemize}
  \item In an experiment reported in Burchinal et al. (\emph{Cognitive
      Development, 1997}) young children were randomly assigned to a
    treatment group and a control group at 6 months of age and their
    cognitive development was measured at 1, 1.5 and 2 years of
    age, on an age-normed scale. (The treatment was exposure to an
    enriched environment.)
  \item Because the treatment began at 6 months of age, we will use
    \code{tos = age - 0.5} (time on study), as the time variable.
    Because the treatment groups are assumed to be homogeneous at the
    beginning of the study, we do not expect to see a \code{trt} main
    effect.  The effect of the treatment, if any, will show up in the
    \code{tos:trt} interaction term.
  \item These data are the initial example in Singer and Willett,
    \emph{Longitudinal Data Analysis} (Oxford, 2003).  Models are fit
    using both SAS and MLWin but the singularity in the estimated
    variance-covariance matrix is never noticed.
  \end{itemize}
\end{frame}
\begin{frame}[fragile]
  \frametitle{Early childhood cognitive data plot}
<<EarlyData,fig=TRUE,echo=FALSE>>=
print(xyplot(cog ~ age | id, Early, type = c('g', 'b'), aspect = 'xy',
             layout = c(29,4), between = list(y = c(0,0.5)),
             xlab = "Age (yr)",
             ylab = "Cognitive development score",
             scales = list(x = list(tick.number = 3, alternating = TRUE,
                           labels = c("1","","2"), at = c(1,1.5,2))),
             par.strip.text = list(cex = 0.6)))
@ 
\end{frame}
\begin{frame}[fragile]
  \frametitle{Some comments on the data plot}
  \begin{itemize}
  \item Notice that, for both groups, the slopes are almost always
    negative.  This is curious if the scale has been age-normed.
  \item Also, the initial measurement at age 1 yr seems high.  If the
    scale is age-normed then we would expect the observations in (at
    least) the control group to have a mean near 100.
  \end{itemize}
<<ttest>>=  
with(subset(Early, age == 1 & trt == 'N'), summary(cog))
t.test(subset(Early, age == 1 & trt == 'N')$cog, mu = 100)$p.value
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Comparative box-and-whisker plots}
<<EarlyBoxW,fig=TRUE,echo=FALSE>>=
print(bwplot(factor(age) ~ cog | trt, Early,
             strip = FALSE, strip.left = TRUE, layout = c(1, 2),
             ylab = "Age (yr) within group",
             xlab = "Cognitive development score"),
      position = c(0, 0, 1, 0.55), more = TRUE)
print(bwplot(trt ~ cog | factor(age), Early,
             strip = FALSE, strip.left = TRUE, layout = c(1, 3),
             ylab = "Group within Age (yr)", xlab = NULL),
      position = c(0, 0.5, 1, 1), more = FALSE)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Density by age within group}
<<EarlyDensity,fig=TRUE,echo=FALSE>>=
print(densityplot(~ cog|trt, Early, plot.points = FALSE,
                  ref = TRUE, layout = c(1,2),
                  xlab = "Cognitive development score",
                  groups = age, auto.key = list(columns = 3),
                  strip = FALSE, strip.left=TRUE))
@ 
\end{frame}


\begin{frame}[fragile]
  \frametitle{Density by group within age}
<<EarlyDensity2,fig=TRUE,echo=FALSE>>=
print(densityplot(~ cog|factor(age), Early, plot.points = FALSE,
                  ref = TRUE, layout = c(1,3),
                  xlab = "Cognitive development score",
                  groups = trt, auto.key = list(columns = 2),
                  strip = FALSE, strip.left=TRUE))
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Initial model fit}
<<Earlymodel>>=  
Early <- within(Early, tos <- age - 0.5)
fm1 <- lmer(cog ~ tos*trt + (tos|id), Early, verbose = TRUE)
@ 
The second parameter being exactly zero is a danger sign.
\end{frame}

\begin{frame}[fragile]
  \frametitle{Summary of initial model}
<<Earlymodelsummary>>=  
print(fm1, corr = FALSE)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Comments on the fitted model}
  \begin{itemize}
  \item As expected, the main effect for treatment group is not
    significant.  Recall that we changed the time scale so zero
    corresponds to the beginning of the intervention.
  \item The only danger sign in the output about singularity of the
    variance-covariance of the random effects is the fact that the
    estimated correlation is exactly -1.
  \item Even after refitting without the main effect for treatment,
    this high negative correlation is present.
  \item If you display the covariance instead of the correlation, it
    is very difficult to see that something is amiss.
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Re-fit without the main effect for treatment}
<<EarlyModel2>>=
print(fm2 <- update(fm1, . ~ . - trt, verbose = FALSE), corr = FALSE)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Conditional modes of the random effects}
<<EarlyCondMode,fig=TRUE,echo=FALSE>>=  
print(plot(ranef(fm2), type = c("g", "p"), aspect = 1)$id)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Handling singularity}
  \begin{itemize}
  \item We can try to redefine the model with fewer random effects or
    with uncorrelated random effects, if such a model would make sense
    and if the model fits are adequate.
  \item In this case, the reduced models with independent random
    effects or with random effects for the intercept only are not adequate.
  \item There is little that can be done here to handle the
    singularity in the model fit, because the singularity is related to
    anomalies in the data.
  \end{itemize}
<<Earlyup>>=
fm2M <- update(fm2, REML = FALSE)
fm3M <- lmer(cog ~ tos + tos:trt + (1|id) + (0+tos|id), Early, REML = FALSE)
fm4M <- lmer(cog ~ tos + tos:trt + (1|id), Early, REML = FALSE)
@ 
\end{frame}
\begin{frame}[fragile]
  \frametitle{Checking possible reduced models}
<<Earlyup>>=
anova(fm3M, fm2M)
anova(fm4M, fm2M)
@ 
\end{frame}