---
title: "Homework 02 Solutions"
date: "2022-09-08"
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)
knitr::opts_chunk$set(fig.width=6, fig.height=5, fig.align="center") 
# library(knitr)
```


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

\frenchspacing


\begin{center}
\Large
36-617: Applied Linear Models\\
Fall 2022\\
Solutions
\end{center}


```{r}

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



# Problem 1  

Let $y = X\beta + \epsilon$, where $y=(y_1,\ldots, y_n)^T$ is an $n\times 1$ column vector, $X$ is an $n\times (p+1)$ matrix whose first column is all 1's, $\beta=(\beta_0,\ldots,\beta_p)^T$ is a $(p+1)\times 1$ column vector, and $\epsilon=(\epsilon_1,\ldots,\epsilon_n)^T \sim N(0,\sigma^2I)$ is an $n\times 1$ random column vector, following a multivariate Normal distribution with mean vector 0 and variance-covariance matrix $\sigma^2I$, where $I$ is the $n\times n$ identity matrix.

## Problem 1(a)

Let 
\[
     U = \left[\begin{array}{c}
               1 \\ 1 \\ \vdots \\ 1
               \end{array}\right]
\]
be an $n\times 1$ column of all 1's.  

\renewcommand{\labelenumi}{\roman{enumi}}

\begin{enumerate}
\item Using the definition of the hat matrix $H$, show $HX=X$
\item Using your result in (i), show $HU=U$
\end{enumerate}

\em

### Part (i)

$HX = [X(X^TX)^{-1}X^T]X = X(X^TX)^{-1}(X^TX) = X$

### Part (ii)

Two ways to think about it:
\begin{itemize}
\item \textbf{In words:} Since $HX=X$ we know that all of the columns of $X$ are unchanged when we premultiply $X$ by $H$.  In particular, since the first column of $X$ is $U$, we know $U$ will be unchanged when we premultiply by $H$: $HU=U$.


\item \textbf{In math:}  Let 
\[
      E_1 = \left[\begin{array}{c}
             1 \\ 0 \\ 0 \\ \vdots \\ 0 \end{array}\right] ,
\]
the $n\times 1$ column vector with a 1 in the first row and 0's in the remaining $n-1$ rows.
Since the first column of $X$ is $U$, we know $XE_1 = U$.  Therefore,
\[
     HU = H(XE_1) = (HX)E_1 = XE_1 = U
\]


\end{itemize}

\rm

## Problem 1(b) 

Use properties of the hat matrix $H=X(X^TX)^{-1}X^T$ to show that the column vector
\[
    y - \hat y = \left[\begin{array}{c}
                       y_1 - \hat y_1 \\
                       y_2 - \hat y_2 \\
                       \vdots \\
                       y_n - \hat y_n
                       \end{array}\right] 
 \]
can be written as $y - \hat y = (I-H)y$.
 
\em
Since $Hy= \hat y$ and $Iy = y$,
\[
    y - \hat y = Iy - Hy = (I-H)y
\]

\rm

## Problem 1(c)

For any $n\times1$ column vector
\[ 
     a =\left[\begin{array}{c}
     a_1 \\ a_2 \\ \vdots \\ a_n
     \end{array}\right] 
\]
we know that $a^TU=\sum_{i=1}^n a_i$.  Use this fact, together with your results in (a) and (b), to show that 
\[
\sum_{i=1}^n (y_i - \hat y_i) = 0 .
\]  


\em

\[
     \sum_{i=1}^n (y_i - \hat y_i) = \underbrace{(y - \hat y)^T}_{1\times n}\underbrace{U}_{n\times 1} = \underbrace{[(I-H)y]^T}_{1\times n}\underbrace{U}_{n\times 1} = \underbrace{y^T}_{1\times n}\underbrace{(I-H)}_{n\times n}\underbrace{U}_{n\times 1} = \underbrace{y^T}_{1\times n}\underbrace{(U-U)}_{n\times 1} = \underbrace{Y^T}_{1\times n}\underbrace{0}_{n\times 1} = \underbrace{0}_{1\times 1}
\]


\rm

