---
title: "Homework 09 Solutions"
date: "12/01/2019"
output: 
  pdf_document:
    keep_tex: true
fig_caption: true
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \usepackage{MnSymbol}
  - \DeclareMathOperator{\Cov}{Cov}
  - \DeclareMathOperator{\Corr}{Corr}
  - \DeclareMathOperator{\Var}{Var}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(cache = FALSE, echo = FALSE, include = TRUE, message = FALSE)
options(scipen = 1, digits = 2)

library(dplyr)
library(ggplot2)
library(lme4)
library(pander)
library(tidyr)
library(scales)

source("residual-functions.r")
```

# 2
We need to assume throughout this problem that the errors $\epsilon_i$ are independent of the quantities $\eta_j$. (Remember, $\epsilon_i$ represents the "unexplainable variation"; see slide 14 of lecture 23.)

## (a)
Since, for two random variables $A, B$, we have that $\Corr(A, B) = 0 \iff \Cov(A, B) = 0$, let's just focus on the covariance. We have:
\begin{align}
  \Cov(y_i, y_{i'}) &= \Cov(\alpha_{j[i]} + \epsilon_i, \alpha_{j[i']} + \epsilon_{i'}) \\
  &= \Cov(\alpha_{j[i]}, \alpha_{j[i']}) + \Cov(\alpha_{j[i]}, \epsilon_{i'}) + \Cov(\epsilon_i, \alpha_{j[i']}) + \Cov(\epsilon_i, \epsilon_{i'}) \label{eq:1a}
\end{align}

Recall that if two random variables $A$ and $B$ are independent, then $f(A)$ and $g(B)$ are independent for any functions $f$ and $g$\footnote{Any measurable functions $f$ and $g$, that is, but you don't really need to worry about this technicality.}. Since $\eta_{j[i]} \upmodels \eta_{j[i']}$, and $\beta_0$ is just a constant, it follows that $\alpha_{j[i]} \upmodels \alpha_{j[i']}$

Furthermore, we have the errors $\epsilon_i$ are independent of everything, by assumption.

Hence each pair of variables in \eqref{eq:1a} is independent, so the covariances are all 0.

## (b)
Substituting $j[i]$ for $j[i']$ in \eqref{eq:1a}, we have 
\begin{align}
  \Cov(y_i, y_{i'}) &= \Cov(\alpha_{j[i]}, \alpha_{j[i]}) + \Cov(\alpha_{j[i]}, \epsilon_{i'}) + \Cov(\epsilon_i, \alpha_{j[i]}) + \Cov(\epsilon_i, \epsilon_{i'}) \\
  &= \Var(\alpha_{j[i]}) + 0 + 0 + 0 \\
  &= \tau^2
\end{align}
Now,
\begin{align}
  \Corr(y_i, y_{i'}) &= \frac{\Cov(y_i, y_{i'})}{\sqrt{\sigma^2_{y_i}\sigma^2_{y_{i'}}}} \label{eq:1b}
\end{align}
Writing the model in variance components form, we have $y_i = \beta_0 + \eta_{j[i]} + \epsilon_i$, which means $y_i \sim N(\beta_0, \sigma^2 + \tau^2)$ (since the sum of two independent normally distributed random variables has variance equal to the sum of the variances). Substituting this into \eqref{eq:1b} yields
\begin{align}
  \Corr(y_i, y_{i'}) &= \frac{\Cov(y_i, y_{i'})}{\sqrt{(\sigma^2 + \tau^2)(\sigma^2 + \tau^2)}} \\
  &= \frac{\tau^2}{\sigma^2 + \tau^2}
\end{align}

## (c)
\begin{align}
  \frac{1}{n_j}\sum Y_i &= \frac{1}{n_j}[\sum\alpha_j + \sum \epsilon_i] \\
                        &= \alpha_j + \frac{1}{n_j}\sum\epsilon_i \\
          \implies \Var(\bar{y}_j) &= \Var(\alpha_j) + \frac{1}{n_j^2}\sum\Var(\epsilon_i) \\
          &= \tau^2 + \frac{\sigma^2}{n_j}
\end{align}

## (d)
Notice that this is equivalent to sampling a group, sampling a bunch of data points from that group, splitting those data points into two sub-groups, and taking the means of each subgroup. We have:
\begin{align}
  \Cov(\bar{y}_j, \bar{y}_j^*) &= \Cov\left(\frac{1}{n_j}\sum_i y_i,  \frac{1}{n_j^*}\sum_k y_k^*\right) \\
  &= \frac{1}{n_j^2}\sum_i\sum_k\Cov(y_i, y_k^*) \\
  &= \frac{1}{n_j^2}n_j^2\tau^2 \\
  &= \tau^2
\end{align}
where in line 2 we assumed that $n_j^* = n_j$, and in line 3 we used that $\Cov(y_i, y_k^*)$ is equivalent to $\Cov(y_i, y_{i'})$ from part (b) above. Now, using the result from part (c), we have
\begin{align}
  \Cov(\bar{y}_j, \bar{y}_j^*) &= \frac{\tau^2}{\sqrt{\sigma^2_{\bar{y}}\sigma^2_{\bar{y}^*}}} \\
  &= \frac{\tau^2}{\tau^2 + \sigma^2/n_j}
\end{align}

# 3
## (a)

```{r 3a}
dat <- read.table("data/cdi.dat")
fit1 <- lmer(per.cap.income ~ pct.hs.grad + (pct.hs.grad | state), data = dat)
print(s1 <- summary(fit1))

