---
title: "Homework 02 Solutions"
date: "2022-02-05"
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)
```

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

# 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 (\blue{\url{https://piazza.com/class/kykc45ull1j7dh?cid=19}}) for general guidance on how to do these\ldots
\rm

## Part 1(a)

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

## Part 1(b)

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

## Part 1(c)

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) \\
    && \Corr(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*}

## Part 1(d)

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

\it
Let $\bar{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) \\
    && \Corr(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*}

## Part 1(e)

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$ $=$ $\bar{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) \\
    && \Corr(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*}

## Part 1(f)

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) \\
    && \Corr(\eta_{0j},\eta_{1j}) ~ = ~ \rho
\end{eqnarray*}

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

## Part 2(a)

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

## Part 2(b)

Which model fits the data better, with or without the cross-level interaction between floor and log.uranium?
 
```{r}
AIC.ml <- function(M) {AIC(update(M,REML=F))}        
BIC.ml <- function(M) {BIC(update(M,REML=F))}        
DIC.ml <- function(M) {extractDIC(update(M,REML=F))} 

mn.radon <- read.table("mn-radon.txt",header=T)

lmer.0 <- lmer(log.radon ~ floor               + (1 + floor | county), data=mn.radon)
lmer.1 <- lmer(log.radon ~ floor + log.uranium + (1 + floor | county), data=mn.radon)
lmer.2 <- lmer(log.radon ~ floor * log.uranium + (1 + floor | county), data=mn.radon)

res <- rbind(AIC=sapply(list(lmer.0, lmer.1,lmer.2), AIC.ml),
             BIC=sapply(list(lmer.0, lmer.1,lmer.2), BIC.ml),
             DIC=sapply(list(lmer.0, lmer.1,lmer.2), DIC.ml))
 
colnames(res) <- c("(w/o log.uranium)", "(w/o interaction)", "(with interaction)")

t(round(res,2))

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

## Part 2(c)

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. 

 
```{r}
nj <- with(mn.radon,unlist(lapply(split(county,county),
                                  function(x) rep(length(x),length(x))))) 

lmer.3 <- lmer(log.radon ~ floor + nj + (1 + floor | county), data=mn.radon)
lmer.4 <- lmer(log.radon ~ floor * nj + (1 + floor | county), data=mn.radon)

res <- rbind(AIC=sapply(list(lmer.1,lmer.2,lmer.3,lmer.4), AIC.ml),
             BIC=sapply(list(lmer.1,lmer.2,lmer.3,lmer.4), BIC.ml),
             DIC=sapply(list(lmer.1,lmer.2,lmer.3,lmer.4), DIC.ml))
 
colnames(res) <- c("lmer.1", "lmer.2", "lmer.3","lmer.4")

t(round(res,2))

```
 

\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




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

 
```{r}
lmer.1 <- lmer(log.radon ~ floor + log.uranium + (1 + floor | county), data=mn.radon)
lmer.1.c0 <- lmer(log.radon ~ floor + log.uranium + (1 + floor || county), data=mn.radon)

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

LL <- function(M) { logLik(update(M,REML=F)) }

LR.stat <- -2*(LL(lmer.1.c0) - LL(lmer.1))

pval <- pchisq(LR.stat,df=1,lower.tail=F)

round(rbind(df=1,LR.stat=LR.stat,pval=pval),2)
```
 

\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

# 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 (\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.

## Part 4(a)

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
```{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)))

## 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)
```

\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
```{r}
g <- ggplot(cd4[first.12,],aes(x=VISIT,y=sqrt.CD4PCT)) +
  facet_wrap( ~ newpid) + 
  geom_point(aes(color=treatmnt))

g
```
\clearpage

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

```{r}
h <- ggplot(cd4[first.12,],aes(x=visage,y=sqrt.CD4PCT)) +
  facet_wrap( ~ newpid, scales="free_x") + 
  geom_point(aes(color=treatmnt))

h
```
 

\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

## Part 4(b)

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
 
```{r}
display(lmer.1 <- lmer(sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1+VISIT|newpid), data=cd4))
```
 

\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

 
```{r}
## fixef(lmer.1)  ## gives the betas for intercept and slopes on VISIT & treatmnt

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

coef.1 <- coef(lmer.1)$newpid[,1:3] ## coef.1[,1] = beta0 + eta0  (random intcpt)
                                    ## coef.1[,2] = beta1 + eta1  (random slope)
                                    ## coef.1[,3] = beta2         (coef on treatmnt)

tx <- with(cd4,sapply(split(treatmnt,newpid), function(x) x[1]))  ## one per child

params <- data.frame(newpid=sort(unique(cd4$newpid)),
                     alpha0=coef.1[,1],
                     alpha1=coef.1[,2],
                     beta2=coef.1[,3],
                     tx 
                     )

g + geom_abline(data=params[params$newpid<=12,],aes(intercept=alpha0 + ifelse(tx==2,beta2,0),
                                                    slope=alpha1,
                                                    color=tx))
```
 
\it
\mbox{}\\
Here's the model summary using \verb|visage| as the measure of time:
\mbox{}\\
\rm
 
