% 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}
\IfFileExists{xurl.sty}{\usepackage{xurl}}{} % add URL line breaks if available
\IfFileExists{bookmark.sty}{\usepackage{bookmark}}{\usepackage{hyperref}}
\hypersetup{
  pdftitle={Homework 02 Solutions},
  hidelinks,
  pdfcreator={LaTeX via pandoc}}
\urlstyle{same} % disable monospaced font for URLs
\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

\title{Homework 02 Solutions}
\author{}
\date{\vspace{-2.5em}2022-02-05}

\begin{document}
\maketitle

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(arm)  }\DocumentationTok{\#\# includes lme4,MASS, Matrix}
\FunctionTok{library}\NormalTok{(ggplot2)}
\end{Highlighting}
\end{Shaded}

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

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

In lecture 06 (crossed random effects, cross-level interactions), six
different models were fitted to the London Schools data, but only the
\texttt{lmer()} model specification was given for each model. For each
of these six models, write out the two-level MLM formulation of the
model. You can use mathematical variable names (\(y_i\), \(x_{3i}\),
etc.) or English variable names (LRT\(_i\), etc.), as long as it is
clear which variables in the data frame are in each of your six MLM
models. Be sure to state all distribution assumptions, and be sure that
the indexing (\(i\), \(j\), \(j[i]\) and so forth) is clear.

\it
In all that follows we will use the definitions::
\begin{eqnarray*}
    i & = & 1\ldots n \mbox{~~(students)}\\
    j & = & 1\ldots J \mbox{~~(schools)}\\
    j[i] & = & \mbox{the school $j$ that student $i$ is in}
\end{eqnarray*}
See my response on Piazza (\textcolor{blue}{\url{https://piazza.com/class/kykc45ull1j7dh?cid=19}}) for general guidance on how to do these\ldots
\rm

\hypertarget{part-1a}{%
\subsection{Part 1(a)}\label{part-1a}}

Random Intercept Only: \verb`Y ~ 1 + (1 | school)` \begin{eqnarray*}
    Y_i & = & \alpha_{0[i]} + \epsilon_i 
    ~~, ~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
    \alpha_{0j} & = & \beta_0 + \eta_{0j}
    ~~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2)
\end{eqnarray*}

\hypertarget{part-1b}{%
\subsection{Part 1(b)}\label{part-1b}}

Random Intercept, Fixed Slope: \verb`Y ~ LRT + (1 | school)`
\begin{eqnarray*}
    Y_i & = & \alpha_{0[i]} + \beta_1 \cdot \mbox{LRT}_i + \epsilon_i 
    ~~, ~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
    \alpha_{0j} & = & \beta_0 + \eta_{0j}
    ~~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2)
\end{eqnarray*}

\hypertarget{part-1c}{%
\subsection{Part 1(c)}\label{part-1c}}

Random Intercept, Random Slope: \verb`Y ~ LRT + (LRT | school)`
\begin{eqnarray*}
    Y_i & = & \alpha_{0[i]} + \alpha_{1j[i]} \cdot \mbox{LRT}_i + \epsilon_i 
    ~~, ~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
    \alpha_{0j} & = & \beta_0 + \eta_{0j} 
    ~~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2) \\
    \alpha_{1j} & = & \beta_1 + \eta_{1j} 
    ~~,~~ \eta_{1j} \stackrel{iid}{\sim} N(0,\tau_1^2) \\
    && \mathop{\mathrm{Corr}}(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*}

\hypertarget{part-1d}{%
\subsection{Part 1(d)}\label{part-1d}}

Random Intercept, Random Slope, CWC:
\verb`Y ~ LRT.cwc + (1 + LRT.cwc | school)`

\it
Let $\overline{LRT}_j$ be the sample average of LRT in school $j$ and let LRT.cwc$_i$ $=$ LRT$_i$ $-$ LRT.cwc$_{j[i]}$.  Then the MLM is
\rm
\begin{eqnarray*}
    Y_i & = & \alpha_{0[i]} + \alpha_{1j[i]} \cdot \mbox{LRT.cwc}_i + \epsilon_i 
    ~~, ~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
    \alpha_{0j} & = & \beta_0 + \eta_{0j} 
    ~~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2) \\
    \alpha_{1j} & = & \beta_1 + \eta_{1j} 
    ~~,~~ \eta_{1j} \stackrel{iid}{\sim} N(0,\tau_1^2) \\
    && \mathop{\mathrm{Corr}}(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*}

