---
title: "Homework 03 Solutions"
date: "2022-02-16"
output: 
  pdf_document:
    keep_tex: true
    fig_caption: true
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \usepackage{hyperref}
  - \usepackage{xcolor}
  - \DeclareMathOperator{\Cov}{Cov}
  - \DeclareMathOperator{\Corr}{Corr}
  - \DeclareMathOperator{\Var}{Var}
  - \DeclareMathOperator{\upmodels}{\perp \!\!\! \perp}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, include = TRUE, message = FALSE)
# library(knitr)
```


```{r}
library(arm)  ## includes lme4, MASS, Matrix
library(ggplot2); theme_set(theme_bw())
library(gridExtra) ## to arrange ggplots...
library(HLMdiag)
library(cAIC4)
```

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


# Problem 1.

Return to the CD4 data from HW02, and consider the model
\begin{verbatim}
lmer.1 <- lmer(sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (1 + VISIT | newpid),
               data=cd4)
\end{verbatim}
Bring together all the residuals (except random effects residuals) that we have talked about, using library(HLMdiag):
```{r}
cd4 <- read.csv("allvar.csv",header=T)

## It's worth doing a little exploration of the data here,
## especially to look at missing values.

apply(cd4,2,function(x) mean(is.na(x)))
```


```{r}
## we only need the variables newpid, VISIT, visage, CD4PCT and treatmnt
## so let's get rid of some (but not all) missing data problems
## by deleting the other variables

cd4 <- with(cd4,data.frame(newpid=newpid,
                           VISIT=VISIT, 
                           visage=visage, 
                           CD4PCT=CD4PCT, 
                           treatmnt=treatmnt))

## and since we still have some missing data, we will just 
## delete the rows that continue to have NA's (not a great
## practice in general, but good enough for this exercise...)
cd4 <- cd4[!apply(cd4,1,function(x) any(is.na(x))),]

## set up the sqrt of the response variable...
cd4$sqrt.CD4PCT <- sqrt(cd4$CD4PCT)

## and for coloring ggplot elements it will better if treatmnt is a factor...
cd4$treatmnt <- as.factor(cd4$treatmnt)

## select the first 12 kids...
first.12 <- (cd4$newpid <= 12)

## Fit the base model for this exercise...
display(lmer.1 <- lmer(sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (1 + VISIT | newpid),
               data=cd4))
```


```{r}
r.1 <- hlm_resid(lmer.1,level=1,include.ls=F)
r.1s <- hlm_resid(lmer.1,level=1,include.ls=F,standardize=T)
r.2 <- hlm_resid(lmer.1,level="newpid",include.ls=F)
r.2s <- hlm_resid(lmer.1,level="newpid",include.ls=F,standardize=T)
names(r.1)
```


```{r}
names(r.1s)
```


```{r}
names(r.2)
```


```{r}
names(r.2s)
```


## Problem 1(a).

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).
```{r}

ggplot(r.1[first.12,],aes(x=.mar.fitted,y=.mar.resid)) +
  facet_wrap( ~ newpid,scales="free_x") +
  geom_point(aes(color=treatmnt)) +
  geom_abline(intercept=0,slope=0)


```

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

## Problem 1(b)

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).
```{r}
ggplot(r.1,aes(x=.mar.fitted,y=.mar.resid)) +
  geom_point(aes(color=treatmnt)) +
  geom_smooth()
```

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}

## Problem 1(c)

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). 

```{r}
ggplot(r.1,aes(x=.fitted,y=.resid)) +
  geom_point(aes(color=treatmnt)) +
  geom_smooth()
```

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}

```{r, fig.height=3}
y.grouped <- with(cd4,sapply(split(sqrt.CD4PCT,newpid),mean))
tx.grouped <- with(cd4,sapply(split(treatmnt,newpid),function(x) x[1]))

g1 <- ggplot(r.2,aes(x=y.grouped,y=.ranef.intercept)) +
  geom_point(aes(color=tx.grouped)) +
  xlab("mean(sqrt.CD4PCT) within child") + 
  ylab("eta0") +
  theme(legend.position="none")

g2 <- ggplot(r.2,aes(x=y.grouped,y=.ranef.visit)) +
  geom_point(aes(color=tx.grouped)) +
  xlab("mean(sqrt.CD4PCT) within child") + 
  ylab("eta1") +
  theme(legend.position="none")

grid.arrange(g1,g2,ncol=2)