```{r}
display(lmer.2 <- lmer(sqrt.CD4PCT ~ 1 + visage + treatmnt + (1+visage|newpid),data=cd4))
```
 
\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
 
```{r}
coef.2 <- coef(lmer.2)$newpid[,1:3]  ## coef.2[,1] = beta0 + eta0  (random intcpt)
                                     ## coef.2[,2] = beta1 + eta1  (random slope)
                                     ## coef.2[,3] = beta2         (coef on treatmnt)

tx <- with(cd4,sapply(split(treatmnt,newpid), function(x) x[1]))  ## one per child

params <- data.frame(newpid=sort(unique(cd4$newpid)),
                     alpha0=coef.2[,1],
                     alpha1=coef.2[,2],
                     beta2=coef.2[,3],
                     tx 
                     )

h + geom_abline(data=params[params$newpid<=12,],aes(intercept=alpha0 + ifelse(tx==2,beta2,0),
                                                    slope=alpha1,
                                                    color=tx))

```

 

## Part 4(c)

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
 
```{r}
display(lmer.3 <- lmer(sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (1+VISIT|newpid), data=cd4))
```
 
\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
 
```{r}
coef.3 <- coef(lmer.3)$newpid[,1:4]  ## coef.3[,1] = beta0 + eta0  (random intcpt)
                                     ## coef.3[,2] = beta1 + eta1  (random slope on time)
                                     ## coef.3[,3] = beta2         (coef on treatmnt)
                                     ## coef.3[,4] = beta3         (coef on interaction)

tx <- with(cd4,sapply(split(treatmnt,newpid), function(x) x[1]))  ## one per child

params <- data.frame(newpid=sort(unique(cd4$newpid)),
                     alpha0=coef.3[,1],
                     alpha1=coef.3[,2],
                     beta2=coef.3[,3],
                     beta3=coef.3[,4],
                     tx 
                     )

g + geom_abline(data=params[params$newpid<=12,],aes(intercept=alpha0 + ifelse(tx==2,beta2,0),
                                                    slope=alpha1 + ifelse(tx==2,beta3,0),
                                                    color=tx))
```
 
\it
\mbox{}\\
Here's the cross level interaction model, using \verb|visage| as time
\mbox{}\\
\rm

 
```{r}
display(lmer.4 <- lmer(sqrt.CD4PCT ~ 1 + visage * treatmnt + (1+visage|newpid),data=cd4))
```
 
\it
Again, no significant effects for \verb|visage|, and a rather large $\hat\tau_0$.

Here's the facet plot...
\rm
 
```{r}
coef.4 <- coef(lmer.4)$newpid[,1:4]  ## coef.4[,1] = beta0 + eta0  (random intcpt)
                                     ## coef.4[,2] = beta1 + eta1  (random slope on time)
                                     ## coef.4[,3] = beta2         (coef on treatmnt)
                                     ## coef.4[,4] = beta3         (coef on interaction)

tx <- with(cd4,sapply(split(treatmnt,newpid), function(x) x[1]))  ## one per child

params <- data.frame(newpid=sort(unique(cd4$newpid)),
                     alpha0=coef.4[,1],
                     alpha1=coef.4[,2],
                     beta2=coef.4[,3],
                     beta3=coef.4[,4],
                     tx 
                     )

h + geom_abline(data=params[params$newpid<=12,],aes(intercept=alpha0 + ifelse(tx==2,beta2,0),
                                                    slope=alpha1 + ifelse(tx==2,beta3,0),
                                                    color=tx))
```

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

lmer.0 <- lmer(sqrt.CD4PCT ~ 1 + VISIT + (1+VISIT|newpid),data=cd4)

res <- rbind(AIC=sapply(list(lmer.0,lmer.1,lmer.2,lmer.3,lmer.4), AIC.ml),
             BIC=sapply(list(lmer.0,lmer.1,lmer.2,lmer.3,lmer.4), BIC.ml),
             DIC=sapply(list(lmer.0,lmer.1,lmer.2,lmer.3,lmer.4), DIC.ml))
 
colnames(res) <- c("lmer.0","lmer.1", "lmer.2", "lmer.3","lmer.4")

t(round(res,2))

```
 
\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
 
```{r}
display(lmer.5 <- lmer(sqrt.CD4PCT ~ 1 + visage + treatmnt + (0+visage|newpid),data=cd4))

coef.5 <- coef(lmer.5)$newpid[,1:3]  ## coef.5[,1] = beta0         (fixed intcpt)
                                     ## coef.5[,2] = beta1 + eta1  (random slope on time)
                                     ## coef.5[,3] = beta2         (coef on treatmnt)
                                  
tx <- with(cd4,sapply(split(treatmnt,newpid), function(x) x[1]))  ## one per child

params <- data.frame(newpid=sort(unique(cd4$newpid)),
                     beta0=coef.5[,1],
                     alpha1=coef.5[,2],
                     beta2=coef.5[,3],
                     tx 
                     )

h + geom_abline(data=params[params$newpid<=12,],aes(intercept=beta0 + ifelse(tx==2,beta2,0),
                                                    slope=alpha1,
                                                    color=tx))


```
 

