% Options for packages loaded elsewhere
\PassOptionsToPackage{unicode}{hyperref}
\PassOptionsToPackage{hyphens}{url}
%
\documentclass[
]{article}
\usepackage{amsmath,amssymb}
\usepackage{lmodern}
\usepackage{iftex}
\ifPDFTeX
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
  \usepackage{textcomp} % provide euro and other symbols
\else % if luatex or xetex
  \usepackage{unicode-math}
  \defaultfontfeatures{Scale=MatchLowercase}
  \defaultfontfeatures[\rmfamily]{Ligatures=TeX,Scale=1}
\fi
% Use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\IfFileExists{microtype.sty}{% use microtype if available
  \usepackage[]{microtype}
  \UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\makeatletter
\@ifundefined{KOMAClassName}{% if non-KOMA class
  \IfFileExists{parskip.sty}{%
    \usepackage{parskip}
  }{% else
    \setlength{\parindent}{0pt}
    \setlength{\parskip}{6pt plus 2pt minus 1pt}}
}{% if KOMA class
  \KOMAoptions{parskip=half}}
\makeatother
\usepackage{xcolor}
\usepackage[margin=1in]{geometry}
\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\usepackage{framed}
\definecolor{shadecolor}{RGB}{248,248,248}
\newenvironment{Shaded}{\begin{snugshade}}{\end{snugshade}}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{0.94,0.16,0.16}{#1}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.77,0.63,0.00}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\BuiltInTok}[1]{#1}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{0.64,0.00,0.00}{\textbf{#1}}}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ImportTok}[1]{#1}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.81,0.36,0.00}{\textbf{#1}}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{#1}}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\usepackage{graphicx}
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth\else\Gin@nat@width\fi}
\def\maxheight{\ifdim\Gin@nat@height>\textheight\textheight\else\Gin@nat@height\fi}
\makeatother
% Scale images if necessary, so that they will not overflow the page
% margins by default, and it is still possible to overwrite the defaults
% using explicit options in \includegraphics[width, height, ...]{}
\setkeys{Gin}{width=\maxwidth,height=\maxheight,keepaspectratio}
% Set default figure placement to htbp
\makeatletter
\def\fps@figure{htbp}
\makeatother
\setlength{\emergencystretch}{3em} % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{-\maxdimen} % remove section numbering
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{hyperref}
\usepackage{xcolor}
\DeclareMathOperator{\Cov}{Cov}
\DeclareMathOperator{\Corr}{Corr}
\DeclareMathOperator{\Var}{Var}
\DeclareMathOperator{\upmodels}{\perp \!\!\! \perp}
\ifLuaTeX
  \usepackage{selnolig}  % disable illegal ligatures
\fi
\IfFileExists{bookmark.sty}{\usepackage{bookmark}}{\usepackage{hyperref}}
\IfFileExists{xurl.sty}{\usepackage{xurl}}{} % add URL line breaks if available
\urlstyle{same} % disable monospaced font for URLs
\hypersetup{
  pdftitle={Homework 08 Solutions},
  hidelinks,
  pdfcreator={LaTeX via pandoc}}

\title{Homework 08 Solutions}
\author{}
\date{\vspace{-2.5em}2022-11-14}

\begin{document}
\maketitle

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(arm)}
\FunctionTok{library}\NormalTok{(ggplot2)}
\FunctionTok{library}\NormalTok{(HLMdiag)    }\DocumentationTok{\#\# easy extraction of MLM residuals}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning: package 'HLMdiag' was built under R version 4.1.2
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(gridExtra)  }\DocumentationTok{\#\# For the grid.arrange() function}

\FunctionTok{library}\NormalTok{(cAIC4)      }\DocumentationTok{\#\# I\textquotesingle{}ve asked you not to bother with cAIC for this assignment, but since}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning: package 'cAIC4' was built under R version 4.1.2
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
                    \DocumentationTok{\#\# it is on the assignment I thought it would be useful for you to see it.}
                    \DocumentationTok{\#\# This library provides a function cAIC() for calculating it.}
\end{Highlighting}
\end{Shaded}

\def\red#1{\textcolor{red}{#1}}
\def\blue#1{\textcolor{blue}{#1}}
\let\bar\overline
\let\ul\underline

\hypertarget{problem-1.}{%
\section{Problem 1.}\label{problem-1.}}

The file allvars.csv contains CD4 percentages (CD4PCT) for a set of 254
young children with HIV who were measured several times over a period of
two years (A CD4 count measures the number of CD4 cells in your blood.
It's used to check the immune system function in people with HIV
\url{https://medlineplus.gov/lab-tests/cd4-lymphocyte-count/}{]}). The
dataset also includes the ages of the children at each measurement. This
is an example of \textit{growth curve data} in which the Level 1
observations are CD4PCT's at each visit, and the Level 2 ``group'\,' is
a child. How the CD4PCT's change over time for each child is the primary
question of interest. Because of skewing, you should replace CD4PCT with
sqrt(CD4PCT) for all of the following questions.

\hypertarget{problem-1a.}{%
\subsection{Problem 1(a).}\label{problem-1a.}}

For the first 12 children only (\texttt{newpid $\leq$ 12}), make a facet
plot of CD4PCT vs VISIT number. Use different plotting point colors for
children who have treatmnt=1 and treatmnt=0 (Presumably, ``treatmnt'\,'
stands for some sort of treatment to improve childrens' HIV status.). Do
the same for CD4PCT vs.~visage (the age at the child on each visit).
(\em Note: If you are using \texttt{ggplot()}, you may wish to set \texttt{scales="free\_x"} in the \texttt{facet\_wrap()} function.\rm)  Which is a better measure of time?

\it
First we read in the data and do a little exploring and munging...
\rm

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{cd4 }\OtherTok{\textless{}{-}} \FunctionTok{read.csv}\NormalTok{(}\StringTok{"allvar.csv"}\NormalTok{,}\AttributeTok{header=}\NormalTok{T)}

\DocumentationTok{\#\# It\textquotesingle{}s worth doing a little exploration of the data here,}
\DocumentationTok{\#\# especially to look at missing values.}

\FunctionTok{apply}\NormalTok{(cd4,}\DecValTok{2}\NormalTok{,}\ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{mean}\NormalTok{(}\FunctionTok{is.na}\NormalTok{(x)))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##       VISIT      newpid       VDATE      CD4PCT         arv      visage 
## 0.000000000 0.000000000 0.000000000 0.142743222 0.103668262 0.109250399 
##    treatmnt      CD4CNT     baseage 
## 0.000000000 0.146730463 0.007177033
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# we only need the variables newpid, VISIT, visage, CD4PCT and treatmnt}
\DocumentationTok{\#\# so let\textquotesingle{}s get rid of some (but not all) missing data problems}
\DocumentationTok{\#\# by deleting the other variables}

\NormalTok{cd4 }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cd4,}\FunctionTok{data.frame}\NormalTok{(}\AttributeTok{newpid=}\NormalTok{newpid,}
                           \AttributeTok{VISIT=}\NormalTok{VISIT, }
                           \AttributeTok{visage=}\NormalTok{visage, }
                           \AttributeTok{CD4PCT=}\NormalTok{CD4PCT, }
                           \AttributeTok{treatmnt=}\NormalTok{treatmnt))}

\DocumentationTok{\#\# and since we still have some missing data, we will just }
\DocumentationTok{\#\# delete the rows that continue to have NA\textquotesingle{}s (not a great}
\DocumentationTok{\#\# practice in general, but good enough for this exercise...)}
\NormalTok{cd4 }\OtherTok{\textless{}{-}}\NormalTok{ cd4[}\SpecialCharTok{!}\FunctionTok{apply}\NormalTok{(cd4,}\DecValTok{1}\NormalTok{,}\ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{any}\NormalTok{(}\FunctionTok{is.na}\NormalTok{(x))),]}

\DocumentationTok{\#\# set up the sqrt of the response variable...}
\NormalTok{cd4}\SpecialCharTok{$}\NormalTok{sqrt.CD4PCT }\OtherTok{\textless{}{-}} \FunctionTok{sqrt}\NormalTok{(cd4}\SpecialCharTok{$}\NormalTok{CD4PCT)}

\DocumentationTok{\#\# and for coloring ggplot elements it will better if treatmnt is a factor...}
\NormalTok{cd4}\SpecialCharTok{$}\NormalTok{treatmnt }\OtherTok{\textless{}{-}} \FunctionTok{as.factor}\NormalTok{(cd4}\SpecialCharTok{$}\NormalTok{treatmnt)}

\DocumentationTok{\#\# select the first 12 kids...}
\NormalTok{first}\FloatTok{.12} \OtherTok{\textless{}{-}}\NormalTok{ (cd4}\SpecialCharTok{$}\NormalTok{newpid }\SpecialCharTok{\textless{}=} \DecValTok{12}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\it
I will show the faceted graph for both \verb|VISIT| and \verb|visage| as the time variable\ldots{}

Here it is using VISIT\ldots \rm

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{g }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(cd4[first}\FloatTok{.12}\NormalTok{,],}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{VISIT,}\AttributeTok{y=}\NormalTok{sqrt.CD4PCT)) }\SpecialCharTok{+}
  \FunctionTok{facet\_wrap}\NormalTok{( }\SpecialCharTok{\textasciitilde{}}\NormalTok{ newpid) }\SpecialCharTok{+} 
  \FunctionTok{geom\_point}\NormalTok{(}\FunctionTok{aes}\NormalTok{(}\AttributeTok{color=}\NormalTok{treatmnt))}

\NormalTok{g}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-3-1.pdf}
\clearpage

\it
\ldots{}and here it is using \verb|visage| as the time variable\ldots

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{h }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(cd4[first}\FloatTok{.12}\NormalTok{,],}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{visage,}\AttributeTok{y=}\NormalTok{sqrt.CD4PCT)) }\SpecialCharTok{+}
  \FunctionTok{facet\_wrap}\NormalTok{( }\SpecialCharTok{\textasciitilde{}}\NormalTok{ newpid, }\AttributeTok{scales=}\StringTok{"free\_x"}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{geom\_point}\NormalTok{(}\FunctionTok{aes}\NormalTok{(}\AttributeTok{color=}\NormalTok{treatmnt))}

\NormalTok{h}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-4-1.pdf}

\it
The graphs look very similar, but \verb|visage|, which is the actual age of the child, would seem to have more physical/medical meaning than \verb|VISIT|, which is just the sequential visit number of each visit.  If you have a reason for preferring \verb|VISIT|, go for it.

I will show the solutions below both ways (it will turn out not to make
much difference!). \rm

\hypertarget{problem-1b.}{%
\subsection{Problem 1(b).}\label{problem-1b.}}

Build a multilevel model with random slopes and intercepts for all the
children, using a measure of time (VISIT or visage, whichever you
answered in part (a)) as a Level 1 predictor, and treatmnt as a Level 2
predictor. Report the estimated fixed effects and variances from the
model, and add the fitted regression lines to your plot from part (a).

\it 
Here's the model, using \verb|VISIT| as the measure of time:
\rm

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{display}\NormalTok{(lmer}\FloatTok{.1} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{+}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid), }\AttributeTok{data=}\NormalTok{cd4))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## lmer(formula = sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1 + VISIT | 
##     newpid), data = cd4)
##             coef.est coef.se
## (Intercept)  4.69     0.13  
## VISIT       -0.03     0.01  
## treatmnt2    0.18     0.18  
## 
## Error terms:
##  Groups   Name        Std.Dev. Corr  
##  newpid   (Intercept) 1.40           
##           VISIT       0.05     -0.10 
##  Residual             0.72           
## ---
## number of obs: 1075, groups: newpid, 251
## AIC = 3140.7, DIC = 3100.6
## deviance = 3113.6
\end{verbatim}

\it
The overall intercept is $\hat\beta_0=4.69$, the overall slope on \verb|VISIT| is $\beta_1=-0.03$ and the overall treatment effect is $\beta_2=0.18$.  The residual variances are $\hat\sigma^2=(0.51)^2=0.26$, $\hat\tau_0^2=(1.40)^2=1.96$ and $\hat\tau_1^2=(0.05)^2=0.0025$, and $\hat\rho$ is a very modest $-0.10$. The only things that seem surprising are

\begin{itemize}
\item $\hat\tau_0 = 1.40$ seems large, considering that the most \verb|sqrt.CD4PCT| values are between zero and 5;
\item The treatmt effect, $\hat\beta_2=0.18$, is not significantly different from zero, even though the scatter plots of the first 12 children would suggest otherwise!
\end{itemize}

Next we overlay regression lines on the scatterplots from part (a). The
regression lines track the data well (and even have different
intercepts, even though the the \verb|treatmnt| effect is nonsignificant
in the model summary! {[}what part of the model could be helping here,
do you suppose??{]})

(It was a little trickier than I expected, to overlay the regression
lines and have them colored differntly for each of the two
\verb|treatmnt| levels, so have a look at the code below also.)

Here are the overlaid regression lines, using \verb|VISIT| as the time
variable: \rm

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# fixef(lmer.1)  \#\# gives the betas for intercept and slopes on VISIT \& treatmnt}

\DocumentationTok{\#\# ranef(lmer.1)  \#\# gives the etas for intercepts and slope on VISIT}

\NormalTok{coef}\FloatTok{.1} \OtherTok{\textless{}{-}} \FunctionTok{coef}\NormalTok{(lmer}\FloatTok{.1}\NormalTok{)}\SpecialCharTok{$}\NormalTok{newpid[,}\DecValTok{1}\SpecialCharTok{:}\DecValTok{3}\NormalTok{] }\DocumentationTok{\#\# coef.1[,1] = beta0 + eta0  (random intcpt)}
                                    \DocumentationTok{\#\# coef.1[,2] = beta1 + eta1  (random slope)}
                                    \DocumentationTok{\#\# coef.1[,3] = beta2         (coef on treatmnt)}

\NormalTok{tx }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cd4,}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{split}\NormalTok{(treatmnt,newpid), }\ControlFlowTok{function}\NormalTok{(x) x[}\DecValTok{1}\NormalTok{]))  }\DocumentationTok{\#\# one per child}

\NormalTok{params }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{newpid=}\FunctionTok{sort}\NormalTok{(}\FunctionTok{unique}\NormalTok{(cd4}\SpecialCharTok{$}\NormalTok{newpid)),}
                     \AttributeTok{alpha0=}\NormalTok{coef}\FloatTok{.1}\NormalTok{[,}\DecValTok{1}\NormalTok{],}
                     \AttributeTok{alpha1=}\NormalTok{coef}\FloatTok{.1}\NormalTok{[,}\DecValTok{2}\NormalTok{],}
                     \AttributeTok{beta2=}\NormalTok{coef}\FloatTok{.1}\NormalTok{[,}\DecValTok{3}\NormalTok{],}
\NormalTok{                     tx }
\NormalTok{                     )}

\NormalTok{g }\SpecialCharTok{+} \FunctionTok{geom\_abline}\NormalTok{(}\AttributeTok{data=}\NormalTok{params[params}\SpecialCharTok{$}\NormalTok{newpid}\SpecialCharTok{\textless{}=}\DecValTok{12}\NormalTok{,],}\FunctionTok{aes}\NormalTok{(}\AttributeTok{intercept=}\NormalTok{alpha0 }\SpecialCharTok{+} \FunctionTok{ifelse}\NormalTok{(tx}\SpecialCharTok{==}\DecValTok{2}\NormalTok{,beta2,}\DecValTok{0}\NormalTok{),}
                                                    \AttributeTok{slope=}\NormalTok{alpha1,}
                                                    \AttributeTok{color=}\NormalTok{tx))}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-6-1.pdf}

\it
\mbox{}\\
Here's the model summary using \verb|visage| as the measure of time:
\mbox{}\\
\rm

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{display}\NormalTok{(lmer}\FloatTok{.2} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ visage }\SpecialCharTok{+}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{visage}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## lmer(formula = sqrt.CD4PCT ~ 1 + visage + treatmnt + (1 + visage | 
##     newpid), data = cd4)
##             coef.est coef.se
## (Intercept)  5.34     0.18  
## visage      -0.23     0.04  
## treatmnt2    0.23     0.18  
## 
## Error terms:
##  Groups   Name        Std.Dev. Corr  
##  newpid   (Intercept) 1.45           
##           visage      0.29     -0.53 
##  Residual             0.75           
## ---
## number of obs: 1075, groups: newpid, 251
## AIC = 3152, DIC = 3119.5
## deviance = 3128.8
\end{verbatim}

\it
The story is much the same using \verb|visage| as a measure time, as using \verb|VISIT| as a measure of time.  Again $\hat\tau_0$ is rather large, and the estimated \verb|treatmnt| effect $\hat\beta_2$ is nonsignificant.

Here's the facets plot for the first 12 children, using \verb|visage|
for time. Again the regression lines are tracking differences in the
children rather well, even though the \verb|treatmnt| effect is
nonsignificant. \rm

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{coef}\FloatTok{.2} \OtherTok{\textless{}{-}} \FunctionTok{coef}\NormalTok{(lmer}\FloatTok{.2}\NormalTok{)}\SpecialCharTok{$}\NormalTok{newpid[,}\DecValTok{1}\SpecialCharTok{:}\DecValTok{3}\NormalTok{]  }\DocumentationTok{\#\# coef.2[,1] = beta0 + eta0  (random intcpt)}
                                     \DocumentationTok{\#\# coef.2[,2] = beta1 + eta1  (random slope)}
                                     \DocumentationTok{\#\# coef.2[,3] = beta2         (coef on treatmnt)}

\NormalTok{tx }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cd4,}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{split}\NormalTok{(treatmnt,newpid), }\ControlFlowTok{function}\NormalTok{(x) x[}\DecValTok{1}\NormalTok{]))  }\DocumentationTok{\#\# one per child}

\NormalTok{params }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{newpid=}\FunctionTok{sort}\NormalTok{(}\FunctionTok{unique}\NormalTok{(cd4}\SpecialCharTok{$}\NormalTok{newpid)),}
                     \AttributeTok{alpha0=}\NormalTok{coef}\FloatTok{.2}\NormalTok{[,}\DecValTok{1}\NormalTok{],}
                     \AttributeTok{alpha1=}\NormalTok{coef}\FloatTok{.2}\NormalTok{[,}\DecValTok{2}\NormalTok{],}
                     \AttributeTok{beta2=}\NormalTok{coef}\FloatTok{.2}\NormalTok{[,}\DecValTok{3}\NormalTok{],}
\NormalTok{                     tx }
\NormalTok{                     )}

\NormalTok{h }\SpecialCharTok{+} \FunctionTok{geom\_abline}\NormalTok{(}\AttributeTok{data=}\NormalTok{params[params}\SpecialCharTok{$}\NormalTok{newpid}\SpecialCharTok{\textless{}=}\DecValTok{12}\NormalTok{,],}\FunctionTok{aes}\NormalTok{(}\AttributeTok{intercept=}\NormalTok{alpha0 }\SpecialCharTok{+} \FunctionTok{ifelse}\NormalTok{(tx}\SpecialCharTok{==}\DecValTok{2}\NormalTok{,beta2,}\DecValTok{0}\NormalTok{),}
                                                    \AttributeTok{slope=}\NormalTok{alpha1,}
                                                    \AttributeTok{color=}\NormalTok{tx))}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-8-1.pdf}

\hypertarget{problem-1c.}{%
\subsection{Problem 1(c).}\label{problem-1c.}}

Try to expand the model by adding an interaction between treatmnt and
your measure of time (Beacuse treatmnt is a level-2 variable and the
time measure is a level-1 variable, this is somethings called a
\textit{cross-level interaction}.). Comment on changes in the estimated
fixed effects and variances for this model, vs.~the model in part (b).
Add the fitted regression lines from this model to your previous facets
plot. Does treatmnt account for an interesting amount of the variation
in slopes (or intercepts) among children in the study?

\it
OK, here we go with the crosslevel interaction model using \verb|VISIT| as the time variable:
\rm

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{display}\NormalTok{(lmer}\FloatTok{.3} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{*}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid), }\AttributeTok{data=}\NormalTok{cd4))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## lmer(formula = sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (1 + VISIT | 
##     newpid), data = cd4)
##                 coef.est coef.se
## (Intercept)      4.71     0.13  
## VISIT           -0.03     0.01  
## treatmnt2        0.14     0.19  
## VISIT:treatmnt2  0.01     0.01  
## 
## Error terms:
##  Groups   Name        Std.Dev. Corr  
##  newpid   (Intercept) 1.40           
##           VISIT       0.05     -0.10 
##  Residual             0.72           
## ---
## number of obs: 1075, groups: newpid, 251
## AIC = 3149.4, DIC = 3093
## deviance = 3113.2
\end{verbatim}

\it
Same story as before: neither the coefficient on \verb|treatmnt| nor the coefficient on the interaction is significantly different from zero, but $\hat\tau_0$ is somewhat large for the scale of the data.

Here are the regression lines plotted over the scatter plots using
\verb|VISIT| for time; again you can see that they are tracking
differences between children that we might have thought would be
detected by the treatment effect. \rm

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{coef}\FloatTok{.3} \OtherTok{\textless{}{-}} \FunctionTok{coef}\NormalTok{(lmer}\FloatTok{.3}\NormalTok{)}\SpecialCharTok{$}\NormalTok{newpid[,}\DecValTok{1}\SpecialCharTok{:}\DecValTok{4}\NormalTok{]  }\DocumentationTok{\#\# coef.3[,1] = beta0 + eta0  (random intcpt)}
                                     \DocumentationTok{\#\# coef.3[,2] = beta1 + eta1  (random slope on time)}
                                     \DocumentationTok{\#\# coef.3[,3] = beta2         (coef on treatmnt)}
                                     \DocumentationTok{\#\# coef.3[,4] = beta3         (coef on interaction)}

\NormalTok{tx }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cd4,}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{split}\NormalTok{(treatmnt,newpid), }\ControlFlowTok{function}\NormalTok{(x) x[}\DecValTok{1}\NormalTok{]))  }\DocumentationTok{\#\# one per child}

\NormalTok{params }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{newpid=}\FunctionTok{sort}\NormalTok{(}\FunctionTok{unique}\NormalTok{(cd4}\SpecialCharTok{$}\NormalTok{newpid)),}
                     \AttributeTok{alpha0=}\NormalTok{coef}\FloatTok{.3}\NormalTok{[,}\DecValTok{1}\NormalTok{],}
                     \AttributeTok{alpha1=}\NormalTok{coef}\FloatTok{.3}\NormalTok{[,}\DecValTok{2}\NormalTok{],}
                     \AttributeTok{beta2=}\NormalTok{coef}\FloatTok{.3}\NormalTok{[,}\DecValTok{3}\NormalTok{],}
                     \AttributeTok{beta3=}\NormalTok{coef}\FloatTok{.3}\NormalTok{[,}\DecValTok{4}\NormalTok{],}
\NormalTok{                     tx }
\NormalTok{                     )}

\NormalTok{g }\SpecialCharTok{+} \FunctionTok{geom\_abline}\NormalTok{(}\AttributeTok{data=}\NormalTok{params[params}\SpecialCharTok{$}\NormalTok{newpid}\SpecialCharTok{\textless{}=}\DecValTok{12}\NormalTok{,],}\FunctionTok{aes}\NormalTok{(}\AttributeTok{intercept=}\NormalTok{alpha0 }\SpecialCharTok{+} \FunctionTok{ifelse}\NormalTok{(tx}\SpecialCharTok{==}\DecValTok{2}\NormalTok{,beta2,}\DecValTok{0}\NormalTok{),}
                                                    \AttributeTok{slope=}\NormalTok{alpha1 }\SpecialCharTok{+} \FunctionTok{ifelse}\NormalTok{(tx}\SpecialCharTok{==}\DecValTok{2}\NormalTok{,beta3,}\DecValTok{0}\NormalTok{),}
                                                    \AttributeTok{color=}\NormalTok{tx))}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-10-1.pdf}

\it
\mbox{}\\
Here's the cross level interaction model, using \verb|visage| as time
\mbox{}\\
\rm

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{display}\NormalTok{(lmer}\FloatTok{.4} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ visage }\SpecialCharTok{*}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{visage}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## lmer(formula = sqrt.CD4PCT ~ 1 + visage * treatmnt + (1 + visage | 
##     newpid), data = cd4)
##                  coef.est coef.se
## (Intercept)       5.34     0.21  
## visage           -0.23     0.05  
## treatmnt2         0.24     0.32  
## visage:treatmnt2  0.00     0.08  
## 
## Error terms:
##  Groups   Name        Std.Dev. Corr  
##  newpid   (Intercept) 1.46           
##           visage      0.30     -0.53 
##  Residual             0.75           
## ---
## number of obs: 1075, groups: newpid, 251
## AIC = 3157.3, DIC = 3116.3
## deviance = 3128.8
\end{verbatim}

\it
Again, no significant effects for \verb|visage|, and a rather large $\hat\tau_0$.

Here's the facet plot\ldots{} \rm

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{coef}\FloatTok{.4} \OtherTok{\textless{}{-}} \FunctionTok{coef}\NormalTok{(lmer}\FloatTok{.4}\NormalTok{)}\SpecialCharTok{$}\NormalTok{newpid[,}\DecValTok{1}\SpecialCharTok{:}\DecValTok{4}\NormalTok{]  }\DocumentationTok{\#\# coef.4[,1] = beta0 + eta0  (random intcpt)}
                                     \DocumentationTok{\#\# coef.4[,2] = beta1 + eta1  (random slope on time)}
                                     \DocumentationTok{\#\# coef.4[,3] = beta2         (coef on treatmnt)}
                                     \DocumentationTok{\#\# coef.4[,4] = beta3         (coef on interaction)}

\NormalTok{tx }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cd4,}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{split}\NormalTok{(treatmnt,newpid), }\ControlFlowTok{function}\NormalTok{(x) x[}\DecValTok{1}\NormalTok{]))  }\DocumentationTok{\#\# one per child}

\NormalTok{params }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{newpid=}\FunctionTok{sort}\NormalTok{(}\FunctionTok{unique}\NormalTok{(cd4}\SpecialCharTok{$}\NormalTok{newpid)),}
                     \AttributeTok{alpha0=}\NormalTok{coef}\FloatTok{.4}\NormalTok{[,}\DecValTok{1}\NormalTok{],}
                     \AttributeTok{alpha1=}\NormalTok{coef}\FloatTok{.4}\NormalTok{[,}\DecValTok{2}\NormalTok{],}
                     \AttributeTok{beta2=}\NormalTok{coef}\FloatTok{.4}\NormalTok{[,}\DecValTok{3}\NormalTok{],}
                     \AttributeTok{beta3=}\NormalTok{coef}\FloatTok{.4}\NormalTok{[,}\DecValTok{4}\NormalTok{],}
\NormalTok{                     tx }
\NormalTok{                     )}

\NormalTok{h }\SpecialCharTok{+} \FunctionTok{geom\_abline}\NormalTok{(}\AttributeTok{data=}\NormalTok{params[params}\SpecialCharTok{$}\NormalTok{newpid}\SpecialCharTok{\textless{}=}\DecValTok{12}\NormalTok{,],}\FunctionTok{aes}\NormalTok{(}\AttributeTok{intercept=}\NormalTok{alpha0 }\SpecialCharTok{+} \FunctionTok{ifelse}\NormalTok{(tx}\SpecialCharTok{==}\DecValTok{2}\NormalTok{,beta2,}\DecValTok{0}\NormalTok{),}
                                                    \AttributeTok{slope=}\NormalTok{alpha1 }\SpecialCharTok{+} \FunctionTok{ifelse}\NormalTok{(tx}\SpecialCharTok{==}\DecValTok{2}\NormalTok{,beta3,}\DecValTok{0}\NormalTok{),}
                                                    \AttributeTok{color=}\NormalTok{tx))}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-12-1.pdf}

\it
It does not appear that \verb|treatmnt| is having a significant effects on the intercepts or the slopes.  We could explore a little further with information criteria:
\rm

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{AIC.ml }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(M) \{}\FunctionTok{AIC}\NormalTok{(}\FunctionTok{update}\NormalTok{(M,}\AttributeTok{REML=}\NormalTok{F))\}}
\NormalTok{BIC.ml }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(M) \{}\FunctionTok{BIC}\NormalTok{(}\FunctionTok{update}\NormalTok{(M,}\AttributeTok{REML=}\NormalTok{F))\}}
\NormalTok{DIC.ml }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(M) \{}\FunctionTok{extractDIC}\NormalTok{(}\FunctionTok{update}\NormalTok{(M,}\AttributeTok{REML=}\NormalTok{F))\}  }\DocumentationTok{\#\# from library(arm)}

\NormalTok{lmer}\FloatTok{.0} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00205385 (tol = 0.002, component 1)
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{res }\OtherTok{\textless{}{-}} \FunctionTok{rbind}\NormalTok{(}\AttributeTok{AIC=}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{list}\NormalTok{(lmer}\FloatTok{.0}\NormalTok{,lmer}\FloatTok{.1}\NormalTok{,lmer}\FloatTok{.2}\NormalTok{,lmer}\FloatTok{.3}\NormalTok{,lmer}\FloatTok{.4}\NormalTok{), AIC.ml),}
             \AttributeTok{BIC=}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{list}\NormalTok{(lmer}\FloatTok{.0}\NormalTok{,lmer}\FloatTok{.1}\NormalTok{,lmer}\FloatTok{.2}\NormalTok{,lmer}\FloatTok{.3}\NormalTok{,lmer}\FloatTok{.4}\NormalTok{), BIC.ml),}
             \AttributeTok{DIC=}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{list}\NormalTok{(lmer}\FloatTok{.0}\NormalTok{,lmer}\FloatTok{.1}\NormalTok{,lmer}\FloatTok{.2}\NormalTok{,lmer}\FloatTok{.3}\NormalTok{,lmer}\FloatTok{.4}\NormalTok{), DIC.ml))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00584617 (tol = 0.002, component 1)
\end{verbatim}

\begin{verbatim}
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00584617 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00584617 (tol = 0.002, component 1)
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{colnames}\NormalTok{(res) }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\StringTok{"lmer.0"}\NormalTok{,}\StringTok{"lmer.1"}\NormalTok{, }\StringTok{"lmer.2"}\NormalTok{, }\StringTok{"lmer.3"}\NormalTok{,}\StringTok{"lmer.4"}\NormalTok{)}

\FunctionTok{t}\NormalTok{(}\FunctionTok{round}\NormalTok{(res,}\DecValTok{2}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##            AIC     BIC     DIC
## lmer.0 3126.56 3156.44 3114.56
## lmer.1 3127.63 3162.49 3113.63
## lmer.2 3142.76 3177.62 3128.76
## lmer.3 3129.21 3169.06 3113.21
## lmer.4 3144.75 3184.60 3128.75
\end{verbatim}

\it
If we restrict our attention to \verb|lmer.1| through \verb|lmer.4|, the winner seems to be \verb|lmer.1|, with model formula \verb:sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1+VISIT|newpid):.  Comparing that with \verb|lmer.0| which just omits \verb|treatmnt|, the simpler \verb|lmer.0| model seems to be favored; all of this is consistent with our earlier findings that treatmt doesn't seem to have a significant effect on either the slope on time, or the intercept, of the MLM's we have built.

\vspace{2em}

\textbf{NOTE (not to be graded):}

What I suspect is going on is that the random part of the intercept,
\(\eta_{0j}\), is accounting for variation between treated and untreated
children, that would otherwise show up in the \verb|treatmnt| effect
coefficient. We could check this out by removing the random part of the
interecept from one of the models and comparing.

Below I have fitted the model with fixed intercept and random slope
(note the \verb|0+visage| term in the lmer call; the 0 removes the
random intercept \(\eta_{0j}\)). As you can see, the coefficient on
\verb|treatmnt| is now (barely) significant, and the model still tracks
the data well.

Comparing all these models is also a good exercise in seeing the
difference between getting a model that does well with prediction, vs
getting a model that might be the ``true'' model. All of the models
we've looked at seem to do well with in-sample prediction (out-of-sample
would also be interesting to look at), but they obviously can't all be
the ``correct/true'' model describing the effect of treatment and time
on \verb|CD4PCT| levels. \rm

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{display}\NormalTok{(lmer}\FloatTok{.5} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ visage }\SpecialCharTok{+}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{0}\SpecialCharTok{+}\NormalTok{visage}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## lmer(formula = sqrt.CD4PCT ~ 1 + visage + treatmnt + (0 + visage | 
##     newpid), data = cd4)
##             coef.est coef.se
## (Intercept)  5.16     0.13  
## visage      -0.20     0.04  
## treatmnt2    0.33     0.14  
## 
## Error terms:
##  Groups   Name   Std.Dev.
##  newpid   visage 0.42    
##  Residual        0.82    
## ---
## number of obs: 1075, groups: newpid, 251
## AIC = 3297.1, DIC = 3267.3
## deviance = 3277.2
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{coef}\FloatTok{.5} \OtherTok{\textless{}{-}} \FunctionTok{coef}\NormalTok{(lmer}\FloatTok{.5}\NormalTok{)}\SpecialCharTok{$}\NormalTok{newpid[,}\DecValTok{1}\SpecialCharTok{:}\DecValTok{3}\NormalTok{]  }\DocumentationTok{\#\# coef.5[,1] = beta0         (fixed intcpt)}
                                     \DocumentationTok{\#\# coef.5[,2] = beta1 + eta1  (random slope on time)}
                                     \DocumentationTok{\#\# coef.5[,3] = beta2         (coef on treatmnt)}
                                  
\NormalTok{tx }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cd4,}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{split}\NormalTok{(treatmnt,newpid), }\ControlFlowTok{function}\NormalTok{(x) x[}\DecValTok{1}\NormalTok{]))  }\DocumentationTok{\#\# one per child}

\NormalTok{params }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{newpid=}\FunctionTok{sort}\NormalTok{(}\FunctionTok{unique}\NormalTok{(cd4}\SpecialCharTok{$}\NormalTok{newpid)),}
                     \AttributeTok{beta0=}\NormalTok{coef}\FloatTok{.5}\NormalTok{[,}\DecValTok{1}\NormalTok{],}
                     \AttributeTok{alpha1=}\NormalTok{coef}\FloatTok{.5}\NormalTok{[,}\DecValTok{2}\NormalTok{],}
                     \AttributeTok{beta2=}\NormalTok{coef}\FloatTok{.5}\NormalTok{[,}\DecValTok{3}\NormalTok{],}
\NormalTok{                     tx }
\NormalTok{                     )}

\NormalTok{h }\SpecialCharTok{+} \FunctionTok{geom\_abline}\NormalTok{(}\AttributeTok{data=}\NormalTok{params[params}\SpecialCharTok{$}\NormalTok{newpid}\SpecialCharTok{\textless{}=}\DecValTok{12}\NormalTok{,],}\FunctionTok{aes}\NormalTok{(}\AttributeTok{intercept=}\NormalTok{beta0 }\SpecialCharTok{+} \FunctionTok{ifelse}\NormalTok{(tx}\SpecialCharTok{==}\DecValTok{2}\NormalTok{,beta2,}\DecValTok{0}\NormalTok{),}
                                                    \AttributeTok{slope=}\NormalTok{alpha1,}
                                                    \AttributeTok{color=}\NormalTok{tx))}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-14-1.pdf}

\hypertarget{problem-2.}{%
\section{Problem 2.}\label{problem-2.}}

Now we will consider in more detail the model below for the CD4 data
set:

\begin{verbatim}
lmer.1 <- lmer(sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (1 + VISIT | newpid),
               data=cd4)
\end{verbatim}

\it 

Bring together all the residuals (except random effects residuals) that
we have talked about, using library(HLMdiag):

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{cd4 }\OtherTok{\textless{}{-}} \FunctionTok{read.csv}\NormalTok{(}\StringTok{"allvar.csv"}\NormalTok{,}\AttributeTok{header=}\NormalTok{T)}

\DocumentationTok{\#\# It\textquotesingle{}s worth doing a little exploration of the data here,}
\DocumentationTok{\#\# especially to look at missing values.}

\FunctionTok{apply}\NormalTok{(cd4,}\DecValTok{2}\NormalTok{,}\ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{mean}\NormalTok{(}\FunctionTok{is.na}\NormalTok{(x)))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##       VISIT      newpid       VDATE      CD4PCT         arv      visage 
## 0.000000000 0.000000000 0.000000000 0.142743222 0.103668262 0.109250399 
##    treatmnt      CD4CNT     baseage 
## 0.000000000 0.146730463 0.007177033
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# we only need the variables newpid, VISIT, visage, CD4PCT and treatmnt}
\DocumentationTok{\#\# so let\textquotesingle{}s get rid of some (but not all) missing data problems}
\DocumentationTok{\#\# by deleting the other variables}

\NormalTok{cd4 }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cd4,}\FunctionTok{data.frame}\NormalTok{(}\AttributeTok{newpid=}\NormalTok{newpid,}
                           \AttributeTok{VISIT=}\NormalTok{VISIT, }
                           \AttributeTok{visage=}\NormalTok{visage, }
                           \AttributeTok{CD4PCT=}\NormalTok{CD4PCT, }
                           \AttributeTok{treatmnt=}\NormalTok{treatmnt))}

\DocumentationTok{\#\# and since we still have some missing data, we will just }
\DocumentationTok{\#\# delete the rows that continue to have NA\textquotesingle{}s (not a great}
\DocumentationTok{\#\# practice in general, but good enough for this exercise...)}
\NormalTok{cd4 }\OtherTok{\textless{}{-}}\NormalTok{ cd4[}\SpecialCharTok{!}\FunctionTok{apply}\NormalTok{(cd4,}\DecValTok{1}\NormalTok{,}\ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{any}\NormalTok{(}\FunctionTok{is.na}\NormalTok{(x))),]}

\DocumentationTok{\#\# set up the sqrt of the response variable...}
\NormalTok{cd4}\SpecialCharTok{$}\NormalTok{sqrt.CD4PCT }\OtherTok{\textless{}{-}} \FunctionTok{sqrt}\NormalTok{(cd4}\SpecialCharTok{$}\NormalTok{CD4PCT)}

\DocumentationTok{\#\# and for coloring ggplot elements it will better if treatmnt is a factor...}
\NormalTok{cd4}\SpecialCharTok{$}\NormalTok{treatmnt }\OtherTok{\textless{}{-}} \FunctionTok{as.factor}\NormalTok{(cd4}\SpecialCharTok{$}\NormalTok{treatmnt)}

\DocumentationTok{\#\# select the first 12 kids...}
\NormalTok{first}\FloatTok{.12} \OtherTok{\textless{}{-}}\NormalTok{ (cd4}\SpecialCharTok{$}\NormalTok{newpid }\SpecialCharTok{\textless{}=} \DecValTok{12}\NormalTok{)}

\DocumentationTok{\#\# Fit the base model for this exercise...}
\FunctionTok{display}\NormalTok{(lmer}\FloatTok{.1} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{*}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{|}\NormalTok{ newpid),}
               \AttributeTok{data=}\NormalTok{cd4))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## lmer(formula = sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (1 + VISIT | 
##     newpid), data = cd4)
##                 coef.est coef.se
## (Intercept)      4.71     0.13  
## VISIT           -0.03     0.01  
## treatmnt2        0.14     0.19  
## VISIT:treatmnt2  0.01     0.01  
## 
## Error terms:
##  Groups   Name        Std.Dev. Corr  
##  newpid   (Intercept) 1.40           
##           VISIT       0.05     -0.10 
##  Residual             0.72           
## ---
## number of obs: 1075, groups: newpid, 251
## AIC = 3149.4, DIC = 3093
## deviance = 3113.2
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{r}\FloatTok{.1} \OtherTok{\textless{}{-}} \FunctionTok{hlm\_resid}\NormalTok{(lmer}\FloatTok{.1}\NormalTok{,}\AttributeTok{level=}\DecValTok{1}\NormalTok{,}\AttributeTok{include.ls=}\NormalTok{F)}
\NormalTok{r}\FloatTok{.1}\NormalTok{s }\OtherTok{\textless{}{-}} \FunctionTok{hlm\_resid}\NormalTok{(lmer}\FloatTok{.1}\NormalTok{,}\AttributeTok{level=}\DecValTok{1}\NormalTok{,}\AttributeTok{include.ls=}\NormalTok{F,}\AttributeTok{standardize=}\NormalTok{T)}
\NormalTok{r}\FloatTok{.2} \OtherTok{\textless{}{-}} \FunctionTok{hlm\_resid}\NormalTok{(lmer}\FloatTok{.1}\NormalTok{,}\AttributeTok{level=}\StringTok{"newpid"}\NormalTok{,}\AttributeTok{include.ls=}\NormalTok{F)}
\NormalTok{r}\FloatTok{.2}\NormalTok{s }\OtherTok{\textless{}{-}} \FunctionTok{hlm\_resid}\NormalTok{(lmer}\FloatTok{.1}\NormalTok{,}\AttributeTok{level=}\StringTok{"newpid"}\NormalTok{,}\AttributeTok{include.ls=}\NormalTok{F,}\AttributeTok{standardize=}\NormalTok{T)}
\FunctionTok{names}\NormalTok{(r}\FloatTok{.1}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] "id"          "sqrt.CD4PCT" "VISIT"       "treatmnt"    "newpid"     
## [6] ".resid"      ".fitted"     ".mar.resid"  ".mar.fitted"
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{names}\NormalTok{(r}\FloatTok{.1}\NormalTok{s)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] "id"              "sqrt.CD4PCT"     "VISIT"           "treatmnt"       
## [5] "newpid"          ".std.resid"      ".fitted"         ".chol.mar.resid"
## [9] ".mar.fitted"
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{names}\NormalTok{(r}\FloatTok{.2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] "newpid"           ".ranef.intercept" ".ranef.visit"
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{names}\NormalTok{(r}\FloatTok{.2}\NormalTok{s)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] "newpid"               ".std.ranef.intercept" ".std.ranef.visit"
\end{verbatim}

\rm

\hypertarget{problem-2a.}{%
\subsection{Problem 2(a).}\label{problem-2a.}}

Make a facets plot of the marginal residuals, as a function of the
marginal fitted values (use \texttt{scales="free\_x"} if needed to make
the plot legible). Explain in a sentence or two why a facets plot is not
very useful for assessing model fit for this problem, whether we look at
the first 12 children, or all 251 children).

\it

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{ggplot}\NormalTok{(r}\FloatTok{.1}\NormalTok{[first}\FloatTok{.12}\NormalTok{,],}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{.mar.fitted,}\AttributeTok{y=}\NormalTok{.mar.resid)) }\SpecialCharTok{+}
  \FunctionTok{facet\_wrap}\NormalTok{( }\SpecialCharTok{\textasciitilde{}}\NormalTok{ newpid,}\AttributeTok{scales=}\StringTok{"free\_x"}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{geom\_point}\NormalTok{(}\FunctionTok{aes}\NormalTok{(}\AttributeTok{color=}\NormalTok{treatmnt)) }\SpecialCharTok{+}
  \FunctionTok{geom\_abline}\NormalTok{(}\AttributeTok{intercept=}\DecValTok{0}\NormalTok{,}\AttributeTok{slope=}\DecValTok{0}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-21-1.pdf}

There are very few data points per child, and so the residual plot
within facet for each child contains almost no information about fit.
(Having said that, there are some hints of misfit, to the extent that
several individual children's residuals are entirely above, or entirely
below, the zero line.) \rm

\hypertarget{problem-2b.}{%
\subsection{Problem 2(b).}\label{problem-2b.}}

Make an ungrouped (that is, no facets) scatter plot of marginal
residuals as a function of marginal fitted values, using the full data
set (not just the first 12 children). Color the points for
\texttt{treatmnt=1} kids and treatmnt=2 kids with different colors.
Overlay a smooth fit (\texttt{geom\_smooth} is the easist to use here).
Explain, in a couple of sentences (optionally with some math):

\begin{itemize}
\item What is causing the dominant structure in this plot, and why that dominant structure is essentially irrelevant for checking the relationship between \texttt{sqrt.CD4PCT} and \texttt{VISIT}; and  
\item What in this plot makes you happy or unhappy about having a linear relationship between  \texttt{sqrt.CD4PCT} and \texttt{VISIT} in the model.
\end{itemize}

\it

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{ggplot}\NormalTok{(r}\FloatTok{.1}\NormalTok{,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{.mar.fitted,}\AttributeTok{y=}\NormalTok{.mar.resid)) }\SpecialCharTok{+}
  \FunctionTok{geom\_point}\NormalTok{(}\FunctionTok{aes}\NormalTok{(}\AttributeTok{color=}\NormalTok{treatmnt)) }\SpecialCharTok{+}
  \FunctionTok{geom\_smooth}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-22-1.pdf}

\rm

Explain, in a couple of sentences (optionally with some math):

\begin{itemize}
\item What is causing the dominant structure in this plot, and why that dominant structure is essentially irrelevant for checking the relationship between \texttt{sqrt.CD4PCT} and \texttt{VISIT};

\it

The dominant structure in the data is the way the residuals are grouped horizontally into 14 groups---7 groups for \texttt{treatmnt=1} and 7 for \texttt{treatmnt=2}.  They are caused by the fact that our main predictor variable is \texttt{VISIT}, which only takes on 7 values. The spacing of the groups is not the same between \texttt{treatmnt=1} and  \texttt{treatmnt=2} because the model contains a \texttt{treatmnt} $\times$ \texttt{VISIT} interaction, so that the slope on \texttt{VISIT} is different for children in the \texttt{treatmnt=1} and  \texttt{treatmnt=2} groups.  

This does not affect our assessment of the fit of the model, since to check the model fit we are interested in vertical patterns (trends, curves, outliers, changing variance, etc.) in the residuals, not horizontal (grouping) patterns.

\rm

and  
\item What in this plot makes you happy or unhappy about having a linear relationship between  \texttt{sqrt.CD4PCT} and \texttt{VISIT} in the model.

\it

The plot actually looks pretty wonderful.  The variance of the residuals looks pretty much constant, as a function of the fitted values, and the trend modeled with the smooth really does look like a horizontal line at 0.

\rm
\end{itemize}

\hypertarget{problem-2c.}{%
\subsection{Problem 2(c).}\label{problem-2c.}}

Make an ungrouped (that is, no facets) scatter plot of conditional
residuals as a function of conditional fitted values, using the full
data set (not just the first 12 children). Color the points for
\texttt{treatmnt=1} kids and \texttt{treatmnt=2} kids with different
colors. Overlay a smooth fit (\texttt{geom\_smooth} is the easist to use
here). Explain, in a couple of sentences (optionally with some math):

\begin{itemize}
\item Why the dominant structure in the marginal residuals is not also present in this plot of conditional residuals
\item What might be causing the trend you see in this plot to be different from the trend in the plot of the marginal residuals.
\end{itemize}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{ggplot}\NormalTok{(r}\FloatTok{.1}\NormalTok{,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{.fitted,}\AttributeTok{y=}\NormalTok{.resid)) }\SpecialCharTok{+}
  \FunctionTok{geom\_point}\NormalTok{(}\FunctionTok{aes}\NormalTok{(}\AttributeTok{color=}\NormalTok{treatmnt)) }\SpecialCharTok{+}
  \FunctionTok{geom\_smooth}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-23-1.pdf}

Explain, in a couple of sentences (optionally with some math):

\begin{itemize}
\item Why the dominant structure in the marginal residuals is not also present in this plot of conditional residuals

\it

The marginal residuals plot shows
\rm
\begin{eqnarray*}
\mbox{.mar.resid}_i & = & \mbox{sqrt.CD4PCT}_i - \hat\beta_0 - \hat\beta_1\cdot\mbox{VISIT}_i - 
\hat\beta_3\cdot\mbox{treatmnt}_{j[i]} - \hat\beta_4\cdot\mbox{VISIT}_i\cdot\mbox{treatmnt}_{j[i]}\\
& \mbox{vs.}& \\
\mbox{.mar.fitted}_i & = & \hat\beta_0 + \hat\beta_1\cdot\mbox{VISIT}_i + 
\hat\beta_3\cdot\mbox{treatmnt}_{j[i]} + \hat\beta_4\cdot\mbox{VISIT}_i\cdot\mbox{treatmnt}_{j[i]}
\end{eqnarray*}
\it
and because \texttt{VISIT} and \texttt{treatmnt} are both discrete, the fitted values will only take on a few discrete values, causing the horizontal grouping in the marginal plot.

The conditional residuals plot shows
\rm
\begin{eqnarray*}
\mbox{.resid}_i & = & \mbox{sqrt.CD4PCT}_i - \hat\beta_0 - \hat\beta_1\cdot\mbox{VISIT}_i - 
\hat\beta_3\cdot\mbox{treatmnt}_{j[i]} - \hat\beta_4\cdot\mbox{VISIT}_i\cdot\mbox{treatmnt}_{j[i]}\\
&& - \eta_{0j[i]} - \eta_{1j[i]}\cdot\mbox{VISIT}_i \\
& \mbox{vs.}  & \\
\mbox{.fitted}_i & = & \hat\beta_0 + \hat\beta_1\cdot\mbox{VISIT}_i +
\hat\beta_3\cdot\mbox{treatmnt}_{j[i]} + \hat\beta_4\cdot\mbox{VISIT}_i\cdot\mbox{treatmnt}_{j[i]}\\
&& + \eta_{0j[i]} + \eta_{1j[i]}\cdot\mbox{VISIT}_i
\end{eqnarray*}
\it
The additional continuous (because the $\eta$'s are continuous) terms $\eta_{0j[i]} + \eta_{1j[i]}\cdot\mbox{VISIT}_i$ essentially "jitter" the fitted values so they no longer show the horizonal grouping structure of the marginal residuals.


\rm

\item What might be causing the trend you see in this plot to be different from the trend in the plot of the marginal residuals.

\it

The trend shows that when the conditional fitted value $\hat y_{cond}$ is larger, it tends to underpredict $y=\mbox{sqrt.CD4PCT}$ (positive residual), and when $\hat y_{cond}$ is smaller, it tends to overpredict $y$ (negative residual).

The only difference between the two plots is the presence of the $\eta_{0j[i]} + \eta_{1j[i]}\cdot\mbox{VISIT}_i$ terms.

In ordinary regression,we know that residuals and fitted values will be uncorrelated: that is why we don't see overall increasing or overall decreasing trends in the marginal residual plots (as far as the math is concerned, marginal residuals and marginal fitted values behave exactly like ordinary regression).  The fact that conditional residuals and conditional fitted values are correlated is due to the values we use for the $\eta$'s: $\hat\eta=E[\eta|\mbox{the data}]$, and "the data" includes $y$.  Since, therefore, the conditional fitted values depend on $y$, it's not surprising that the residuals and fitted values are correlated.

The plots below explore the association between the $\hat\eta$'s and $y$.  We can see that, although there's not much association between $\hat\eta_{1j}$ and $y$, there is quite a strong association between $\hat\eta_{0j}$ and $y$. 

\rm

\end{itemize}

\vspace*{1em}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{y.grouped }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cd4,}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{split}\NormalTok{(sqrt.CD4PCT,newpid),mean))}
\NormalTok{tx.grouped }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cd4,}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{split}\NormalTok{(treatmnt,newpid),}\ControlFlowTok{function}\NormalTok{(x) x[}\DecValTok{1}\NormalTok{]))}

\NormalTok{g1 }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(r}\FloatTok{.2}\NormalTok{,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{y.grouped,}\AttributeTok{y=}\NormalTok{.ranef.intercept)) }\SpecialCharTok{+}
  \FunctionTok{geom\_point}\NormalTok{(}\FunctionTok{aes}\NormalTok{(}\AttributeTok{color=}\NormalTok{tx.grouped)) }\SpecialCharTok{+}
  \FunctionTok{xlab}\NormalTok{(}\StringTok{"mean(sqrt.CD4PCT) within child"}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylab}\NormalTok{(}\StringTok{"eta0"}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{theme}\NormalTok{(}\AttributeTok{legend.position=}\StringTok{"none"}\NormalTok{)}

\NormalTok{g2 }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(r}\FloatTok{.2}\NormalTok{,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{y.grouped,}\AttributeTok{y=}\NormalTok{.ranef.visit)) }\SpecialCharTok{+}
  \FunctionTok{geom\_point}\NormalTok{(}\FunctionTok{aes}\NormalTok{(}\AttributeTok{color=}\NormalTok{tx.grouped)) }\SpecialCharTok{+}
  \FunctionTok{xlab}\NormalTok{(}\StringTok{"mean(sqrt.CD4PCT) within child"}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylab}\NormalTok{(}\StringTok{"eta1"}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{theme}\NormalTok{(}\AttributeTok{legend.position=}\StringTok{"none"}\NormalTok{)}

\FunctionTok{grid.arrange}\NormalTok{(g1,g2,}\AttributeTok{ncol=}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-24-1.pdf}

\hypertarget{problem-2d.}{%
\subsection{Problem 2(d).}\label{problem-2d.}}

Use standardized residuals and standardized random effects estimates to
assess the normality of \(\epsilon_i\), \(\eta_{0j}\) and \(\eta_{1j}\)
in the fitted model, and to check for any outliers. Include qq plots for
each, and accompany each plot with a sentence or two describing what is
good or bad in that plot.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{g1 }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(r}\FloatTok{.1}\NormalTok{s,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{sample=}\NormalTok{.std.resid)) }\SpecialCharTok{+} 
  \FunctionTok{geom\_qq}\NormalTok{() }\SpecialCharTok{+}
  \FunctionTok{geom\_qq\_line}\NormalTok{() }\SpecialCharTok{+}
  \FunctionTok{xlim}\NormalTok{(}\SpecialCharTok{{-}}\FloatTok{3.5}\NormalTok{,}\FloatTok{3.5}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylim}\NormalTok{(}\SpecialCharTok{{-}}\DecValTok{8}\NormalTok{,}\DecValTok{8}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{xlab}\NormalTok{(}\StringTok{"Normal Quantiles"}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylab}\NormalTok{(}\StringTok{"Data Quantiles"}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{ggtitle}\NormalTok{(}\StringTok{"Standardized conditional residuals"}\NormalTok{)}

\NormalTok{g2 }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(r}\FloatTok{.1}\NormalTok{s,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{sample=}\NormalTok{.chol.mar.resid)) }\SpecialCharTok{+} 
  \FunctionTok{geom\_qq}\NormalTok{() }\SpecialCharTok{+}
  \FunctionTok{geom\_qq\_line}\NormalTok{() }\SpecialCharTok{+}
  \FunctionTok{xlim}\NormalTok{(}\SpecialCharTok{{-}}\FloatTok{3.5}\NormalTok{,}\FloatTok{3.5}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylim}\NormalTok{(}\SpecialCharTok{{-}}\DecValTok{8}\NormalTok{,}\DecValTok{8}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{xlab}\NormalTok{(}\StringTok{"Normal Quantiles"}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylab}\NormalTok{(}\StringTok{"Data Quantiles"}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{ggtitle}\NormalTok{(}\StringTok{"Cholesky marginal residuals"}\NormalTok{)}

\NormalTok{g3 }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(r}\FloatTok{.2}\NormalTok{s,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{sample=}\NormalTok{.std.ranef.intercept)) }\SpecialCharTok{+} 
  \FunctionTok{geom\_qq}\NormalTok{() }\SpecialCharTok{+}
  \FunctionTok{geom\_qq\_line}\NormalTok{() }\SpecialCharTok{+}
  \FunctionTok{xlim}\NormalTok{(}\SpecialCharTok{{-}}\FloatTok{3.5}\NormalTok{,}\FloatTok{3.5}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylim}\NormalTok{(}\SpecialCharTok{{-}}\DecValTok{8}\NormalTok{,}\DecValTok{8}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{xlab}\NormalTok{(}\StringTok{"Normal Quantiles"}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylab}\NormalTok{(}\StringTok{"Data Quantiles"}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{ggtitle}\NormalTok{(}\StringTok{"Standardized eta0"}\NormalTok{)}

\NormalTok{g4 }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(r}\FloatTok{.2}\NormalTok{s,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{sample=}\NormalTok{.std.ranef.visit)) }\SpecialCharTok{+} 
  \FunctionTok{geom\_qq}\NormalTok{() }\SpecialCharTok{+}
  \FunctionTok{geom\_qq\_line}\NormalTok{() }\SpecialCharTok{+}
  \FunctionTok{xlim}\NormalTok{(}\SpecialCharTok{{-}}\FloatTok{3.5}\NormalTok{,}\FloatTok{3.5}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylim}\NormalTok{(}\SpecialCharTok{{-}}\DecValTok{8}\NormalTok{,}\DecValTok{8}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{xlab}\NormalTok{(}\StringTok{"Normal Quantiles"}\NormalTok{) }\SpecialCharTok{+} 
  \FunctionTok{ylab}\NormalTok{(}\StringTok{"Data Quantiles"}\NormalTok{) }\SpecialCharTok{+}
  \FunctionTok{ggtitle}\NormalTok{(}\StringTok{"Standardized eta1"}\NormalTok{)}

\FunctionTok{grid.arrange}\NormalTok{(g1,g2,g3,g4,}\AttributeTok{ncol=}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning: Removed 1 rows containing missing values (geom_point).
\end{verbatim}

\includegraphics{hw08-solutions_files/figure-latex/unnamed-chunk-25-1.pdf}
\it

The ``middle'' part of all four plots seem to follow the normal
distribution fairly well, but the tails fail in various ways:

\begin{itemize}

\item The standardized conditional and marginal residuals both have longer tails than the Normal distribution.  There seem to be two low outliers and two high outliers among these residuals.

\item The Cholesky marginal residual seem a little more extreme that the standardized conditional residuals, which make sense since $X\beta + Z\eta$ should do a better job predicting $y$ than $X\beta$ alone.   There seem to be two low outliers and two high outliers among these residuals.  It would be useful to know if these are the same data points as in the first QQ plot.

\item The standardized $\hat\eta_0$'s have a long left tail and slighltly short right tail.  In fact the left tail for $\eta_0$ is the most extreme of all four plots.  There don't appear to be any clear outliers among the $\hat\eta_0$'s.

\item The distribution of $\hat\eta_1$ is closest to Normal, with a right tail that is only a little longer than the Normal, and a left tail that is the least extreme of all four plots.  There may be one low outlier among the $\hat\eta_1$'s.

\end{itemize}

\rm

\hypertarget{problem-3.}{%
\section{Problem 3.}\label{problem-3.}}

Continuing with the CD4 data\ldots

\begin{quote}\em

{\bf Note:} Because I am not covering cAIC very much in class I have asked you not to compare models with cAIC (that is, no credit lost if you do not consider cAIC).  However, since it is on the assignment, I wanted to (a) make some notes on it here; and (b) show what the results using cAIC would be in the solutions below.

Here are notes on cAIC:
\begin{itemize}

\item Earlier in the semester we worked with a quantity called ``CAIC''.  Since this quantity is called ``cAIC'' you might think these are the same.  They are actually not at all the same!
\begin{itemize}
\item CAIC is only appropriate for regular lm() models.  It is a version of AIC that has been adjusted for small sample size:
\begin{eqnarray*}
AIC & = & -2logLik(M) + 2 (k+1) ~~~ \mbox{, where $k=df=p+1$.} \\
CAIC & = & -2logLik(M) + 2\frac{n+2}{n-k}(k+1)
\end{eqnarray*}
\item cAIC is designed for lmer() models. Like DIC, it is a version of AIC that has been adjusted for reduced degrees of freedom due to the presence of random effects.  The mathematics is a bit much for the amount of time/space that we have in this class, but here's the overview from 10,000 feet:
\begin{description}

\item[DIC:] $DIC = -2logLik(M) + 2 k_{eff}$
\begin{itemize}
\item $LogLik(M)$ based on marginal model $f(\underline{Y}|\underline{\beta}, \underline{\omega}, \sigma^2)$
\item $k_{eff}$ is estimated from the curvature of the likelihood $f(\underline{Y}|\underline{\beta}, \underline{\omega}, \sigma^2)$, which is driven by the size of the $\tau^2$'s
\end{itemize}

\item[cAIC:] $cAIC = -2logLik(M) + 2 k_{eff}$
\begin{itemize}
\item $logLik(M)$ based on the conditional model $f(\underline{Y}|\underline{\eta}, \underline{\beta}, \underline{\omega}, \sigma^2)$
\item $k_{eff}$ is a different ``model curvature'' estimate
\end{itemize}

\end{description}

\end{itemize}

\item You can use the cAIC() function from library(cAIC4) to calculate cAIC.  

\item Since cAIC and DIC are both adjustments to AIC to account for random $\eta$'s, when there are no $\eta$'s in the problem, both cAIC and DIC reduce to just AIC.

\end{itemize}

In the solutions below I have just integrated work with cAIC with all the other work on this problem.  

\end{quote}

\hypertarget{problem-3a.}{%
\subsection{Problem 3(a).}\label{problem-3a.}}

Make a table giving values of AIC, BIC, DIC, and cAIC:

\begin{verbatim}
sqrt.CD4PCT ~ 1 + visage +            (1+visage|newpid)
sqrt.CD4PCT ~ 1 + visage + treatmnt + (1+visage|newpid)
sqrt.CD4PCT ~ 1 + visage * treatmnt + (1+visage|newpid)
sqrt.CD4PCT ~ 1 + VISIT  +            (1+VISIT|newpid)
sqrt.CD4PCT ~ 1 + VISIT  + treatmnt + (1+VISIT|newpid)
sqrt.CD4PCT ~ 1 + VISIT  * treatmnt + (1+VISIT|newpid)
\end{verbatim}

Comment briefly on any similarities or differences in how the different
criteria choose fixed effects.

\it

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{lmer}\FloatTok{.1} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ visage }\SpecialCharTok{+}\NormalTok{            (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{visage}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\NormalTok{lmer}\FloatTok{.2} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ visage }\SpecialCharTok{+}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{visage}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\NormalTok{lmer}\FloatTok{.3} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ visage }\SpecialCharTok{*}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{visage}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\NormalTok{lmer}\FloatTok{.4} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT  }\SpecialCharTok{+}\NormalTok{            (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\NormalTok{lmer}\FloatTok{.5} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT  }\SpecialCharTok{+}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00584617 (tol = 0.002, component 1)
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{lmer}\FloatTok{.6} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT  }\SpecialCharTok{*}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}

\NormalTok{DIC }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(M) \{}
         \ControlFlowTok{if}\NormalTok{(}\FunctionTok{class}\NormalTok{(M)}\SpecialCharTok{==}\StringTok{"lm"}\NormalTok{) \{ x }\OtherTok{\textless{}{-}} \FunctionTok{AIC}\NormalTok{(M) \}}
         \ControlFlowTok{else}\NormalTok{ \{ x }\OtherTok{\textless{}{-}} \FunctionTok{unname}\NormalTok{(}\FunctionTok{extractDIC}\NormalTok{(M)) \}}
         \FunctionTok{return}\NormalTok{(x)}
\NormalTok{       \}}
  
\NormalTok{CAIC }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(M) \{}\FunctionTok{cAIC}\NormalTok{(M)}\SpecialCharTok{$}\NormalTok{caic\}}

\NormalTok{IC }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(M) \{}
\NormalTok{    x }\OtherTok{\textless{}{-}}\FunctionTok{c}\NormalTok{(}\AttributeTok{AIC=}\FunctionTok{AIC}\NormalTok{(M),}\AttributeTok{BIC=}\FunctionTok{BIC}\NormalTok{(M),}\AttributeTok{DIC=}\FunctionTok{DIC}\NormalTok{(M),}\AttributeTok{cAIC=}\FunctionTok{CAIC}\NormalTok{(M))}
    \FunctionTok{return}\NormalTok{(x)}
\NormalTok{\}}

\NormalTok{IC.table }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(...) \{}
\NormalTok{    Mlist }\OtherTok{\textless{}{-}} \FunctionTok{list}\NormalTok{(...)}
\NormalTok{    x }\OtherTok{\textless{}{-}} \FunctionTok{suppressWarnings}\NormalTok{(}\FunctionTok{lapply}\NormalTok{(Mlist,IC))}
\NormalTok{    x }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\FunctionTok{matrix}\NormalTok{(}\FunctionTok{unlist}\NormalTok{(x),}\AttributeTok{ncol=}\DecValTok{4}\NormalTok{,}\AttributeTok{byrow=}\NormalTok{T))}
    \FunctionTok{names}\NormalTok{(x) }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\StringTok{"AIC"}\NormalTok{,}\StringTok{"BIC"}\NormalTok{,}\StringTok{"DIC"}\NormalTok{,}\StringTok{"cAIC"}\NormalTok{)}