\textrm{\textbf{Aside:} Hence $y$ and $\hat y$ have the same sample mean $\bar y$, and the residuals vs fitted plot should always be ``balanced'', in some sense, around the horizontal line at $\hat e = 0$.}


# Problem 2

Let $y = X\beta + \epsilon$, where $y=(y_1,\ldots, y_n)^T$ is an $n\times 1$ column vector, $X$ is an $n\times (p+1)$ matrix whose first column is all 1's, $\beta=(\beta_0,\ldots,\beta_p)^T$ is a $(p+1)\times 1$ column vector, and $\epsilon=(\epsilon_1,\ldots,\epsilon_n)^T \sim N(0,\sigma^2I)$ is an $n\times 1$ random column vector, following a multivariate Normal distribution with mean vector 0 and variance-covariance matrix $\sigma^2I$, where $I$ is the $n\times n$ identity matrix.

## Problem 2(a)

Use properties of the hat matrix $H=X(X^TX)^{-1}X^T$ and the multivariate Normal distribution as discussed in class, to show\footnote{Very similar to problem 1(b).}
\[
			\hat e \sim N(0,(I-H)\sigma^2)
\]
where $\hat e$ is the column vector $\hat e = y - \hat y$.

\em

We will use several facts from lecture. First, let $W \sim N(\mu, \Sigma)$ be an arbitrary multivariate normal random vector and $A$ a matrix whose second dimension is the same as the length of $W$, so that $AW$ is defined. Then we have
\begin{align}
  & AW \sim N(A\mu, A\Sigma A^T) \label{eq:linear_normal}
\end{align}
Additionally, we have from lecture that
\begin{align}
               & HX\beta = X\beta \label{eq:idempotent1} \\
  & (I - H)(I - H)^T = (I - H)(I - H) = (I - H) \label{eq:idempotent2} \\
                 & y \sim N(X\beta, \sigma^2 I)  \label{eq:y_dist}
\end{align}
Let $\hat{y}$ be the $n \times 1$ column vector of fitted values.
\begin{align*}
  \hat{e} &= y - \hat{y} \\
          &= y - Hy \\
          &= (I - H)y \\
          &\sim N((I - H)X\beta, (I-H)\sigma^2 I (I-H)^T) \tag{by facts \eqref{eq:linear_normal} and \eqref{eq:y_dist}} \\
          &= N(0, \sigma^2 (I - H)(I - H)^T) \tag{by fact \eqref{eq:idempotent1}} \\
          &= N(0, \sigma^2 (I - H)) \tag{by fact \eqref{eq:idempotent2}}
\end{align*} 


\rm


## Problem 2(b) 

Let $H$ be the hat matrix for the multivariate regression model $y=X\beta + \epsilon$ as in part (a), and let $H_1$ be the hat matrix for the intercept-only model 
$y=\beta_0  + \epsilon$. 

\renewcommand{\labelenumi}{\roman{enumi}}

\begin{enumerate}
\item Show that the fitted values $\hat y $ for the intercept-only model is an $n\times 1$ column vector, all of whose entries are $\bar y$, that is, 
$$
			\hat y  = \left[\begin{array}{c} \hat y _1 \\ \vdots \\ \hat y _n\end{array}\right]
			= \left[\begin{array}{c} \bar y \\ \vdots \\ \bar y\end{array}\right]
			\eqno{(*)}
$$
(where the first ``='' is the definition of $\hat y $ and the second ``='' is what I want you to show).

\item Find a simple expression, in terms of (some or all of) $y$, $I$, $H$ and $H_1$, for the sample covariance
\[
    \mbox{Cov}(y,\hat y) = \frac1n {\sum}_{i=1}^n (y_i - \bar y)(\hat y_i - \bar y).
\]
\textit{(\textbf{ Hint:} We can rewrite $\mbox{Cov}(y,\hat y) = \frac1n (y - \bar y)^T(\hat y - \bar y)$, where $\hat y$ is the column vector of fitted values from $y=X\beta + \epsilon$ and, abusing notation slightly, $\bar y$ is the column vector in $(*)$ above, i.e. the fitted values from the intercept-only model $y = \beta_0  + \epsilon $.)}
\item Continue along the lines of the calculations in part (ii) to show that the sample correlation between $y$ and $\hat y$ can be written as
\[
	    \mbox{Corr}(y,\hat y) = \sqrt{\frac{SS_{reg}}{SST}}
\] 
and hence $R^2$ for the regression model $y = X\beta + \epsilon$ really is the squared correlation between $y$ and $\hat y$:
\[
			R^2 = \mbox{Corr}(y,\hat y)^2~.
\]

