---
title: "Homework 07 Solutions"
date: "2022-11-04"
output: 
  pdf_document:
    keep_tex: true
    fig_caption: true
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \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(lme4)
library(ggplot2)
```


# Problem 1.

## Problem 1a.

On July 23, 2006, the Houston Chronicle published an article entitled "Reading:
First-grade standard too tough for many". The article claimed in part that "more
students (across Texas) are having to repeat first grade. Experts attribute the
increase partially to an increase in poverty."  The data in the file  `houston.txt` is derived from this study.  Here is a quick look:
```{r}
houston <- read.table("houston.txt",header=T)

attach(houston)
houston$y <- Pct.Repeating.1st.Grade
houston$x <- Pct.Low.income.students
houston$Year <- factor(Year)
detach()

summary(houston)
```

\textit{Here, I use ordinary linear models to this data with} `lm()` \textit{to answer each question.  There is code below to make plots, but they are not very informative, so I commented them out to save space.  Basically the residuals for each of the models are a bit right-skewed, but not enough to motivate me to make a transformation and complicate interpretation of the models.}

## Problem 1a.

Is an increase in the percentage of low income students associated with an
increase in the percentage of students repeating first grade?

\textit{We just have to regress} `y` \textit{on} `x` \textit{and see if the coefficient on} `x` \textit{is significantly different from zero.}
```{r}
summary(lm.1 <- lm(y ~ x, data=houston))
## par(mfrow=c(2,2))
## plot(lm.1)
```
\textit{From the summary we see that the coefficient on $x$ is significantly different from zero, which suggests that there is an association between the percentage of low income students and the percentage of students repeating first grade.}



## Problem 1b.

Has there been an increase in the percentage of students repeating first grade
between 1994--1995 and 2004–-2005?

\textit{Similar to part (a), we regress on the categorical variable} `Year` \textit{and check for a significant coefficient.}
```{r}
summary(lm.2 <- lm(y ~ Year, data=houston))
## par(mfrow=c(2,2))
## plot(lm.2)

```
\textit{The coefficient on the} `Year2004` \textit{dummy variable is not significantly different from zero, so it does not look like there is a real increase in the percentage of students repeating first grade
between 1994--1995 and 2004–-2005.}

## Problem 1c.

Is the association (if any) between the percentage of students repeating first grade and the
percentage of low-income students different between 1994--1995 and 2004--2005?

\textit{The association is quantified by the slope of the regression line between} `y` \textit{and} `x`\textit{.  To see if the association changes, we should see if there is a significant difference in the slopes of the regression lines for the two years. We can do this by putting an interaction between} `Year` \textit{and} `x` \textit{in the model, and checking to see if the coefficient on the interaction is significantly different from zero (since that coefficient will be the difference in slopes between the two years.)}
```{r}
summary(lm.3 <- lm(y ~ x*Year, data=houston))
## par(mfrow=c(2,2))
## plot(lm.3)

```
\textit{We can see from the summary that the coefficient for the interaction is not significantly different from zero, so the slope doesn't change significantly from 1994 to 2004.  Actually the main effect for} `Year2004` \textit{is not significant either, so it doesn't even look like the intercept changes between years.}


# Problem 2

In the data for problem 1, it is plausible that there would be some variation among Counties in percentage of students repeating first grade.  We will use a multilevel model to explore the variation across counties, \textit{using the 2004 data \underline{\textbf{only}}}.

## Problem 2a.

Extract from the file `houston.txt` a data frame containing just the data from 2004.  Report the following for the extracted data frame:
```{r}
houston04 <- houston[houston$Year==2004,]