\hypertarget{part-1e}{%
\subsection{Part 1(e)}\label{part-1e}}

Random Intercept, individual and group predictor:
\verb`Y ~ LRT + sch.avg + (1 + LRT | school)`

\it
Using the same definitions as in the previous part, let \mbox{\rm sch.avj}$_j$ $=$ $\overline{LRT}_j$. Then the MLM is
\rm
\begin{eqnarray*}
    Y_i & = & \alpha_{0[i]} + \alpha_{1j[i]} \cdot \mbox{LRT}_i + \beta_2 \cdot \mbox{sch.avg}_{j[i]} + \epsilon_i 
    ~~, ~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
    \alpha_{0j} & = & \beta_0 + \eta_{0j} 
    ~~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2) \\
    \alpha_{1j} & = & \beta_1 + \eta_{1j} 
    ~~,~~ \eta_{1j} \stackrel{iid}{\sim} N(0,\tau_1^2) \\
    && \mathop{\mathrm{Corr}}(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*}

\hypertarget{part-1f}{%
\subsection{Part 1(f)}\label{part-1f}}

Random Intercept, Cross-Level Interaction:
\verb`Y ~ LRT * sch.avg + (1 + LRT | school)`

\it
With all the same definitions as above, the MLM is
\rm
\begin{eqnarray*}
    Y_i & = & \alpha_{0[i]} + \alpha_{1j[i]} \cdot \mbox{LRT}_i + \beta_2 \cdot \mbox{sch.avg}_{j[i]} 
    + \beta_3 \cdot \mbox{LRT}_i \cdot \mbox{sch.avg}_{j[i]} + \epsilon_i 
    ~~, ~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
    \alpha_{0j} & = & \beta_0 + \eta_{0j} 
    ~~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2) \\
    \alpha_{1j} & = & \beta_1 + \eta_{1j} 
    ~~,~~ \eta_{1j} \stackrel{iid}{\sim} N(0,\tau_1^2) \\
    && \mathop{\mathrm{Corr}}(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*}

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

Return to the model on slide 10 of lecture 05 (random slopes,
correlation and no centering). Use the AIC.ml, BIC.ml, and DIC.lm
functions defined at the end of lecture 06 (crossed random effects,
cross-level interactions) as needed for the following quesitons.

\hypertarget{part-2a}{%
\subsection{Part 2(a)}\label{part-2a}}

Write the multilevel model that has a cross-level interaction between
log.uranium and floor. Interpret this model: what does it say about how
log.uranium and floor influence log.radon in each house, and how is it
different from the model without the cross-level interaction?

\it

Here is the MLM with the cross-level interaction: \begin{eqnarray*}
    y_i & = & \alpha_{0[i]} + \alpha_{1j[i]} \cdot x_i + \beta_2 \cdot u_{j[i]} 
    + \beta_3 \cdot x_i \cdot u_{j[i]} + \epsilon_i 
    ~~, ~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
    \alpha_{0j} & = & \beta_0 + \eta_{0j} 
    ~~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2) \\
    \alpha_{1j} & = & \beta_1 + \eta_{1j} 
    ~~,~~ \eta_{1j} \stackrel{iid}{\sim} N(0,\tau_1^2) \\
    && \mathop{\mathrm{Corr}}(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*} where \(y_i = \log(radon)\) reading for the \(i^{th}\)
house, \(x_i\) \(=\) 0 if the radon level was measured in the basement;
1 if it was measured on the first floor, and \(j[i]\) is the county
\(j\) that house \(i\) is in.

And here is the MLM without the cross-level interaction:
\begin{eqnarray*}
    y_i & = & \alpha_{0[i]} + \alpha_{1j[i]} \cdot x_i + \beta_2 \cdot u_{j[i]} + \epsilon_i 
    ~~, ~~ \epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
    \alpha_{0j} & = & \beta_0 + \eta_{0j} 
    ~~,~~ \eta_{0j} \stackrel{iid}{\sim} N(0,\tau_0^2) \\
    \alpha_{1j} & = & \beta_1 + \eta_{1j} 
    ~~,~~ \eta_{1j} \stackrel{iid}{\sim} N(0,\tau_1^2) \\
    && \mathop{\mathrm{Corr}}(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*} In the model with the cross-level interaction,