\end{enumerate}

\em

### Part (i)
In the intercept-only model, the design matrix $X$ is an $n\times 1$ matrix in
which each entry is 1. This can be used to derive the hat matrix $H_1$:
\begin{align*}
  X^TX &= [n]_{1\times 1} \tag{the $1\times 1$ matrix whose only entry is $n$} \\
  \implies (X^TX)^{-1} &= [1/n]_{1\times 1} \tag{the $1\times 1$ matrix whose only entry is $1/n$} \\
  \implies X(X^TX)^{-1}X^T &= \frac{1}{n}XX^T \\
                           &= \frac{1}{n} \begin{bmatrix}
                           1 & \ldots & 1 \\
                           \vdots & \ddots & \vdots \\
                           1 & \ldots & 1 \end{bmatrix}_{n\times n} \\
                           &= H_1 \tag{by the definition of the hat matrix}
\end{align*}
For any $i$ in $1, 2, \ldots, n$, we derive $\hat{y}_i$ by multiplying the $i$th row of $H_1$ by the vector $y$. In other words, for the intercept-only model, $\hat{y}_i = \frac{1}{n}\sum_{i=1}^n y_i = \bar{y}$ for all $i$. 

### Part (ii)
To avoid using $\bar{y}$ to denote both a scalar and a vector, let's immediately use the fact that $H_1y$ is equal to the $n\times 1$ vector $\begin{bmatrix} \bar{y} & \bar{y} & \ldots & \bar{y}\end{bmatrix}^T$, as shown in the previous problem. Rewriting the covariance in vector form, per the problem suggestion, we have:
\begin{align*}
  \Cov(y, \hat{y}) &= \frac1n {\sum}_{i=1}^n (y_i - \bar y)(\hat y_i - \bar y) \\
                   &= \frac{1}{n}(y - H_1y)^T(\hat{y} - H_1y) \\
                   &= \frac{1}{n}(y - H_1y)^T(Hy - H_1y) \\
                   &= \frac{1}{n}[(I - H_1)y]^T(H - H_1)y \\
                   &= \frac{1}{n}y^T(I - H_1)^T(H - H_1)y \tag{since $(AB)^T = B^T A^T$} \\
                   &= \frac{1}{n}y^T(H - H_1 - H_1^TH + H_1^TH_1)y \\
                   &= \frac{1}{n}y^T(H - H_1H)y \tag{since $H_1^T = H_1$ and $H_1H_1 = H_1$} \\
                   &= \frac{1}{n}y^T(H - H_1)y \tag{since $H_1H = HH_1 = H_1$}
\end{align*} 

### Part (iii)
Recall that the sample correlation is the standardized sample covariance, defined as
\begin{align*}
  \\widehat{\Corr}(y, \hat{y}) &= \frac{\hat{\Cov}(y, \hat{y})}{\sqrt{ \hat{\sigma}_y^2 \hat{\sigma}_{\hat{y}}^2 }} \\
  &= \frac{\sum_{i=1}^n (y_i - \bar{y})(\hat{y}_i - \bar{y})} {\sqrt{\sum_{i=1}^n (y_i - \bar{y})^2 \sum_{i=1}^n (\hat{y}_i - \bar{y})^2}} \\
  &= \frac{(y-H_1y)^T(\hat{y}-H_1y)}{\sqrt{(y-H_1y)^T(y-H_1y)(\hat{y}-H_1y)^T(\hat{y}-H_1y)}} \\
  &= \frac{y^T(H-H_1)y}{\sqrt{y^T(I-H_1)^T(I-H_1)y y^T(H-H_1)^T(H-H_1)y}} \tag{using the previous problem} \\
  &= \frac{\sqrt{y^T(H-H_1)y}}{\sqrt{y^T(I-H_1)^T(I-H_1)y}} \\
  &= \frac{\sqrt{y^T(H-H_1)y}}{\sqrt{y^T(I-H_1)y}} \tag{using a version of fact \eqref{eq:idempotent2}} \\
  &= \sqrt{\frac{SS_{reg}}{SST}} \\
  &= \sqrt{R^2} \tag{by the definition of $R^2$}