attach(houston04)
```
*  Number of groups, $J$.  (Each County is a group.)
```{r}
cat("J =",J <- length(unique(County)))
```
* For each group $j$, number of observations $n_j$.
```{r}
nj <- c(table(County))
cat("nj =\n")
nj
```
*  Total number of observations, $n$. (Each District is an observation.)
```{r}
cat("n =",n <- sum(nj))
detach()
```

Use only data from this data frame in the rest of this problem.

## Problem 2b.

Use `lmer()` from `library(lme4)` to fit the multilevel model
\begin{equation}
        \left.
        \begin{array}{rcl}
	y_i & = & \alpha_{j[i]} + \epsilon_i, 
		~\epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
            \alpha_j & = & \beta_0 + \eta_j, 
		~\eta_j \stackrel{iid}{\sim} N(0,\tau^2)
        \end{array}
        \right\}
\end{equation}
where $i$ indexes Districts, $j$ indexes Counties, and $y_i$ is the percentage of students repeating first grade in the $i^{th}$ District.  Provide: a `summary()` of your model, and report the following quantities from your summary:
```{r}
summary(lmer.1 <- lmer(y ~ 1 + (1|County), data=houston04))
```

\textit{We can read these values off the summary:}

* $\hat\sigma^2=14.068$;

* $\hat\tau^2=1.028$;

* $\hat\beta_0=6.7283$.

## Problem 2c.  

Use `library(ggplot2)` (or another R package if you prefer) to make a facet plot of $y$=Pct.Repeating.1st.Grade vs $x$=Pct.Low.income.students, with a different facet for each County.  Add to the facets

* A black horizontal line indicating the average of $y$ in the whole data set (ignoring groups)
* A red horizontal line indicating the average of $y$ in each County
* A green horizontal line indicating the "random effect" estimates $\hat\alpha_j = \hat\eta_j + \hat\beta_0$ for each County, as indicated in equation (1). 

and submit your plot.  (Does your plot illustrate the "shrinkage" phenomenon of MLM's?)

\textit{Here's the code to make the plot...}
```{r}
params <- with(houston04,
                 data.frame(
                    sort(unique(County)),
                    rep(mean(y),J),
                    0,
                    sapply(split(y,County),mean),
                    0,
                    coef(lmer.1)$County[,1], ## gives \beta_0 + \eta_j
                    0
                 )
          )
names(params) = c("County","int0","slo0","int1","slo1","int2","slo2")

g <- ggplot(houston04,aes(x=x,y=y)) + 
  facet_wrap( ~ County) + 
  geom_point()

g + geom_abline(data=params,aes(intercept=int0,slope=slo0),color="black") +
  geom_abline(data=params,aes(intercept=int1,slope=slo1),color="red") +
  geom_abline(data=params,aes(intercept=int2,slope=slo2),color="green") 
```

\textit{The facet plot does indeed show the "shrinkage" phenomenon: The fitted county means from the MLM (the green lines) are between the grand mean (the black lines) that would be estimated from} `lm(y ~ 1)` \textit{and the raw county means (the red lines) that would be estimated from} `lm(y ~ County)`.

## Problem 2d.

Fit the ordinary linear model `lm(y ~ County)`, where the coefficients on the County dummy variables are constrained to sum to zero.  Provide the `summary()` of your model, and decide from the output whether we need to keep the `County` variable in the model (this is also an informal way to decide whether to keep the multilevel model that you fitted in part (a)).
```{r}
cty <- as.factor(houston04$County)
contrasts(cty) <- contr.sum(J)

summary(lm.1.04 <- lm(y ~ cty,data=houston04))

```

\textit{Only one of the} `County` \textit{coefficients is even close to significantly different from zero.}


\textit{The F statistic, 1.642 on 7 and 53 df, tests the fitted model against the intercept-only model} `lm(y ~ 1)` \textit{The p-value is  0.1439, which suggests there is not enough evidence in the data to reject the intercept-only model and keep} `County` \textit{in the model.}

\textit{We would get exactly the same result comparing the models directly with the} `anova()` \textit{function:}

```{r}
lm.0.04 <- lm(y ~ 1, data=houston04)
anova(lm.0.04,lm.1.04)

