\documentclass{chapman}

%%% copy Sweave.sty definitions

%%% keeps `sweave' from adding `\usepackage{Sweave}': DO NOT REMOVE
%\usepackage{Sweave} 


\RequirePackage[T1]{fontenc}
\RequirePackage{graphicx,ae,fancyvrb}
\IfFileExists{upquote.sty}{\RequirePackage{upquote}}{}
\usepackage{relsize}

\DefineVerbatimEnvironment{Sinput}{Verbatim}{}
\DefineVerbatimEnvironment{Soutput}{Verbatim}{fontfamily=courier,
                                              fontshape=it,
                                              fontsize=\relsize{-1}}
\DefineVerbatimEnvironment{Scode}{Verbatim}{}
\newenvironment{Schunk}{}{}

%%% environment for raw output
\newcommand{\SchunkRaw}{\renewenvironment{Schunk}{}{}
    \DefineVerbatimEnvironment{Soutput}{Verbatim}{fontfamily=courier,
                                                  fontshape=it,
                                                  fontsize=\small}
    \rawSinput
}

%%% environment for labeled output
\newcommand{\nextcaption}{}
\newcommand{\SchunkLabel}{
  \renewenvironment{Schunk}{\begin{figure}[ht] }{\caption{\nextcaption}
  \end{figure} }
  \DefineVerbatimEnvironment{Sinput}{Verbatim}{frame = topline}
  \DefineVerbatimEnvironment{Soutput}{Verbatim}{frame = bottomline, 
                                                samepage = true,
                                                fontfamily=courier,
                                                fontshape=it,
                                                fontsize=\relsize{-1}}
}


%%% S code with line numbers
\DefineVerbatimEnvironment{Sinput}
{Verbatim}
{
%%  numbers=left
}

\newcommand{\numberSinput}{
    \DefineVerbatimEnvironment{Sinput}{Verbatim}{numbers=left}
}
\newcommand{\rawSinput}{
    \DefineVerbatimEnvironment{Sinput}{Verbatim}{}
}