\(u_{j[i]}\) affects the slope on \(x_i\), which we can see by just
collecting terms that involve a factor of \(x_i\): \[
  \mbox{slope on $x_i$} = \beta_1 + \beta_2\cdot u_{j[i]} + \eta_{1j[i]}
\] whereas in the model without the cross-level interaction, the slope
on \(x_i\) doesn't depend on \(u_{j[i]}\) at all: \[
  \mbox{slope on $x_i$} = \beta_1 + \eta_{1j[i]}
\] \rm

\hypertarget{part-2b}{%
\subsection{Part 2(b)}\label{part-2b}}

Which model fits the data better, with or without the cross-level
interaction between floor and log.uranium?

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

\NormalTok{mn.radon }\OtherTok{\textless{}{-}} \FunctionTok{read.table}\NormalTok{(}\StringTok{"mn{-}radon.txt"}\NormalTok{,}\AttributeTok{header=}\NormalTok{T)}

\NormalTok{lmer}\FloatTok{.0} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(log.radon }\SpecialCharTok{\textasciitilde{}}\NormalTok{ floor               }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1} \SpecialCharTok{+}\NormalTok{ floor }\SpecialCharTok{|}\NormalTok{ county), }\AttributeTok{data=}\NormalTok{mn.radon)}
\NormalTok{lmer}\FloatTok{.1} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(log.radon }\SpecialCharTok{\textasciitilde{}}\NormalTok{ floor }\SpecialCharTok{+}\NormalTok{ log.uranium }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1} \SpecialCharTok{+}\NormalTok{ floor }\SpecialCharTok{|}\NormalTok{ county), }\AttributeTok{data=}\NormalTok{mn.radon)}
\NormalTok{lmer}\FloatTok{.2} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(log.radon }\SpecialCharTok{\textasciitilde{}}\NormalTok{ floor }\SpecialCharTok{*}\NormalTok{ log.uranium }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1} \SpecialCharTok{+}\NormalTok{ floor }\SpecialCharTok{|}\NormalTok{ county), }\AttributeTok{data=}\NormalTok{mn.radon)}

\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{), AIC.ml),}
             \AttributeTok{BIC=}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{list}\NormalTok{(lmer}\FloatTok{.0}\NormalTok{, lmer}\FloatTok{.1}\NormalTok{,lmer}\FloatTok{.2}\NormalTok{), BIC.ml),}
             \AttributeTok{DIC=}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{list}\NormalTok{(lmer}\FloatTok{.0}\NormalTok{, lmer}\FloatTok{.1}\NormalTok{,lmer}\FloatTok{.2}\NormalTok{), DIC.ml))}
 
\FunctionTok{colnames}\NormalTok{(res) }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\StringTok{"(w/o log.uranium)"}\NormalTok{, }\StringTok{"(w/o interaction)"}\NormalTok{, }\StringTok{"(with interaction)"}\NormalTok{)}

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

\begin{verbatim}
##                        AIC     BIC     DIC
## (w/o log.uranium)  2173.11 2202.05 2161.11
## (w/o interaction)  2131.60 2165.37 2117.60
## (with interaction) 2135.22 2173.81 2119.22
\end{verbatim}

\it
As you can see, all three information criteria like \verb|log.uranium| in the model, and they all prefer the model with a main effect, better than the model with the interaction. 
\rm

\hypertarget{part-2c}{%
\subsection{Part 2(c)}\label{part-2c}}

Construct a new Level 2 variable nj that is equal to the number of
houses sampled in each county. Like log.uranium, this variable will have
the same value for each house in the county. Fit two models like the one
on slide 10 of lecture 05, but with nj replacing log.uranium: one model
without a cross-level interaction between nj and floor, and one model
with the cross-level interaction. Of the four models (with nj or with
log.uranium, and with or without the cross-level interaction) which
one(s) fit the data best? Accompany your results with a short paragraph
interpreting them.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{nj }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(mn.radon,}\FunctionTok{unlist}\NormalTok{(}\FunctionTok{lapply}\NormalTok{(}\FunctionTok{split}\NormalTok{(county,county),}
                                  \ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{rep}\NormalTok{(}\FunctionTok{length}\NormalTok{(x),}\FunctionTok{length}\NormalTok{(x))))) }