```


# Problem 3.

\it

In the multilevel model $(*)$ for data $y_i$, $i=1,\ldots, n$, arranged into $J$ groups, $j=1,\ldots, J$, where
each group $j$ has $n_j$ observations,
$$
        \left.
        \begin{array}{rcl}
	y_i & = & \alpha_{j[i]} + \epsilon_i, 
		~\epsilon_i \stackrel{iid}{\sim} N(0,\sigma^2)\\
            \alpha_j & = & \beta_0 + \eta_j, 
		~\eta_j \stackrel{iid}{\sim} N(0,\tau^2)
        \end{array}
        \right\}~,
        \tag{*}
$$
we also assume that  the $\epsilon$'s and $\eta$'s are independent of each other.  

## Problem 3a.

We calculate
\begin{align}
  \Cov(y_i, y_{i'}) &= \Cov(\beta_0 + \eta_{j[i]} + \epsilon_i, \beta_0 + \eta_{j[i']} + \epsilon_{i'}) \\
  &= \Cov(\eta_{j[i]} + \epsilon_i,  \eta_{j[i']} + \epsilon_{i'}) 
  \label{e:no-betas} \\                
  &= \Cov(\eta_{j[i]}, \eta_{j[i']}) + \Cov(\eta_{j[i]}, \epsilon_{i'}) + \Cov(\epsilon_i, \eta_{j[i']}) + \Cov(\epsilon_i, \epsilon_{i'}) \label{eq:1a} \\
  &= \Cov(\eta_{j[i]}, \eta_{j[i']}) + 0 + 0 + 0 \label{e:zeros}
\end{align}                               
where \eqref{e:no-betas} follows because adding the constant $\beta_0$ does not change the covariance, and \eqref{e:zeros} follows because the $\epsilon$'s are independent of each other and the $\eta$'s are independent of the $\epsilon$'s.

If $j[i]\neq j[i']$ then $\Cov(\eta_{j[i]}, \eta_{j[i']})=0$ too, because different $\eta$'s are also independent of each other, and so $\Corr(\eta_{j[i]}, \eta_{j[i']})=0$ also.

## Problem 3b.

If $j[i]=j[i']$ in \eqref{eq:1a}, then the $\eta$'s in \eqref{e:zeros} are the same and so
\begin{align}
  \Cov(\eta_{j[i]}, \eta_{j[i']}) + 0 + 0 + 0 
  &= \Var(\eta_{j[i]}) + 0 + 0 + 0 \\
  &= \tau^2\label{eq:1b}
\end{align}

Finally, when $j[i]=j[i']$,
\begin{align}
  \Corr(y_i, y_{i'}) 
  &= \frac{\Cov(y_i, y_{i'})}{\sqrt{\sigma^2_{y_i}\sigma^2_{y_{i'}}}} \\
  &= \frac{\tau^2}{\sqrt{(\sigma^2 + \tau^2)(\sigma^2 + \tau^2)}} \label{e:almost} \\
  &= \frac{\tau^2}{\sigma^2 + \tau^2}
\end{align}
where the numerator in \eqref{e:almost} comes from \eqref{eq:1b} and the terms in the denominator follows from 
\begin{align}
   \Var(y_i) &= \Var(\beta_0 + \eta_j[i] + \epsilon_i) \\
   &= 0 + \tau^2 + \sigma^2
\end{align}


## Problem 3c.
\def\i{_{\{i:j[i]=j\}}}
\begin{align}
  \frac{1}{n_j}\sum\i Y_i &= \frac{1}{n_j}\left[\sum\i\beta_0 + \sum\i\eta_j + \sum\i \epsilon_i\right] \\
                        &= \beta_0 + \eta_j + \frac{1}{n_j}\sum\i\epsilon_i 
\end{align}
which means
\begin{align}
   \Var(\bar{y}_j) &= \Var(\eta_j) + \frac{1}{n_j^2}\sum\i\Var(\epsilon_i) \\
          &= \tau^2 + \frac{\sigma^2}{n_j} \label{e:bar-var}
\end{align}

## Problem 3d.

Using the same kinds of ideas as above, we can calculate
\begin{align}
  \Var(\bar{y}_j) &= \tau^2 + \frac{\sigma^2}{n_j} \mbox{~~(from \eqref{e:bar-var})} \\
  \Var(\bar{y}_j^*) &= \tau^2 + \frac{\sigma^2}{n_j} \mbox{~~(from \eqref{e:bar-var} again)} \\
  \Cov(\bar{y}_j,\bar{y}_j^*) &= \Cov\left(\beta_0 + \eta_j + \frac{1}{n_j}\sum\i\epsilon_i, \beta_0 + \eta_j + \frac{1}{n_j}\sum\i\epsilon_i^*\right) \\
  &=\tau^2 \mbox{~~(similarly to \eqref{eq:1b})}
\end{align}
and therefore
\begin{align}
\Corr(\bar{y}_j,\bar{y}_j^*) &=  \frac{\Cov(\bar{y}_j,\bar{y}_j^*)}{\sqrt{\Var(\bar{y}_j)\Var(\bar{y}_j^*)}} \\
  &= \frac{\tau^2}{\sqrt{(\tau^2 + \sigma^2/n_j)(\tau^2 + \sigma^2/n_j)}}\\
  &= \frac{\tau^2}{\tau^2 + \sigma^2/n_j}
\end{align}


\rm