beta0 <- s1[["coefficients"]][1, 1]
beta1 <- s1[["coefficients"]][2, 1]
se.beta0 <- s1[["coefficients"]][1, 2]
se.beta1 <- s1[["coefficients"]][2, 2]

tau0 <- VarCorr(fit1)[[1]][1, 1]
tau1 <- VarCorr(fit1)[[1]][2, 2]

as.data.frame(VarCorr(fit1))
sigma <- as.data.frame(VarCorr(fit1))[4, 5]^2
rho <- as.data.frame(VarCorr(fit1))[3, 5]
```

The parameter estimates are as follows.

**Random effects estimates**
\begin{align*}
    \widehat{\sigma^2} &= `r format(sigma, scientific = FALSE, big.mark = ",")` \\
    \widehat{\tau_0^2} &= `r format(tau0, scientific = FALSE, big.mark = ",")` \\
    \widehat{\tau_1^2} &= `r format(tau1, scientific = FALSE, big.mark = ",")` \\
    \widehat{\mathrm{Corr}}(\eta_0,\eta_1) &= `r format(rho, digits = 5)`
\end{align*}

**Fixed effects estimates**
\begin{align*}
    \widehat{\beta_0} &= `r format(beta0, scientific = FALSE, big.mark = ",")` \\
    \mathrm{SE}\left(\widehat{\beta_0}\right) &= `r format(se.beta0, scientific = FALSE, big.mark = ",")` \\
    \widehat{\beta_1} &= `r beta1` \\
    \mathrm{SE}\left(\widehat{\beta_1}\right) &= `r se.beta1` \\
\end{align*}

If we fit the model calling `lmer` and setting the option within `control = lmerControl(optimizer = 'bobyqa')` we get a similar but slightly different model.

```{r 3a2}
summary(lmer(per.cap.income ~ pct.hs.grad + (pct.hs.grad | state), data = dat,
     control = lmerControl(optimizer = 'bobyqa')))

```


## (b)

```{r 3b, fig.height=7.5, fig.cap="Problem 3b: HS graduation rate vs. county level per capita income with regression lines.\\label{f:3b_plot}", warning=FALSE}

ranef.states <- ranef(fit1)[["state"]]
names(ranef.states) <- c("intercept_ranef", "slope_ranef")
ranef.states$state <- row.names(ranef.states)

fit.pooled <- lm(per.cap.income ~ pct.hs.grad, data = dat)

state_lm <- function(df){
    mod <- lm(per.cap.income ~ pct.hs.grad, data = df)
    return(data.frame(intercept_unpooled = mod[["coefficients"]][1], slope_unpooled = mod[["coefficients"]][2]))
}
lm.df <- group_by(dat, state) %>% do(., state_lm(.))
lm.df$intercept_pooled <- coef(fit.pooled)[1]
lm.df$slope_pooled <- coef(fit.pooled)[2]
lm.df <- left_join(lm.df, ranef.states, by = "state")
lm.df$intercept_betas <- beta0
lm.df$slope_betas <- beta1
lm.df$slope_alphas <- lm.df$slope_ranef + lm.df$slope_betas
lm.df$intercept_alphas <- lm.df$intercept_ranef + lm.df$intercept_betas

lm.df <- gather(lm.df, value = est, key = coef, -state) %>%
    separate(coef, into = c("term", "model"), sep = "_") %>%
    spread(value = est, key = term) %>%
    filter(model != "ranef")

ggplot(dat) + geom_point(aes(pct.hs.grad, per.cap.income)) +
    geom_abline(data = lm.df, aes(intercept = intercept, slope = slope, color = model, group = model, linetype = model), alpha = 0.7) +
    facet_wrap(~state, ncol = 5) +
    scale_y_continuous(label = comma) +
    xlab("Percentage High School Graduates") + ylab("Per Capitca Income ($)") +
    theme_minimal() +
    theme(legend.position = "bottom")
