% 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/GLMM,include=TRUE}
\setkeys{Gin}{width=\textwidth}
<<preliminaries,echo=FALSE,results=hide>>=
options(width = 70, show.signif.stars = FALSE)
data(Contraception, package = "mlmRev")
library(lattice)
library(Matrix)
library(lme4)
#lattice.options(default.theme = standard.theme(color = FALSE))
@ 

\begin{frame}[fragile]
  \frametitle{The Contraception data set}
  \begin{itemize}
  \item One of the data sets in the \code{"mlmRev"} package, derived
    from data files available on the multilevel modelling web site, is
    from a fertility survey of women in Bangladesh.
  \item We consider a binary response - whether the woman currently
    uses artificial contraception.
  \item Covariates included the woman's age (on a centered scale),
    the number of live children she has, whether she lives in an urban
    or rural setting, and the district in which she lives.
  \item These data are quite unbalanced with regard to the covariates
    (some districts have only $2$ observations, some have nearly
    $120$).
  \item We should bear in mind that the binary responses have low
    per-observation information content (exactly one bit per
    observation).  Districts with few observations will not contribute
    strongly to estimates of random effects.
  \item Within-district plots will be too rough.  We can examine the
    influence of some of the covariates by plotting scatterplot
    smoother curves of the response versus age by other covariates.
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Contraception use versus age by urban and livch}
<<Contra1,fig=TRUE,echo=FALSE>>=
print(xyplot(ifelse(use == "Y", 1, 0) ~ age|urban, Contraception,
             groups = livch, type = c("g", "smooth"),
             auto.key = list(space = "top", points = FALSE,
             lines = TRUE, columns = 4),
             ylab = "Proportion", xlab = "Centered age"))
@ 
\end{frame}

\begin{frame}
  \frametitle{Comments on the data plot}
  \begin{itemize}
  \item On the multilevel modelling web site they compare various
    forms of multilevel software fitting a model that is linear in
    \code{age} to these data.  The model is clearly inappropriate.
  \item The form of the curves suggests at least a quadratic in
    \code{age}.  Once you see the plot it is obvious why this should
    be so.  (They don't say what age correspond to 0 on this scale but
    my guess is about 25 years of age.)
  \item The urban versus rural differences may be additive.
  \item It appears that the \code{livch} factor could be dichotomized
    into ``0'' versus ``1 or more''.
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Preliminary model fit}
<<cm1,echo=FALSE>>=
print(cm1 <- glmer(use ~ age + I(age^2) + urban + livch + (1|district), 
                  Contraception, binomial), corr = FALSE)
@ 
\end{frame}

\begin{frame}
  \frametitle{Comments on the model fit}
  \begin{itemize}
  \item There is a highly significant quadratic term in \code{age}.
  \item The linear term in \code{age} is not significant but we retain
    it because the \code{age} scale has been centered at an arbitrary
    (and unknown) value.
  \item The \code{urban} factor is highly significant (as indicated by
    the plot).
  \item Levels of \code{livch} greater than 0 are significantly
    different from 0 but may not be different from each other.
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{Interpreting coefficient estimates}
  \begin{itemize}
  \item We are using the logit link, which is the natural link for the
    Bernoulli.
  \item For the logit link the coefficients apply to the linear
    predictor of the log-odds.
  \item The intercept is the predicted log-odds for a woman with
    centered age of 0 (we expect this means an age of 25 or so), not
    in an urban environment and with 0 live children.  As an odds
    ratio this is $e^{-1.035}=0.355$ and a probability of $0.262$.
  \item For a dichotomous factor like \code{urban} the coefficient
    $0.697$ is the increase in the log-odds for urban versus rural
    when other covariates are held fixed.
  \item This corresponds to multiplication of the odds-ratio by
    $e^{0.697}=2.008$.  The predicted probability of contraception use
    for a woman with centered age 0 in an urban environment with no
    live children is $0.416$ (it is the odds-ratio, not the
    probability, that is multiplied by $2.008$)
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Consider dichotomizing livch to 0/1+}
<<Contra2,fig=TRUE,echo=FALSE,height=7.2>>=
Contraception$ch <- factor(Contraception$livch != 0, labels = c("N","Y"))
print(xyplot(ifelse(use == "Y", 1, 0) ~ age|urban, Contraception,
             groups = ch, type = c("g", "smooth"),
             auto.key = list(space = "top", points = FALSE,
             lines = TRUE, columns = 2),
             ylab = "Proportion", xlab = "Centered age"))
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Reduced model with dichotomized livch}
<<cm2,echo=TRUE>>=
print(cm2 <- glmer(use ~ age + I(age^2) + urban + ch + (1|district),
                  Contraception, binomial), corr = FALSE)