```

\clearpage

## Problem 1(d)

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.

```{r}
g1 <- ggplot(r.1s,aes(sample=.std.resid)) + 
  geom_qq() +
  geom_qq_line() +
  xlim(-3.5,3.5) + 
  ylim(-8,8) +
  xlab("Normal Quantiles") + 
  ylab("Data Quantiles") +
  ggtitle("Standardized conditional residuals")

g2 <- ggplot(r.1s,aes(sample=.chol.mar.resid)) + 
  geom_qq() +
  geom_qq_line() +
  xlim(-3.5,3.5) + 
  ylim(-8,8) +
  xlab("Normal Quantiles") + 
  ylab("Data Quantiles") +
  ggtitle("Cholesky marginal residuals")

g3 <- ggplot(r.2s,aes(sample=.std.ranef.intercept)) + 
  geom_qq() +
  geom_qq_line() +
  xlim(-3.5,3.5) + 
  ylim(-8,8) +
  xlab("Normal Quantiles") + 
  ylab("Data Quantiles") +
  ggtitle("Standardized eta0")

g4 <- ggplot(r.2s,aes(sample=.std.ranef.visit)) + 
  geom_qq() +
  geom_qq_line() +
  xlim(-3.5,3.5) + 
  ylim(-8,8) +
  xlab("Normal Quantiles") + 
  ylab("Data Quantiles") +
  ggtitle("Standardized eta1")

grid.arrange(g1,g2,g3,g4,ncol=2)
```
\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


## Problem 2.

Continuing with the CD4 data\ldots

## Problem 2(a)

Make a table giving values of AIC, BIC, DIC, and cAIC (you compared two of these on the last assignment,using just AIC, BIC and DIC):
```{r}
lmer.1 <- lmer(sqrt.CD4PCT ~ 1 + visage +            (1+visage|newpid),data=cd4,REML=F)
lmer.2 <- lmer(sqrt.CD4PCT ~ 1 + visage + treatmnt + (1+visage|newpid),data=cd4,REML=F)
lmer.3 <- lmer(sqrt.CD4PCT ~ 1 + visage * treatmnt + (1+visage|newpid),data=cd4,REML=F)
lmer.4 <- lmer(sqrt.CD4PCT ~ 1 + VISIT  +            (1+VISIT|newpid),data=cd4,REML=F)
lmer.5 <- lmer(sqrt.CD4PCT ~ 1 + VISIT  + treatmnt + (1+VISIT|newpid),data=cd4,REML=F)
lmer.6 <- lmer(sqrt.CD4PCT ~ 1 + VISIT  * treatmnt + (1+VISIT|newpid),data=cd4,REML=F)

DIC <- function(M) {
         if(class(M)=="lm") { x <- AIC(M) }
         else { x <- unname(extractDIC(M)) }
         return(x)
       }
  
CAIC <- function(M) {cAIC(M)$caic}

IC <- function(M) {
    x <-c(AIC=AIC(M),BIC=BIC(M),DIC=DIC(M),cAIC=CAIC(M))
    return(x)
}

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

IC.table(lmer.1,lmer.2,lmer.3,lmer.4,lmer.5,lmer.6)
```

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.5 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

## Problem 2(b)

Make a table giving values of  AIC, BIC, DIC, and cAIC for the following models:
```{r}
lm.7 <- lm(sqrt.CD4PCT ~ 1 + VISIT + treatmnt,data=cd4)
lmer.8 <- lmer(sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1|newpid),data=cd4,REML=F)
lmer.9 <- lmer(sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (0+VISIT|newpid),data=cd4,REML=F)
lmer.10 <- lmer(sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1+VISIT|newpid),data=cd4,REML=F)

IC.table(lm.7,lmer.8,lmer.9,lmer.10)
```

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 BIC 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

## Problem 2(c)

Repeat part (b) but with the interaction \verb|VISIT * treatmnt| in each model instead of just the main effects \verb|VISIT + treatmnt|.
```{r}
lm.11 <- lm(sqrt.CD4PCT ~ 1 + VISIT * treatmnt,data=cd4)
lmer.12 <- lmer(sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (1|newpid),data=cd4,REML=F)
lmer.13 <- lmer(sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (0+VISIT|newpid),data=cd4,REML=F)
lmer.14 <- lmer(sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (1+VISIT|newpid),data=cd4,REML=F)

IC.table(lm.11,lmer.12,lmer.13,lmer.14)
```

\it

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