\NormalTok{    models }\OtherTok{\textless{}{-}} \FunctionTok{sapply}\NormalTok{(Mlist,formula)}
\NormalTok{    models }\OtherTok{\textless{}{-}} \FunctionTok{abbreviate}\NormalTok{(}\FunctionTok{substr}\NormalTok{(models,}\DecValTok{14}\NormalTok{,}\DecValTok{200}\NormalTok{),}\DecValTok{30}\NormalTok{)}
    \FunctionTok{rownames}\NormalTok{(x) }\OtherTok{\textless{}{-}}\NormalTok{ models}
    \FunctionTok{return}\NormalTok{(x)}
\NormalTok{\}}

\FunctionTok{IC.table}\NormalTok{(lmer}\FloatTok{.1}\NormalTok{,lmer}\FloatTok{.2}\NormalTok{,lmer}\FloatTok{.3}\NormalTok{,lmer}\FloatTok{.4}\NormalTok{,lmer}\FloatTok{.5}\NormalTok{,lmer}\FloatTok{.6}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                                     AIC      BIC      DIC     cAIC
## 1+visage+(1+visage|newpid)     3142.447 3172.327 3130.447 2697.868
## 1+visag+treatmnt+(1+visg|nwpd) 3142.755 3177.616 3128.755 2698.827
## 1+visag*treatmnt+(1+visg|nwpd) 3144.755 3184.595 3128.755 2699.562
## 1+VISIT+(1+VISIT|newpid)       3126.563 3156.443 3114.563 2642.798
## 1+VISIT+tretmnt+(1+VISIT|nwpd) 3127.625 3162.486 3113.625 2643.307
## 1+VISIT*tretmnt+(1+VISIT|nwpd) 3129.215 3169.056 3113.215 2644.733
\end{verbatim}

\rm

Comment briefly on any similarities or differences in how the different
criteria choose fixed effects.

\it

Here's an overall sum of the best (minimizing the criterion) models for
each criterion, with second and third place winners as well:

\begin{description}

\item[AIC:] Best model is lmer.4, with lmer.5 a close second.

\item[BIC:] Best model is lmer.4, with lmer.5 a substantially less close second.

\item[DIC:] Best model is lmer.6, with lmer.4 and lmer.5 almost indistinguishably close.

\item[cAIC:] Best model is lmer.4, with lmer.5 and lmer.6 almost indistinguishably close.

\end{description}

All the criteria prefer having \texttt{VISIT} as a predictor rather than
\texttt{visage}. (This kind of surprises me because visage has more
physical meaning, but so be it\ldots) Consistent with our findings on
earlier hw that \text{treatmnt} did not seem like an important
predictor, \texttt{lmer.4} (without \texttt{treatmnt} in the model) is
first or second with all the criteria, with \texttt{lmer.5} (main effect
for \texttt{treatmt}, but no interaction) a strong second or third with
all the criteria.

\rm

\hypertarget{problem-3b.}{%
\subsection{Problem 3(b).}\label{problem-3b.}}

Now let's look at the random effects in the model

\begin{verbatim}
sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1+VISIT|newpid)
\end{verbatim}