@
\end{frame}

\begin{frame}[fragile]
  \frametitle{Comparing the model fits}
  \begin{itemize}
  \item A likelihood ratio test can be used to compare these nested models.
  \end{itemize}
<<anovac>>=  
anova(cm2, cm1)
@ 
\begin{itemize}
\item The large p-value indicates that we would not reject \code{cm2}
  in favor of \code{cm1} hence we prefer the more parsimonious \code{cm2}.
\item The plot of the scatterplot smoothers according to live children
  or none indicates that there may be a difference in the age pattern
  between these two groups.
\end{itemize}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Allowing age pattern to vary with ch}
<<cm3,echo=FALSE>>=
print(cm3 <- glmer(use ~ age*ch + I(age^2) + urban + (1|district),
              Contraception, binomial), corr = FALSE)
@
\end{frame}

\begin{frame}[fragile]
  \frametitle{Prediction intervals on the random effects}
<<ContraCat,fig=TRUE,echo=FALSE,height=7.5>>=
print(dotplot(ranef(cm3, post = TRUE),strip=FALSE)$district)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Extending the random effects}
  \begin{itemize}
  \item We may want to consider allowing a random effect for
    urban/rural by district. This is complicated by the fact the many
    districts only have rural women in the study
  \end{itemize}
<<urbanRural,echo=FALSE>>=  
cat(head(capture.output(xtabs(~urban+district, Contraception)),7),sep='\n')
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Including a random effect for urban by district}
<<cm4,echo=FALSE>>=
print(cm4 <- glmer(use ~ age*ch + I(age^2) + urban + (urban|district),
                   Contraception, binomial), corr = FALSE)
@
\end{frame}

\begin{frame}[fragile]
  \frametitle{Significance of the additional random effect}
<<anovacm4>>=  
anova(cm4,cm3)
@ 
\begin{itemize}
\item The additional random effect is highly significant in this test.
\item Most of the prediction intervals still overlap zero.
\item A scatterplot of the random effects shows several random effects
  vectors falling along a straight line.  These are the districts with
  all rural women or all urban women.
\end{itemize}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Prediction intervals for the bivariate random effects}
<<ContraCat2,fig=TRUE,echo=FALSE>>=
print(dotplot(ranef(cm4, post = TRUE))[["district"]])
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Scatter plot of the conditional modes}
<<ContraSc,fig=TRUE,echo=FALSE>>=
print(plot(ranef(cm4), type = c("g","p"), aspect = 1)[["district"]])
@ 
\end{frame}

\begin{frame}
  \frametitle{Conclusions from the example}
  \begin{itemize}
  \item Again, carefully plotting the data is enormously helpful in
    formulating the model.
  \item Observational data tend to be unbalanced and have many more
    covariates than data from a designed experiment.  Formulating a
    model is often more difficult than in a designed experiment.
  \item A generalized linear model family, typically
    \code{binomial} or \code{poisson}, is specified as the
    \code{family} argument in the call to \code{glmer}.
  \item We use likelihood-ratio tests and z-tests in the model
    building.
  \end{itemize}
\end{frame}