\documentclass[]{article}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
\else % if luatex or xelatex
  \ifxetex
    \usepackage{mathspec}
  \else
    \usepackage{fontspec}
  \fi
  \defaultfontfeatures{Ligatures=TeX,Scale=MatchLowercase}
\fi
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
% use microtype if available
\IfFileExists{microtype.sty}{%
\usepackage{microtype}
\UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\usepackage[margin=1in]{geometry}
\usepackage{hyperref}
\hypersetup{unicode=true,
            pdftitle={Homework 09 Solutions},
            pdfborder={0 0 0},
            breaklinks=true}
\urlstyle{same}  % don't use monospace font for urls
\usepackage{graphicx,grffile}
\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}
\IfFileExists{parskip.sty}{%
\usepackage{parskip}
}{% else
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
}
\setlength{\emergencystretch}{3em}  % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{0}
% Redefines (sub)paragraphs to behave more like sections
\ifx\paragraph\undefined\else
\let\oldparagraph\paragraph
\renewcommand{\paragraph}[1]{\oldparagraph{#1}\mbox{}}
\fi
\ifx\subparagraph\undefined\else
\let\oldsubparagraph\subparagraph
\renewcommand{\subparagraph}[1]{\oldsubparagraph{#1}\mbox{}}
\fi

%%% Use protect on footnotes to avoid problems with footnotes in titles
\let\rmarkdownfootnote\footnote%
\def\footnote{\protect\rmarkdownfootnote}

%%% Change title format to be more compact
\usepackage{titling}

% Create subtitle command for use in maketitle
\providecommand{\subtitle}[1]{
  \posttitle{
    \begin{center}\large#1\end{center}
    }
}

\setlength{\droptitle}{-2em}

  \title{Homework 09 Solutions}
    \pretitle{\vspace{\droptitle}\centering\huge}
  \posttitle{\par}
    \author{}
    \preauthor{}\postauthor{}
      \predate{\centering\large\emph}
  \postdate{\par}
    \date{12/01/2019}

\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{MnSymbol}
\DeclareMathOperator{\Cov}{Cov}
\DeclareMathOperator{\Corr}{Corr}
\DeclareMathOperator{\Var}{Var}

\begin{document}
\maketitle

\section{2}\label{section}

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

\subsection{(a)}\label{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.

\subsection{(b)}\label{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}

\subsection{(c)}\label{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}

\subsection{(d)}\label{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}

\section{3}\label{section-1}

\subsection{(a)}\label{a-1}

\begin{verbatim}
## Linear mixed model fit by REML ['lmerMod']
## Formula: per.cap.income ~ pct.hs.grad + (pct.hs.grad | state)
##    Data: dat
## 
## REML criterion at convergence: 8293
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.146 -0.583 -0.085  0.454  4.321 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  state    (Intercept) 7949288  2819.4        
##           pct.hs.grad    3676    60.6   -1.00
##  Residual             8115992  2848.9        
## Number of obs: 440, groups:  state, 48
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  -4920.7     1699.0    -2.9
## pct.hs.grad    297.2       23.1    12.9
## 
## Correlation of Fixed Effects:
##             (Intr)
## pct.hs.grad -0.984
## convergence code: 0
## boundary (singular) fit: see ?isSingular
\end{verbatim}

\begin{verbatim}
##        grp        var1        var2    vcov sdcor
## 1    state (Intercept)        <NA> 7949288  2819
## 2    state pct.hs.grad        <NA>    3676    61
## 3    state (Intercept) pct.hs.grad -170937    -1
## 4 Residual        <NA>        <NA> 8115992  2849
\end{verbatim}

The parameter estimates are as follows.

\textbf{Random effects estimates}

\begin{align*}
    \widehat{\sigma^2} &= 8,115,992 \\
    \widehat{\tau_0^2} &= 7,949,288 \\
    \widehat{\tau_1^2} &= 3,676 \\
    \widehat{\mathrm{Corr}}(\eta_0,\eta_1) &= -0.99999
\end{align*}

\textbf{Fixed effects estimates}

\begin{align*}
    \widehat{\beta_0} &= -4,921 \\
    \mathrm{SE}\left(\widehat{\beta_0}\right) &= 1,699 \\
    \widehat{\beta_1} &= 297.23 \\
    \mathrm{SE}\left(\widehat{\beta_1}\right) &= 23.05 \\
\end{align*}

If we fit the model calling \texttt{lmer} and setting the option within
\texttt{control\ =\ lmerControl(optimizer\ =\ \textquotesingle{}bobyqa\textquotesingle{})}
we get a similar but slightly different model.

\begin{verbatim}
## Linear mixed model fit by REML ['lmerMod']
## Formula: per.cap.income ~ pct.hs.grad + (pct.hs.grad | state)
##    Data: dat
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 8292
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -3.053 -0.575 -0.096  0.453  4.244 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  state    (Intercept) 23008246 4797          
##           pct.hs.grad     7397   86     -1.00
##  Residual              8064024 2840          
## Number of obs: 440, groups:  state, 48
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  -4455.3     1791.2   -2.49
## pct.hs.grad    291.3       24.8   11.73
## 
## Correlation of Fixed Effects:
##             (Intr)
## pct.hs.grad -0.987
## convergence code: 0
## boundary (singular) fit: see ?isSingular
\end{verbatim}

\subsection{(b)}\label{b-1}

\begin{figure}
\centering
\includegraphics{hw09_solutions_files/figure-latex/3b-1.pdf}
\caption{Problem 3b: HS graduation rate vs.~county level per capita
income with regression lines.\label{f:3b_plot}}
\end{figure}

See Figure \ref{f:3b_plot}.

\subsection{(c)}\label{c-1}

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.

\subsection{(d)}\label{d-1}

\includegraphics{hw09_solutions_files/figure-latex/3d-1.pdf}
\includegraphics{hw09_solutions_files/figure-latex/3d-2.pdf}

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.

\subsection{(e)}\label{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 \texttt{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.

\begin{figure}
\centering
\includegraphics{hw09_solutions_files/figure-latex/3e-1.pdf}
\caption{Problem 3e: QQ plot of multilevel model
residuals.\label{f:3e_qq}}
\end{figure}

\begin{figure}
\centering
\includegraphics{hw09_solutions_files/figure-latex/3e_marg-1.pdf}
\caption{Problem 3e: Facet plot of marginal residuals.\label{f:3e_marg}}
\end{figure}

\begin{figure}
\centering
\includegraphics{hw09_solutions_files/figure-latex/3e_cond-1.pdf}
\caption{Problem 3e: Facet plot of conditional
residuals.\label{f:3e_cond}}
\end{figure}

\begin{figure}
\centering
\includegraphics{hw09_solutions_files/figure-latex/3e_ranef-1.pdf}
\caption{Problem 3e: Facet plot of random effects. Blue line shows mean
state random effect.\label{f:3e_ranef}}
\end{figure}


\end{document}