Again, make a table giving values of AIC, BIC, DIC, and cAIC for the
following models:

\begin{verbatim}
sqrt.CD4PCT ~ 1 + VISIT + treatmnt 
sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1|newpid)
sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (0+VISIT|newpid)
sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1+VISIT|newpid)
\end{verbatim}

\it

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{lm}\FloatTok{.7} \OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{+}\NormalTok{ treatmnt,}\AttributeTok{data=}\NormalTok{cd4)}
\NormalTok{lmer}\FloatTok{.8} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{+}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\NormalTok{lmer}\FloatTok{.9} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{+}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{0}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\NormalTok{lmer}\FloatTok{.10} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{+}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00584617 (tol = 0.002, component 1)
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{IC.table}\NormalTok{(lm}\FloatTok{.7}\NormalTok{,lmer}\FloatTok{.8}\NormalTok{,lmer}\FloatTok{.9}\NormalTok{,lmer}\FloatTok{.10}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                                     AIC      BIC      DIC     cAIC
## 1 + VISIT + treatmnt           4009.810 4029.730 4009.810 4009.810
## 1+VISIT+treatmnt+(1|newpid)    3153.879 3178.779 3143.879 2728.147
## 1+VISIT+tretmnt+(0+VISIT|nwpd) 3688.619 3713.520 3678.619 3490.226
## 1+VISIT+tretmnt+(1+VISIT|nwpd) 3127.625 3162.486 3113.625 2643.307
\end{verbatim}

\rm

You'll fit the first model with lm(), and the others with lmer().
Comment briefly on any similarities or differences in how the different
criteria choose random effects. (Note that only DIC and cAIC have a
strong theoretical justification here).

\it

There's a convergence warning for one of the models, which you can (and
should) pursue using the methods in the R notes from lecture.