\end{align*} 
And therefore, $R^2 = \widehat{\Corr}(y, \hat{y})^2$. 


\rm

## Problem 2(c)

Show that $\hat e$ and $\hat y$ have sample correlation 0, and hence a scatter plot of $\hat e$ vs $\hat y$ should show no increasing or decreasing overall trend.

\em

The sample correlation is 0 iff the sample covariance is 0, so let's just concern
ourselves with the sample covariance and not worry about the denominator in the
correlation.

Recall that the residuals sum to 0 by construction, which also means that their sample mean is 0. Hence:
\begin{align*}
  \hat{\Cov}(\hat{e}, \hat{y}) &= \frac{1}{n}\sum_{i=1}^n (\hat{e_i} - 0)(\hat{y_i} - \bar{y}) \\
  &= \frac{1}{n}(\hat{e})^T(\hat{y} - \bar{y}) \\
  &= \frac{1}{n}(y - \hat{y})^T(\hat{y} - \bar{y}) \\
  &= \frac{1}{n}[(I-H)y]^T(H-H_1)y \\
  &= \frac{1}{n}y^T(I-H)(H-H_1)y \tag{since $(I-H)^T = (I-H)$} \\
  &= \frac{1}{n}y^T(H-H_1 - HH + HH_1)y \\
  &= \frac{1}{n}y^T(H - H_1 - H + H_1)y \tag{since $HH_1 = H_1$} \\
  &= 0
\end{align*}



\rm