%%% R / System symbols
\newcommand{\R}{\textsf{R}}
\newcommand{\rR}{{R}}
\renewcommand{\S}{\textsf{S}}
\newcommand{\SPLUS}{\textsf{S-PLUS}}
\newcommand{\rSPLUS}{{S-PLUS}}
\newcommand{\SPSS}{\textsf{SPSS}}
\newcommand{\EXCEL}{\textsf{Excel}}
\newcommand{\ACCESS}{\textsf{Access}}
\newcommand{\SQL}{\textsf{SQL}}
%%\newcommand{\Rpackage}[1]{\hbox{\rm\textit{#1}}}
%%\newcommand{\Robject}[1]{\hbox{\rm\texttt{#1}}}
%%\newcommand{\Rclass}[1]{\hbox{\rm\textit{#1}}}
%%\newcommand{\Rcmd}[1]{\hbox{\rm\texttt{#1}}}
\newcommand{\Rpackage}[1]{\index{#1 package@{\fontseries{b}\selectfont #1} package} {\fontseries{b}\selectfont #1}}
\newcommand{\rpackage}[1]{{\fontseries{b}\selectfont #1}}
\newcommand{\Robject}[1]{\texttt{#1}}
\newcommand{\Rclass}[1]{\index{#1 class@\textit{#1} class}\textit{#1}}
\newcommand{\Rcmd}[1]{\index{#1 function@\texttt{#1} function}\texttt{#1}}
\newcommand{\Roperator}[1]{\texttt{#1}}
\newcommand{\Rarg}[1]{\texttt{#1}}
\newcommand{\Rlevel}[1]{\texttt{#1}}


%%% other symbols
\newcommand{\file}[1]{\hbox{\rm\texttt{#1}}}
%%\newcommand{\stress}[1]{\index{#1}\textit{#1}} 
\newcommand{\stress}[1]{\textit{#1}} 
\newcommand{\booktitle}[1]{\textit{#1}} %%'

%%% Math symbols
\usepackage{amstext}
\usepackage{amsmath}
\newcommand{\E}{\mathsf{E}}   
\newcommand{\Var}{\mathsf{Var}}   
\newcommand{\Cov}{\mathsf{Cov}}   
\newcommand{\Cor}{\mathsf{Cor}}   
\newcommand{\x}{\mathbf{x}}   
\newcommand{\y}{\mathbf{y}}   
\renewcommand{\a}{\mathbf{a}}
\newcommand{\W}{\mathbf{W}}   
\newcommand{\C}{\mathbf{C}}   
\renewcommand{\H}{\mathbf{H}}   
\newcommand{\X}{\mathbf{X}}   
\newcommand{\B}{\mathbf{B}}   
\newcommand{\V}{\mathbf{V}}   
\newcommand{\I}{\mathbf{I}}   
\newcommand{\D}{\mathbf{D}}   
\newcommand{\bS}{\mathbf{S}}   
\newcommand{\N}{\mathcal{N}}   
\renewcommand{\L}{L}   
\renewcommand{\P}{\mathsf{P}}   
\newcommand{\K}{\mathbf{K}}
\newcommand{\m}{\mathbf{m}}
\newcommand{\argmin}{\operatorname{argmin}\displaylimits}
\newcommand{\argmax}{\operatorname{argmax}\displaylimits}

\newcommand{\bx}{\mathbf{x}}
\newcommand{\bbeta}{\mathbf{\beta}}


%%% links
\usepackage{hyperref}

\hypersetup{%
  pdftitle = {A Handbook of Statistical Analyses Using R (3rd Edition)},
  pdfsubject = {Book},
  pdfauthor = {Torsten Hothorn and Brian S. Everitt},
  colorlinks = {black},
  linkcolor = {black},
  citecolor = {black},
  urlcolor = {black},
  hyperindex = {true},
  linktocpage = {true},
}


%%% captions & tables
%% <FIXME>: conflics with figure definition in chapman.cls
%%\usepackage[format=hang,margin=10pt,labelfont=bf]{caption}
%% </FIMXE>
\usepackage{longtable}
\usepackage[figuresright]{rotating}

%%% R symbol in chapter 1
\usepackage{wrapfig}

%%% Bibliography
\usepackage[round,comma]{natbib}
\renewcommand{\refname}{References \addcontentsline{toc}{chapter}{References}}
\citeindexfalse

%%% texi2dvi complains that \newblock is undefined, hm...
\def\newblock{\hskip .11em plus .33em minus .07em}

%%% Example sections
\newcounter{exercise}[chapter]
\setcounter{exercise}{0}
\newcommand{\exercise}{\stepcounter{exercise} \item{Ex.~\arabic{chapter}.\arabic{exercise} }}


%% URLs
\newcommand{\curl}[1]{\begin{center} \url{#1} \end{center}}

%%% for manual corrections
%\renewcommand{\baselinestretch}{2}

%%% plot sizes
\setkeys{Gin}{width=0.95\textwidth}

%%% color
\usepackage{color}

%%% hyphenations
\hyphenation{drop-out}
\hyphenation{mar-gi-nal}

%%% new bidirectional quotes need 
\usepackage[utf8]{inputenc}

%\usepackage{setspace}
\definecolor{sidebox_todo}{rgb}{1,1,0.2}
\newcommand{\todo}[1]{
        \hspace{0pt}%
        \marginpar{% 
                        \fcolorbox{black}{sidebox_todo}{%
                                \parbox{\marginparwidth} {

\raggedright\sffamily\footnotesize{TODO: #1}%
                                }
                        }%
        }
}
\begin{document}

%% Title page

\title{A Handbook of Statistical Analyses Using \R{} --- 3rd Edition}

\author{Torsten Hothorn and Brian S. Everitt}

\maketitle
%%\VignetteIndexEntry{Chapter Analyzing Longitudinal Data I}
%%\VignetteDepends{lme4,multcomp}
\setcounter{chapter}{12}


\SweaveOpts{prefix.string=figures/HSAUR,eps=FALSE,keep.source=TRUE} 

<<setup, echo = FALSE, results = hide>>=
rm(list = ls())
s <- search()[-1]
s <- s[-match(c("package:base", "package:stats", "package:graphics", "package:grDevices",
                "package:utils", "package:datasets", "package:methods", "Autoloads"), s)]
if (length(s) > 0) sapply(s, detach, character.only = TRUE)
if (!file.exists("tables")) dir.create("tables")
if (!file.exists("figures")) dir.create("figures")
set.seed(290875)
options(prompt = "R> ", continue = "+  ",
    width = 63, # digits = 4, 
    show.signif.stars = FALSE,
    SweaveHooks = list(leftpar = function() 
        par(mai = par("mai") * c(1, 1.05, 1, 1)),
        bigleftpar = function()
        par(mai = par("mai") * c(1, 1.7, 1, 1))))
HSAURpkg <- require("HSAUR3")
if (!HSAURpkg) stop("cannot load package ", sQuote("HSAUR3"))
rm(HSAURpkg)
 ### </FIXME> hm, R-2.4.0 --vanilla seems to need this
a <- Sys.setlocale("LC_ALL", "C")
 ### </FIXME>
book <- TRUE
refs <- cbind(c("AItR", "DAGD", "SI", "CI", "ANOVA", "MLR", "GLM", 
                "DE", "RP", "GAM", "SA", "ALDI", "ALDII", "SIMC", "MA", "PCA", 
                "MDS", "CA"), 1:18)
ch <- function(x) {
    ch <- refs[which(refs[,1] == x),]
    if (book) {
        return(paste("Chapter~\\\\ref{", ch[1], "}", sep = ""))
    } else {
        return(paste("Chapter~", ch[2], sep = ""))
    }
}
if (file.exists("deparse.R"))
    source("deparse.R")

setHook(packageEvent("lattice", "attach"), function(...) {
    lattice.options(default.theme = 
        function()
            standard.theme("pdf", color = FALSE))
    })
@

\pagestyle{headings}
<<singlebook, echo = FALSE>>=
book <- FALSE
@

<<ALDI-setup, echo = FALSE, results = hide>>=
library("lme4")
library("multcomp")
residuals <- function(object) {
    y <- getME(object, 'y')
    y - fitted(object)
}
@

\chapter[Analyzing Longitudinal Data I]{Analyzing Longitudinal Data I: 
         Computerized Delivery of Cognitive
         Behavioral Therapy -- Beat the Blues \label{ALDI}}

\section{Introduction}



\section{Analyzing Longitudinal Data}


\section{Analysis Using \R{}}


\begin{figure}
\begin{center}
<<ALDI-plot-BtheB, echo = TRUE, fig = TRUE, height = 4>>=
data("BtheB", package = "HSAUR3")
layout(matrix(1:2, nrow = 1))
ylim <- range(BtheB[,grep("bdi", names(BtheB))], 
              na.rm = TRUE)
tau <- subset(BtheB, treatment == "TAU")[,
    grep("bdi", names(BtheB))]
boxplot(tau, main = "Treated as Usual", ylab = "BDI", 
        xlab = "Time (in months)", names = c(0, 2, 3, 5, 8), 
        ylim = ylim)
btheb <- subset(BtheB, treatment == "BtheB")[,
    grep("bdi", names(BtheB))]
boxplot(btheb, main = "Beat the Blues", ylab = "BDI", 
        xlab = "Time (in months)", names = c(0, 2, 3, 5, 8), 
        ylim = ylim)
@
\caption{Boxplots for the repeated measures by treatment group for the 
         \Robject{BtheB} data. \label{ALDI:boxplots}}
\end{center}
\end{figure}


We shall fit both random intercept and random intercept and 
slope models to the data including the baseline BDI values
(\Robject{pre.bdi}), \Robject{treatment} group, 
\Robject{drug}, and \Robject{length} as fixed effect covariates. Linear mixed 
effects models are fitted in \R{} by using the \Rcmd{lmer} function 
contained in the \Rpackage{lme4} package
\citep{PKG:lme4,HSAUR:PinheiroBates2000,HSAUR:Bates2005},
but an essential first step is to rearrange the data from the `wide form' in which they appear in the \Robject{BtheB} data frame %%'
into the `long form' in which each separate repeated measurement %%'
and associated covariate values appear as a separate row in a
\Rclass{data.frame}. 
This rearrangement can be made using the following code:
<<ALDI-long-BtheB, echo = TRUE>>=
data("BtheB", package = "HSAUR3")
BtheB$subject <- factor(rownames(BtheB))
nobs <- nrow(BtheB)
BtheB_long <- reshape(BtheB, idvar = "subject", 
    varying = c("bdi.2m", "bdi.3m", "bdi.5m", "bdi.8m"),
    direction = "long") 
BtheB_long$time <- rep(c(2, 3, 5, 8), rep(nobs, 4)) 
@
such that the data are now in the form (here shown for the first three
subjects)
<<ALDI-showlong-BtheB, echo = TRUE>>=
subset(BtheB_long, subject %in% c("1", "2", "3"))
@
The resulting \Rclass{data.frame} \Robject{BtheB\_long} contains a number 
of missing values 
\index{Missing values}
and in applying the \Rcmd{lmer} function these 
will be dropped. But notice it is only the missing values 
that are removed, \stress{not} participants that have at least 
one missing value. All the available data is used in the model 
fitting process. The \Rcmd{lmer} function is used in a similar way to 
the \Rcmd{lm} function met in \Sexpr{ch("MLR")}
with the addition of a random term to identify the source of the 
repeated measurements, here \Robject{subject}. We can fit the two 
models (\ref{ALDI:ModelA}) and (\ref{ALDI:ModelB}) and test which is most 
appropriate using
<<ALDI-fit-BtheB, echo = TRUE>>=
library("lme4")
BtheB_lmer1 <- lmer(bdi ~ bdi.pre + time + treatment + drug + 
    length + (1 | subject), data = BtheB_long, 
    REML = FALSE, na.action = na.omit)
BtheB_lmer2 <- lmer(bdi ~ bdi.pre + time + treatment + drug + 
    length + (time | subject), data = BtheB_long, 
    REML = FALSE, na.action = na.omit)   
anova(BtheB_lmer1, BtheB_lmer2)
@


\renewcommand{\nextcaption}{\R{} output of the linear mixed-effects model fit
                            for the \Robject{BtheB} data.
                            \label{ALDI-BtheB-summary}}
\SchunkLabel
<<ALDI-summary-BtheB, echo = TRUE>>=
summary(BtheB_lmer1)
@
\SchunkRaw

The \Rcmd{summary} method for \Rclass{lmer} objects doesn't print $p$-values
for Gaussian mixed models because the degrees of freedom of the $t$ 
reference distribution are not obvious. However, one can rely on the 
asymptotic normal distribution for computing univariate $p$-values for the 
fixed effects using the \Rcmd{cftest} function from package \Rpackage{multcomp}. 
The asymptotic $p$-values are given in Figure~\ref{ALDI-BtheB-summary-p}.

\renewcommand{\nextcaption}{\R{} output of the asymptotic $p$-values for linear mixed-effects model fit
                            for the \Robject{BtheB} data.
                            \label{ALDI-BtheB-summary-p}}
\SchunkLabel
<<ALDI-summary-BtheB-p, echo = TRUE>>=
cftest(BtheB_lmer1)
@
\SchunkRaw

We can check the assumptions of the final model fitted to
the \Robject{BtheB} data, i.e., the normality of the random effect terms 
and the residuals, by first using the \Rcmd{ranef} method 
to \stress{predict} the former and the \Rcmd{residuals} method  
to calculate the differences between the observed data values 
and the fitted values, and then using normal probability plots 
on each. How the random effects are predicted is explained briefly in 
Section~\ref{ALDI:predictrandom}. 
The necessary \R{} code to obtain the effects, residuals, 
and plots is shown with Figure~\ref{ALDI:qqplots}.
There appear to be no large departures from linearity in either plot.
\begin{figure}
\begin{center}
<<ALDI-qqnorm-BtheB, echo = TRUE, fig = TRUE, height = 4, width = 7>>=
layout(matrix(1:2, ncol = 2))
qint <- ranef(BtheB_lmer1)$subject[["(Intercept)"]]
qres <- residuals(BtheB_lmer1)
qqnorm(qint, ylab = "Estimated random intercepts",
       xlim = c(-3, 3), ylim = c(-20, 20),
       main = "Random intercepts")
qqline(qint)
qqnorm(qres, xlim = c(-3, 3), ylim = c(-20, 20),
       ylab = "Estimated residuals",
       main = "Residuals")
qqline(qres)
@
\caption{Quantile-quantile plots of predicted random intercepts and
residuals for the random intercept model \Robject{BtheB\_lmer1} fitted to the
\Robject{BtheB} data. \label{ALDI:qqplots}}
\end{center}
\end{figure}


\begin{figure}
\begin{center}
<<ALDI-dropout, echo = TRUE, fig = TRUE>>=
bdi <- BtheB[, grep("bdi", names(BtheB))]
plot(1:4, rep(-0.5, 4), type = "n", axes = FALSE, 
     ylim = c(0, 50), xlab = "Months", ylab = "BDI")
axis(1, at = 1:4, labels = c(0, 2, 3, 5))
axis(2)
for (i in 1:4) {
    dropout <- is.na(bdi[,i + 1])
    points(rep(i, nrow(bdi)) + ifelse(dropout, 0.05, -0.05), 
           jitter(bdi[,i]), pch = ifelse(dropout, 20, 1))
}
@
\caption{Distribution of BDI values for patients
         that do (circles) and do not (bullets) attend the next scheduled visit. 
         \label{ALDI-dropout}}
\end{center}
\end{figure}



\bibliographystyle{LaTeXBibTeX/refstyle}
\bibliography{LaTeXBibTeX/HSAUR}
\end{document}