All four information criteria strongly favor at least a random intercept
(\verb|lm.7| and \verb|lmer.9| are very strongly disfavored). Among the
remaining two models, the model with random slope and random intercept
is strongly favored by all the criteria.

\rm

\hypertarget{problem-3c.}{%
\subsection{Problem 3(c).}\label{problem-3c.}}

Repeat part (b) but with the interaction \verb|VISIT * treatmnt| in each
model instead of just the main effects \verb|VISIT + treatmnt|.

\it

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{lm}\FloatTok{.11} \OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{*}\NormalTok{ treatmnt,}\AttributeTok{data=}\NormalTok{cd4)}
\NormalTok{lmer}\FloatTok{.12} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{*}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\NormalTok{lmer}\FloatTok{.13} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{*}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{0}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}
\NormalTok{lmer}\FloatTok{.14} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(sqrt.CD4PCT }\SpecialCharTok{\textasciitilde{}} \DecValTok{1} \SpecialCharTok{+}\NormalTok{ VISIT }\SpecialCharTok{*}\NormalTok{ treatmnt }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{+}\NormalTok{VISIT}\SpecialCharTok{|}\NormalTok{newpid),}\AttributeTok{data=}\NormalTok{cd4,}\AttributeTok{REML=}\NormalTok{F)}