```

See Figure \ref{f:3b_plot}.

## (c)

Examining Figure \ref{f:3b_plot} we see that the pooled regression and mlm $\beta$'s agree almost exactly, which makes sense given that they are both models considering all of the data.  There is a lot more variance in the unpooled regression lines since they are estimated individually for each state.  We also note that for states with only one observation we are unable to estimate the state level regression.  The mlm $\alpha$'s fall somewhere between the mlm $\beta$'s and the unpooled regression lines.  In states with very little data they agree with the mlm $\beta$'s (e.g. HI, DE) while in states with a lot of data (e.g. CT, NJ) they agree more closely with the unpooled regression line.

## (d)

```{r 3d, fig.width=3.5, fig.cap="Problem 3d: Variance-covariance plot for mlm from a.\\label{f:3d_plot}"}
par(mfrow = c(1, 2))
image(Sigma.y(fit1), main = "Unsorted")
image(Sigma.y(fit1)[order(dat$state),][ , order(dat$state)], main = "Sorted by State")
```

We plot the unsorted variance-covariance matrix in the left panel of Figure \ref{f:3d_plot} and the same matrix with rows and columns sorted by state in the right panel.  The right panel of Figure \ref{f:3d_plot} makes is clear that the estimated covariance between counties in the same state is high and zero (by assumption) between counties in different states. Larger squares correspond to states with more observed counties.

## (e)

QQ plots are shown in Figure \ref{f:3e_qq} while the marginal residuals, conditional residuals, and random effects are shown in Figures \ref{f:3e_marg}, \ref{f:3e_cond}, and \ref{f:3e_ranef} respectively.  The QQ plots all show a heavy right tail as well as a fewer large (in absolute value) residuals in the left tail.  The marginal residuals generally appear mean zero but suggest that the deviations from normality may be coming primarily from CA and NJ which do not appear mean zero when considered individually.  It is worth noting that these states both contain a lot of observations and have relatively large random effects which is what makes them so noticeable.  The conditional residuals appear good for the most part but there may be some trend in PA although it is hard to know how to interpret this since they are plotted against an index and not a covariate. The random effects appear fine with the residuals clustered around the state mean in all cases.

Given the deviation from normality observed in the qq plots we should consider a transformation of `per.cap.income`, perhaps a log transform, to try to address the heavy tails.  The conditional residuals suggest that we might be able to make some improvements in how the random effects are modeled but that these are likely to be marginal improvements affecting only a few states.

```{r 3e, fig.height=3, fig.cap="Problem 3e: QQ plot of multilevel model residuals.\\label{f:3e_qq}"}
resids <- data.frame(id = dat$id,
                     resid_marginal = r.marg(fit1)[, 1],
                     resid_conditional = r.cond(fit1),
                     resid_random = r.reff(fit1)[, 1],
                     resid_cholesky = r.chol(fit1),
                     yhat_marg = yhat.marg(fit1),
                     yhat_cond = yhat.cond(fit1),
                     yhat_rand = yhat.reff(fit1)) %>%
left_join(., dat, by = "id")

par(mfrow = c(1, 3))
qqnorm(resids$resid_marginal, main = "Marginal Residuals")
qqline(resids$resid_marginal)
qqnorm(resids$resid_conditional, main = "Conditional Residuals")
qqline(resids$resid_conditional)
qqnorm(resids$resid_random, main = "Random Residuals")
qqline(resids$resid_random)
```

```{r 3e_marg, fig.height=7.5, fig.cap="Problem 3e: Facet plot of marginal residuals.\\label{f:3e_marg}", warning=FALSE}
ggplot(resids, aes(id, resid_marginal)) +
    geom_point() +
    geom_hline(yintercept = 0) +
    facet_wrap(~state, ncol = 5) +
    scale_y_continuous(label = comma) +
    xlab("") + ylab("Marginal Residual") +
    theme_minimal()
```

```{r 3e_cond, fig.height=7.5, fig.cap="Problem 3e: Facet plot of conditional residuals.\\label{f:3e_cond}", warning=FALSE}
ggplot(resids, aes(id, resid_conditional)) +
    geom_point() +
    geom_hline(yintercept = 0) +
    facet_wrap(~state, ncol = 5) +
    scale_y_continuous(label = comma) +
    xlab("") + ylab("Conditional Residuals") +
    theme_minimal()
```

```{r 3e_ranef, fig.height=7.5, fig.cap="Problem 3e: Facet plot of random effects. Blue line shows mean state random effect.\\label{f:3e_ranef}", warning=FALSE}
state.means <- group_by(resids, state) %>% summarise_all(., mean)

ggplot(resids, aes(id, resid_random)) +
    geom_point() +
    geom_hline(yintercept = 0) +
    geom_hline(data = state.means, aes(yintercept = resid_random), color = "blue") +
    facet_wrap(~state, ncol = 5) +
    scale_y_continuous(label = comma) +
    xlab("") + ylab("Random Effects") +
    theme_minimal()
```
