% 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/interactions,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(Machines, package = "MEMSS")
#lattice.options(default.theme = standard.theme(color = FALSE))
@ 

\begin{frame}
  \frametitle{Interactions of covariates and grouping factors}
  \begin{itemize}
  \item For longitudinal data, having a random effect for the slope
    w.r.t. time by subject is reasonably easy to understand.
  \item Although not generally presented in this way, these random
    effects are an interaction term between the grouping factor for the
    random effect (subject) and the time covariate.
  \item We can also define interactions between discrete covariates in the
    fixed-effects terms and a random-effects grouping factor.  However,
    there is more than one way to define such an interaction.
  \item Different ways of expressing such interactions lead to
    different numbers of random effects.
  \item Models with interactions defined in different ways have levels
    of complexity, affecting both their expressive power and the
    ability to estimate all the parameters in the model.
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{Machines data}
  \begin{itemize}
  \item Milliken and Johnson (1989) provide (probably artificial) data
    on an experiment to measure productivity according to the machine
    being used for a particular operation.
  \item In the experiment, a sample of six different operators used
    each of the three machines on three occasions --- a total of nine
    runs per operator.
  \item These three machines were the specific machines of interest
    and we model their effect as a fixed-effect term.
  \item The operators represented a sample from the population of
    potential operators.  We model this factor, (\code{Worker}), as a
    random effect.
  \item This is a replicated ``subject/stimulus'' design with a fixed
    set of stimuli that are themselves of interest.  (In other
    situations the stimuli may be a sample from a population of
    stimuli.)
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Machines data plot}
<<Machinesplot,fig=TRUE,echo=FALSE>>=
print(dotplot(reorder(Worker, score) ~ score, Machines,
              groups = Machine, 
              xlab = "Quality and productivity score", 
              ylab = "Worker", type = c("p","a"), 
              auto.key = list(columns = 3, lines = TRUE), 
              jitter.y = TRUE))
@ 
\end{frame}

\begin{frame}
  \frametitle{Comments on the data plot}
  \begin{itemize}
  \item There are obvious differences between the scores on different machines.
  \item It seems likely that \code{Worker} will be a significant
    random effect, especially when considering the low variation
    within replicates.
  \item There also appears to be a significant \code{Worker:Machine}
    interaction.  \code{Worker} 6 has a very different pattern
    w.r.t. machines than do the others.
  \item We can approach the interaction in one of two ways: define
    simple, scalar random effects for \code{Worker} and for the
    \code{Worker:Machine} interaction or define vector-valued random
    effects for \code{Worker}
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Random effects for subject and subject/stimulus}
<<MachinesModel1>>=
print(fm1 <- lmer(score ~ Machine + (1|Worker) + (1|Worker:Machine),
                  Machines), corr = FALSE)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Vector-valued random effects by subject}
<<MachinesModel1>>=
print(fm2 <- lmer(score ~ Machine + (0+Machine|Worker), Machines),
      corr = FALSE)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Comparing the model fits}
  \begin{itemize}
  \item Although not obvious from the specifications, the model fits
    are nested.  If the variance-covariance matrix for the
    vector-valued random effects has a special form, called
    \Emph{compound symmetry}, the model reduces to model \code{fm1}.
  \item The p-value from this comparison is borderline significant.
  \end{itemize}
<<Machinesanova>>=
fm2M <- update(fm2, REML = FALSE)
fm1M <- update(fm1, REML = FALSE)
anova(fm2M, fm1M)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Model comparisons eliminating the unusual combination}
  \begin{itemize}
  \item In a case like this we may want to check if a single, unusual
    combination (\code{Worker} 6 using \code{Machine} ``B'') causes
    the more complex model to appear necessary.  We eliminate that
    unusual combination.
  \end{itemize}
<<MachineswoB6>>=  
Machines1 <- subset(Machines, Worker != "6" | Machine != 'B')
xtabs(~ Machine + Worker, Machines1)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Machines data after eliminating the unusual combination}
<<Machines1plot,fig=TRUE,echo=FALSE>>=
print(dotplot(reorder(Worker, score) ~ score, Machines1,
              groups = Machine, 
              xlab = "Quality and productivity score", 
              ylab = "Worker", type = c("p","a"), 
              auto.key = list(columns = 3, lines = TRUE), 
              jitter.y = TRUE))
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Model comparisons without the unusual combination}
<<fm1aM>>=  
fm1aM <- lmer(score ~ Machine + (1|Worker) + (1|Worker:Machine), Machines1, REML = FALSE)
fm2aM <- lmer(score ~ Machine + (0 + Machine|Worker), Machines1, REML = FALSE)
anova(fm2aM, fm1aM)
@   
\end{frame}
\begin{frame}
  \frametitle{Trade-offs when defining interactions}
  \begin{itemize}
  \item It is important to realize that estimating scale parameters (i.e.
    variances and covariances) is considerably more difficult than
    estimating location parameters (i.e. means or fixed-effects coefficients).
  \item A vector-valued random effect term having $q_i$ random effects
    per level of the grouping factor requires $q_i(q_i+1)/2$
    variance-covariance parameters to be estimated.  A simple, scalar
    random effect for the interaction of a ``random-effects'' factor
    and a ``fixed-effects'' factor requires only 1 additional
    variance-covariance parameter.
  \item Especially when the ``fixed-effects'' factor has a moderate to
    large number of levels, the trade-off in model complexity argues
    against the vector-valued approach.
  \item One of the major sources of difficulty in using the
    \code{lme4} package is the tendency to overspecify the number of
    random effects per level of a grouping factor.
  \end{itemize}
\end{frame}