\FunctionTok{IC.table}\NormalTok{(lm}\FloatTok{.11}\NormalTok{,lmer}\FloatTok{.12}\NormalTok{,lmer}\FloatTok{.13}\NormalTok{,lmer}\FloatTok{.14}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                                     AIC      BIC      DIC     cAIC
## 1 + VISIT * treatmnt           4010.343 4035.244 4010.343 4010.343
## 1+VISIT*treatmnt+(1|newpid)    3155.181 3185.062 3143.181 2729.893
## 1+VISIT*tretmnt+(0+VISIT|nwpd) 3690.150 3720.030 3678.150 3491.763
## 1+VISIT*tretmnt+(1+VISIT|nwpd) 3129.215 3169.056 3113.215 2644.733
\end{verbatim}

The story here is essentially the same as in part (b): The model with
random slope and random intercept is strongly favored by all four
criteria.

Not part of the required answer here, but, comparing parts (b) and (c)
it seems like the cross-level interaction is not needed.

\rm

\hypertarget{problem-4.}{%
\section{Problem 4.}\label{problem-4.}}

IMRAD and IDMRAD.

\hypertarget{problem-4a.}{%
\subsection{Problem 4(a).}\label{problem-4a.}}

\textit{(You do not need to turn anything in for this part.)} Compare
the papers

\begin{itemize}

\item \verb|COVID breakthrough rates in England - PIIS1473309921004606.pdf| and\\ \verb|COVID breakthrough... (appendix) 1-s2.0-S1473309921004606-mmc1.pdf|

\item \verb|menu pricing IDMRAD paper with technical appendix.pdf|
\end{itemize}

The COVID paper is a published IMRAD paper with a separate technical
appendix. The ``menu pricing'\,' paper is a simple IDMRAD paper (with
technical appendix attached) that you can use as a model for your final
project paper. Note similarities and differences. In particular, when
you write your paper, the technical appendix should show all the R code
and results that you used to answer the questions posed for the data
analysis (A small portion of the results in the appendix actually makes
it into the main paper). This is a little different from the published
COVID paper, which really just provides some supplementary details in
the technical appendix that didn't fit in the main paper.

\it

There really is nothing to turn in here. I just want you to see

\begin{itemize}

\item An example of a ``data'' section (present in the ``menu pricing'' paper but not in the COVID paper)

\item Difference in the content of the technical appendices of the two papers
\begin{itemize}
\item The technical appendix of the COVID paper doesn't contain much more computational detail than is in the main paper.  It is mainly other supplemental information for the paper.

\item The technical appendix for the ``menu pricing'' paper is much more all of the technical detail and computation that justifies or ``backs up'' the assertions in the main part of the paper.  With the more computational technical appendix I will be able to see how your computation corresponds to what you have written in the main body of the paper. 
\end{itemize}
Note also all of the cross-references in the main body of the paper to specific sections and pages of the technical appendix, in \underline{both} papers.  This is really needed so that I do not have to leaf through your appendix to find the right computation each time I want to check an assertion you made in the main part of your paper.   \bf Please include lots and lots of cross-references to specific pages of the technical appendix when you write your IDMRAD paper!\rm\it

\item The liberal use of subsections in the COVID paper to make material easy to find, for a reader who is in a hurry.  \bf Please use subsections to organize the material in your IDMRAD paper to make things easy to find for a reader in a hurry---whether it is me, or the TA, or your peer reviewers!\rm\it

\end{itemize}

\rm

\hypertarget{problem-3b.-1}{%
\subsection{Problem 3(b).}\label{problem-3b.-1}}

\textit{(You do need to turn something in for this part.)} Refer again
to the two papers in part (a). The difference between an IMRAD paper and
an IDMRAD paper is that the ``Data'\,' part of the \textbf{Methods}
section of an IMRAD paper is taked out of \textbf{Methods} and made into
its own \textbf{Data} section between \textbf{Introduction} and
\textbf{Methods}. \textbackslash{[}-0.75em{]}

Find the \textbf{Methods} section of the COVID paper, several pages in,
and within the \textbf{Methods} section find the subsections that
describe the data (rather than describing methods of analysis). Copy and
paste these sections (including subsection titles) into a separate
\textbf{Data} section, and turn that section in as an answer to this
question.

\it

The
\bf Data\rm\it{} section should consist of the following three subsections of the \bf Methods\rm\it{} section of this paper:

\begin{itemize}

\item Study design and participants
\item Risk factor variable definitions
\item Disease severity, duration, and symptom definitions

\end{itemize}

The full data section would look like this:

\bf Data.\rm\it

\underline{Study design and participants}

This prospective, community-based, nested, case-control study used data
from UK-based, adult (\(\geq\) -18 years) participants of the COVID
Symptom Study logged through a free mobile phone app developed by ZOE
(London, UK) and King's College London (London, UK).\(^{22}\) The app
was launched in the UK on March 24, 2020, and by July 4, 2021, had
nearly 4·5 million unique participants providing data by self-report or
proxy report. At registration, each participant reported baseline
demographic information (eg, age, sex, ethnicity, weight, height, and
health-care worker status), geographical location, and information on
health risk factors, including comorbidities, lifestyle, frailty, visits
to hospital, and adherence to mask-wearing guidance. Participants were
encouraged to self-report any of 32 pre specified symptoms (appendix p
1) daily, providing prospective longitudinal information on incident
symptoms. Those with new symptoms were prompted to book and take a
SARS-CoV-2 test. All users were encouraged to record any SARS-CoV-2
testing results (whether prompted by the app or otherwise), and, from
Dec 11, 2020, any SARS-CoV-2 vaccination and subsequent
symptoms.\(^{22}\) Users with missing or inconsistent information were
excluded from our analysis. The inclusion process for cases and controls
is shown in the appendix (p 14).

Cases had received a first or second dose of a COVID-19 vaccine since
Dec 8, 2020; had either a positive RT-PCR test or lateral flow antigen
test (LFAT) at least 14 days after their first vaccination (but before
their second; cases 1) or a positive RT-PCR test or LFAT at least 7 days
after their second vaccination (cases 2); and had no positive SARS-CoV-2
test before vaccination. If more than one positive test result was
reported, only the first positive test was selected. To identify risk
factors for postvaccination infection, we selected two control groups
among the vaccinated (since Dec 8, 2020) UK-based adult users of the
COVID Symptom Study app who had not tested positive for SARS-CoV-2
before vaccination: a control group of users reporting a negative RT-PCR
test or LFAT at least 14 days after their first vaccination but before
their second (controls 1) and a control group of users reporting a
negative RT-PCR test or LFAT at least 7 days after their second
vaccination (controls 2). Controls 1 and controls 2 were matched (1:1)
with cases 1 and cases 2, respectively, by use of the date of the
postvaccination COVID-19 test, health-care worker status, and sex. If
multiple negative tests were reported, the last test date was used for
matching.

To compare the disease profile of SARS-CoV-2 infection before and after
vaccination, we sub-selected participants from cases 1 and cases 2 who
had used the app for at least 14 consecutive days after testing positive
for SARS-CoV-2 (denoted as cases 3 and cases 4, respectively), so that
symptoms of infection could be assessed. Controls for the disease
profile analysis were those who reported a SARS-CoV-2-positive RT-PCR
test or LFAT, were unvaccinated until data censoring, and had used the
app for at least 14 consecutive days after the test. Among these users,
two groups were formed: controls 3 and controls 4, matching (1:1) with
cases 3 and cases 4, respectively, by the date of the positive COVID-19
test, health-care worker status, sex, body-mass index (BMI), and age.
Individuals in all case and control groups who did not report an RT-PCR
or LFAT test after Dec 8, 2020, were excluded.

For all control groups, we used a matching algorithm based on minimum
Euclidean distance\(^{23}\) between the vectors of the covariates, with
age, BMI, and the date of the test as numerical variables, and sex as a
binary variable multiplied by 100 to ensure balance between covariate
strengths. We considered health-care worker status, coded as a a
categorical variable in the app, as a numerical variable (0=not a
health-care worker; 1=health-care worker who does not interact with
patients; 2=health-care worker who does not treat patients;
3=health-care worker who interacts with patients; 4=health-care worker
who treats patients). Participants could only choose one of these
options.

All app users provided digital informed consent for data usage for
COVID-19-related research. In the UK, the app and the study were
approved by King's College London's ethics committee (REMAS number
18210; reference LRS-19/20--18210).

\underline{Risk factor variable definitions}

For this analysis, the outcome variable was case status (a self-reported
positive RT-PCR test or LFAT for SARS-CoV-2). We considered a
priori-defined risk factors for SARS-CoV-2 infection based on previous
evidence in unvaccinated individuals:16--18 age; BMI; self-reported
comorbidities (ie, cancer, diabetes, asthma, lung disease, heart
disease, and kidney disease), analysed individually as binary variables;
dependency level (frailty) assessed by the PRISMA-7 questionnaire, which
is embedded in app registration,\(^{24,25}\) as a binary variable
(PRISMA-7 score \(\geq\) -3 defined as frail and \textless3 defined as
not frail);\(^{26}\) local area Index of Multiple Deprivation (IMD; a
score ranging from 1 {[}most deprived{]} to 10 {[}least deprived{]}
estimating relative locality deprivation derived from postal code and
lower layer super output area) divided into low IMD (1--3), intermediate
IMD (4--7), and high IMD (8--10) groups;\(^{27}\) and four healthy
lifestyle factors (no current smoking, no obesity {[}BMI \textless30
kg/m²{]}, physical activity at least once per week {[}non-sedentary{]},
and a healthy diet pattern; appendix p 17). We also calculated a healthy
lifestyle score on the basis of these four lifestyle factors,\(^{28}\)
by which participants received 1 point for each healthy lifestyle factor
and the sum of the scores gave a total healthy lifestyle score ranging
from 0 to 4, with higher scores indicating a healthier lifestyle
(appendix p 17).

\underline{Disease severity, duration, and symptom definitions}

To compare the disease profile in vaccinated versus unvaccinated
individuals testing positive for SARS-CoV-2, we assessed disease
severity (asymptomatic or sympto- matic; more than five symptoms or five
or fewer symptoms reported in the first week of illness;20 and
self-reported presentation to hospital or no hospital presentation),
illness duration (duration of \textless28 days or \(\geq\) -28 days),
and individual symptom reports. Vaccination status was the exposure. For
cases, controls 3, and controls 4, symptoms were considered within a
window between 3 days before the COVID-19 test date and up to 14 days
after the test date (appendix p 1). This window was used because it
might have taken up to 3 days to request an RT-PCR test and receive a
result following symptom onset, and symptoms can occur up to 14 days
following SARS-CoV-2 exposure.\(^{29}\) \rm

\end{document}