\textrm{\textbf{Aside:} It is interesting to note that the results of problem \#1 (c) and problem \#2 (b) and (c) did not depend at all on the normality assumption, or even on whether the model fits the data well or not.  In other words, even if the fit is terrible, it will still be true that the sum of $\hat e$ is zero, $R^2 = [Corr(y,\hat y)]^2$, and the sample correlation between $\hat e$ and $\hat y$ is zero.  The only thing that matters is that the model have an intercept, i.e. the $X$ matrix should have a column of 1's.}

## Problem 3

[Based on Gelman \& Hill. Ch 3, \#1, p. 49] The file \texttt{ pyth.dat}, in the same folder as this hw,  contains outcome \texttt{ y} and inputs \texttt{ x1}, \texttt{ x2} for 40 data points, with a further 20 points with the inputs but no observed outcome (for this problem we will ignore these last 20 points). Save the file to your working directory and read it into R using the \texttt{ read.table()} function.

## Problem 3(a)

Fit the two models
\begin{eqnarray*}  
	\textbf{ M1:~} \texttt{ y} & = & \beta_0 + \beta_1 \texttt{ x1} + \varepsilon \\
	\textbf{ M2:~} \texttt{ y} & = & \beta_0 + \beta_1 \texttt{ x2} + \varepsilon 
\end{eqnarray*}
Which model provides a better fit for \texttt{ y}?  Why?

\em

```{r}
gh.data <- read.table("pyth.dat",header=T)
gh.data <- gh.data[apply(gh.data,1,function(x) {!any(is.na(x))}),]
str(gh.data)

M1 <- lm(y ~ x1, data=gh.data)
M2 <- lm(y ~ x2, data=gh.data)

summary(M1)
summary(M2)
##
par(mfrow=c(2,2))
plot(M1)
plot(M2)
```

The $R^2$ is much lower for M1 (0.1535) than for M2 (0.8787).  Neither set of residual diagnostic plots looks great: the residual vs fitted and scale-location plots somewhat favor M1, and the QQ plots somewhat favor M2. The Cook's distances are a bit better for M1 also.   For $R^2$ and normality of residuals, I prefer M2.


\rm

## Problem 3(b)

Construct new variables \texttt{ y2 = y$^2$}, \texttt{ x12 = x1$^2$}, and \texttt{ x22 = x2$^2$} and fit the models
\begin{eqnarray*}  
	\textbf{ M3:~} \texttt{ y2} & = & \beta_0 + \beta_1 \texttt{ x12} + \varepsilon \\
	\textbf{ M4:~} \texttt{ y2} & = & \beta_0 + \beta_1 \texttt{ x22} + \varepsilon 
\end{eqnarray*}
Compare the fits of these two models to the models in part (a).  Which fits best?  Why?

\em

```{r}
attach(gh.data)
y2 <- y^2
x12 <- x1^2
x22 <- x2^2
detach()
M3 <- lm(y2 ~ x12,data=gh.data)
M4 <- lm(y2 ~ x22,data=gh.data)

summary(M3)
summary(M4)

par(mfrow=c(2,2))
plot(M3)
plot(M4)
```

Model M4 has the highest $R^2$ (0.9397), and has residuals vs fitted and scale-location plots that are at least as good as any of the others; on the other hand, we seem to be losing normality of the residuals.  Nevertheless I prefer M4 so far.


\rm

## problem 3(c)

Fit both of the models
\begin{eqnarray*}
	\textbf{ M5:~} \texttt{ ~y} & = & \beta_0 + \beta_1 \texttt{ x1} + \beta_2 \texttt{ x2} + \varepsilon \\
	\textbf{ M6:~} \texttt{ y2} & = & \beta_0 + \beta_1 \texttt{ x12}  + \beta_2 \texttt{ x22} + \varepsilon 
\end{eqnarray*}
Compare these to the earlier models.  Which fits best?  Why?

\em


```{r}
M5 <- lm(y ~ x1 + x2, data=gh.data)
M6 <- lm(y2 ~ x12 + x22, data=gh.data)

summary(M5)
summary(M6)

par(mfrow=c(2,2))
plot(M5)
plot(M6)
```

Putting both x1 and x2 in the model for y really improved the model: M5 has an $R^2$ of 0.9724, and both predictors are significant (have coefficient estimates significantly different from zero).  However the residual diagnostic plots don't look great; in particular it seems like the residuals have a lot of right skew, and leverage seems tio increase with the size of the standardized resduals.

M6 is really winning, though: $R^2=1$, and the QQ plot shows good agreement between the residuals and the normal distribution.  There seems to be some evidence for non-constant variance though: the residuals vs fitted plot fans out as the fitted valuess increase, and the scale-location plot tells a similar story.  On the other hand, only one data point seems to have a concerning Cook's distance.

Based on all of this I like M6 best.  Looking at the estimated coefficients for M6, I notice something interesting: $\hat\beta_0$ is indistinguishable from 0, and both $\hat\beta_1$ and $\hat\beta_2$ equal 1, to at least two decimal places (even if we compute the 95\% CI's!).


\rm


## Problem 3(d)

Can you find a simple, recognizable function \texttt{ x3 = (something involving both x1 and x2)}, so that 
\[
		\textbf{ M7:~} \texttt{ y} ~ = ~ \beta_0 + \beta_1 \texttt{ x3} + \varepsilon
\]
provides a fit comparable to the best fitting models above?  What is going on?


\em

In problem 3(c) we saw that the model M6 was very nearly
\[
    y^2 = x1^2 + x2^2 + \varepsilon
\]
If we ignore $\varepsilon$, take square roots, and put $\varepsilon$ and some ``unknown'' regression coefficients back in, we get a model like
\[
    y = \beta_0 + \beta_1\sqrt{x1^2 + x2^2} + \varepsilon
\]
i.e., $y$ is the distance to the origin from some points $(x1,x2)$ in Cartesian space.

Let's try fitting this model:
```{r}
M7 <- lm(y ~ I(sqrt(x1^2 + x2^2)),data=gh.data)
summary(M7)
par(mfrow=c(2,2))
plot(M7)
```

This seems to confirm our suspicions!  $R^2=1$, the estimated regression coefficients are essentially $\hat\beta_0 = 0$ and $\hat\beta_1=1$, and the diagnostic plots looks great: 
\begin{itemize}
\item The residual vs fitted plot shows little vertical structure.
\item The QQ plot shows good adherence to normality.
\item The scale-location plot is consistent with constant-variance residuals.
\item None of the data points has Cook's distance above 0.5.
\end{itemize}


\rm
