\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}