\NormalTok{lmer}\FloatTok{.3} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(log.radon }\SpecialCharTok{\textasciitilde{}}\NormalTok{ floor }\SpecialCharTok{+}\NormalTok{ nj }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1} \SpecialCharTok{+}\NormalTok{ floor }\SpecialCharTok{|}\NormalTok{ county), }\AttributeTok{data=}\NormalTok{mn.radon)}
\NormalTok{lmer}\FloatTok{.4} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(log.radon }\SpecialCharTok{\textasciitilde{}}\NormalTok{ floor }\SpecialCharTok{*}\NormalTok{ nj }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1} \SpecialCharTok{+}\NormalTok{ floor }\SpecialCharTok{|}\NormalTok{ county), }\AttributeTok{data=}\NormalTok{mn.radon)}

\NormalTok{res }\OtherTok{\textless{}{-}} \FunctionTok{rbind}\NormalTok{(}\AttributeTok{AIC=}\FunctionTok{sapply}\NormalTok{(}\FunctionTok{list}\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{.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{.1}\NormalTok{,lmer}\FloatTok{.2}\NormalTok{,lmer}\FloatTok{.3}\NormalTok{,lmer}\FloatTok{.4}\NormalTok{), DIC.ml))}
 
\FunctionTok{colnames}\NormalTok{(res) }\OtherTok{\textless{}{-}} \FunctionTok{c}\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.1 2131.60 2165.37 2117.60
## lmer.2 2135.22 2173.81 2119.22
## lmer.3 2169.47 2203.23 2155.47
## lmer.4 2171.26 2209.85 2155.26
\end{verbatim}

\it
The two models with \verb|nj| do quite a bit worse, on all three fit measures, than the models with \verb|log.uranium|.  And, as before, the model with only a main effect for \verb|log.uranium| an no interaction is preferred by all of AIC, BIC, and DIC.
\rm

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

Although it is usual to keep the correlation between random effects in a
multilevel model, there are times when we'd like to formally establish
whether the data support having the correlation in the model or not. One
way to do that is with a likelihood ratio test. Recall that the
likelihood ratio test for testing model \(M_0\) vs model \(M_1\), where
\(M_0\) is obtained from \(M_1\) by linear constraint(s), computes the
log-likelihood ratio test statistic \[
       -2[\log(likelihood_{M_0}) - \log(likelihood_{M_1})]
\] where the likelihoods are evaluated at the maximum likelihood values
of the parameters (e.g.~as fitted by \texttt{lmer()}), and compares the
value with a chi-squared distribution with degrees of freedom equal to
the number of linear constraints. A statistically significant result
rejects \(M_0\) in favor of \(M_1\) and conversely a nonsignificant
result suggests that there is not enough evidence in the data to reject
\(M_0\).

For the model on slide 10 of lecture 05 (random slopes, correlation and
no centering), use a likelihood ratio test to determine whether there is
sufficient evidence in the data to keep the correlation \(\rho\) between
\(\eta_{0j}\) and \(\eta_{1j}\) in the model.
\textit{(Note: For reasons that we will talk about next week, you need to set \texttt{REML=FALSE} in \texttt{lmer()} when you fit $M_0$ and $M_1$ for the likelihood ratio test.  You will also find the \texttt{logLik()} function in R useful for this exercise.)}
Be sure to report the number of degrees of freedom of the test, the
value of the test statistic, and the criterion you are using to
determine statistical significance.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{lmer}\FloatTok{.1} \OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(log.radon }\SpecialCharTok{\textasciitilde{}}\NormalTok{ floor }\SpecialCharTok{+}\NormalTok{ log.uranium }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1} \SpecialCharTok{+}\NormalTok{ floor }\SpecialCharTok{|}\NormalTok{ county), }\AttributeTok{data=}\NormalTok{mn.radon)}
\NormalTok{lmer.}\FloatTok{1.}\NormalTok{c0 }\OtherTok{\textless{}{-}} \FunctionTok{lmer}\NormalTok{(log.radon }\SpecialCharTok{\textasciitilde{}}\NormalTok{ floor }\SpecialCharTok{+}\NormalTok{ log.uranium }\SpecialCharTok{+}\NormalTok{ (}\DecValTok{1} \SpecialCharTok{+}\NormalTok{ floor }\SpecialCharTok{||}\NormalTok{ county), }\AttributeTok{data=}\NormalTok{mn.radon)}

\DocumentationTok{\#\# Since there is only one linear constraint needed (rho = 0) to get from lmer.1 }
\DocumentationTok{\#\# to lmer.1.c0, we will be referring the likelihood ratio test statistic to a }
\DocumentationTok{\#\# chi{-}squared distribution on with df=1.}

\NormalTok{LL }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(M) \{ }\FunctionTok{logLik}\NormalTok{(}\FunctionTok{update}\NormalTok{(M,}\AttributeTok{REML=}\NormalTok{F)) \}}

\NormalTok{LR.stat }\OtherTok{\textless{}{-}} \SpecialCharTok{{-}}\DecValTok{2}\SpecialCharTok{*}\NormalTok{(}\FunctionTok{LL}\NormalTok{(lmer.}\FloatTok{1.}\NormalTok{c0) }\SpecialCharTok{{-}} \FunctionTok{LL}\NormalTok{(lmer}\FloatTok{.1}\NormalTok{))}

\NormalTok{pval }\OtherTok{\textless{}{-}} \FunctionTok{pchisq}\NormalTok{(LR.stat,}\AttributeTok{df=}\DecValTok{1}\NormalTok{,}\AttributeTok{lower.tail=}\NormalTok{F)}

\FunctionTok{round}\NormalTok{(}\FunctionTok{rbind}\NormalTok{(}\AttributeTok{df=}\DecValTok{1}\NormalTok{,}\AttributeTok{LR.stat=}\NormalTok{LR.stat,}\AttributeTok{pval=}\NormalTok{pval),}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##         [,1]
## df      1.00
## LR.stat 0.20
## pval    0.65
\end{verbatim}

\it 
Clearly, for any reasonable cutoff like 0.01 or 0.05, the p-value here is huge, so there is not enough evidence in the data to reject the null hypothesis that the model without correlation \verb|lmer.1.c0| fits adequately, vs the alternative \verb|lmer.1|.
\rm

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

The file allvars.csv contains CD4 percentages (\verb|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
(\textcolor{blue}{\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 \verb|CD4PCT|`s at each visit, and the Level 2
``group'' is a child. How the \verb|CD4PCT|'s change over time for each
child is the primary question of interest. Because of skewing, you
should replace \verb|CD4PCT| with sqrt(\verb|CD4PCT|) for all of the
following questions.

\hypertarget{part-4a}{%
\subsection{Part 4(a)}\label{part-4a}}

For the first 12 children only (\texttt{newpid $\leq$ 12}), make a facet
plot of \verb|CD4PCT| vs \verb|VISIT| number. Use different plotting
point colors for children who have \verb|treatmnt|=1 and
\verb|treatmnt|=0\footnote{Presumably, ``treatmnt'' stands for some sort of treatment to improve childrens' HIV status.}.
Do the same for \verb|CD4PCT|
vs.~\verb|visage|\footnote{If you are using \texttt{ggplot()}, you may wish to set \texttt{scales=''free\_x''} in the \texttt{facet\_wrap()} function.}
(the age of the child on each visit). 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{hw02-solutions_files/figure-latex/unnamed-chunk-6-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{hw02-solutions_files/figure-latex/unnamed-chunk-7-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{part-4b}{%
\subsection{Part 4(b)}\label{part-4b}}

Build a multilevel model with random slopes and intercepts for all the
children, using a measure of time (\verb|VISIT| or \verb|visage|,
whichever you answered in part (a)) as a Level 1 predictor, and
\verb|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{hw02-solutions_files/figure-latex/unnamed-chunk-9-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{hw02-solutions_files/figure-latex/unnamed-chunk-11-1.pdf}

\hypertarget{part-4c}{%
\subsection{Part 4(c)}\label{part-4c}}

Try to expand the model by adding a (cross-level) interaction between
\verb|treatmnt| and your measure of time. 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 \verb|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{hw02-solutions_files/figure-latex/unnamed-chunk-13-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{hw02-solutions_files/figure-latex/unnamed-chunk-15-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{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{hw02-solutions_files/figure-latex/unnamed-chunk-17-1.pdf}

\end{document}
