---
title: "Homework 04 Solutions"
date: "9/28/2020"
output: 
  pdf_document:
    keep_tex: true
    fig_caption: true
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \DeclareMathOperator{\Cov}{Cov}
  - \DeclareMathOperator{\Corr}{Corr}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, include = TRUE, message = FALSE)
# library(knitr)
library(GGally)
library(reshape2)
library(tidyverse)
```

# 1.
## (1a)
We will utilize 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*} \qed

## (1b)
### (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$:
\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 \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$ by the vector $y$. In other words, $\hat{y}^*_i = \frac{1}{n}\sum_{i=1}^n y_i = \bar{y}$ for all $i$. \qed

### (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}) &= \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*} \qed

### (iii)
Recall that the sample correlation is the standardized sample covariance, defined as
\begin{align*}
  \hat{\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}} \\
  &= R^2
\end{align*} \qed

## (1c)
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*}  \qed

# 2.
## (2a)
```{r 2a}
beauty <- read.csv("../ProfEvaltnsBeautyPublic.csv")
model1 <- lm(courseevaluation ~ btystdave + female, data = beauty)
model2 <- lm(courseevaluation ~ btystdave*female, data = beauty)
coef1 <- coef(model1)  # Get model parameters
coef2 <- coef(model2)  # Get model parameters
```

### (i)
```{r 2a_i, fig.height=4, fig.cap="Problem 2(a)i: Course evaluations against beauty ratings, with no interaction (left) and with an interaction (right). \\label{fig:scatter_plots}"}
par(mfrow = c(1, 2), oma = c(2, 0, 1, 0), mar = c(4, 4, 1, 2))
# First model
plot(courseevaluation ~ btystdave, data = beauty)
abline(a = coef1["(Intercept)"], b = coef1["btystdave"], col = "red")  # line for males
abline(a = coef1["(Intercept)"] + coef1["female"], 
       b = coef1["btystdave"], col = "blue")  # line for females

# Second model
plot(courseevaluation ~ btystdave, data = beauty)
abline(a = coef2["(Intercept)"], b = coef2["btystdave"], col = "red")  # line for males
abline(a = coef2["(Intercept)"] + coef2["female"], 
       b = coef2["btystdave"] + coef2["btystdave:female"], col = "blue")  # line for females

legend(0.5, 2.6, legend = c("Males", "Females"), 
       col = c("red", "blue"), lty = c(1, 1), cex = 0.6)
```

The red lines represent the predicted values for males, and the blue lines represent predicted values for females. In the first plot, the intercept for the
red line is the `(Intercept)` coefficient, while the intercept for the blue line is the `(Intercept)` coefficient plus the `female` coefficient. The slope for both
lines is the `btystdave` coefficient.

In the second plot, the intercept for the red line is the `(Intercept)` coefficient. The intercept for the blue line again is the `(Intercept)` coefficient
plus the `female` coefficient. The slope for the red line is the `btystdave`
coefficient, while the slope for the blue line is the `btystdave` coefficient
plus the `btystdave:female` coefficient.

### (ii)
```{r 2a_ii_1, fig.cap="Problem 2(a)ii: Diagnostic plots for the first model\\label{fig:diagnostic1}", fig.height=4}
par(mfrow=c(2, 2), mar = c(4, 4.5, 2, 1))
plot(model1)
```

In Figure \ref{fig:diagnostic1}, there don't appear to be any worrying trends
between the fitted values and the residuals. However, the Q-Q plot shows some
non-normality in the residuals, on the right side of the plot; and the bottom
right plot appears to show that the residuals are left skewed, as there are a
good number of standardized residuals with values below -2. This is apparent also looking back at the scatter plots in Figure \ref{fig:scatter_plots}.

```{r 2a_ii_2, fig.cap="Problem 2(a)ii: Diagnostic plots for the second model\\label{fig:diagnostic2}", fig.height=4}
par(mfrow=c(2, 2), mar = c(4, 4.5, 2, 1))
plot(model2)
```

Figure \ref{fig:diagnostic2} for the second model, with the interaction effect, shows essentially the same results as for the first model. It's possible that the fit could be improved by transforming the predictor(s) and/or the outcome.

### (iii)
#### Here's the summary for the first model with no interaction:
```{r 2a_iii_model1}
summary(model1)
```

#### Here's the summary for the second model, which has the interaction:
```{r 2a_iii_model2}
summary(model2)
```

#### Here's a partial F-test, computed with `anova(model1, model2)`, where `model1` is the model with no interaction:
```{r 2a_iii_f-test}
anova(model1, model2)
```

The coefficient estimates for `btystdave` and `female` are quite large in both models with respect to their estimated standard errors, which is reflected in the fact that their associated p-values are quite small: the three stars (***) mark that they are less than 0.001. Notice that the estimates differ between the two models. This is not surprising; in general, adding a new variable to a model will change the estimates of all the coefficients, unless that new variable is completely uncorrelated with the previous variables. Since the "new" variable here is an interaction between the other two variables, it is naturally correlated with those variables, so we should expect the estimates to change.

The adjusted $R^2$ values for both models are quite small, indicating that the predictors do not account for much of the total variance. The addition of the interaction term did not increase the adjusted $R^2$ by a non-trivial amount.

The estimate of `btystdave:female`, the interaction term, is not significantly different from 0 at the conventional $\alpha = 0.05$ level. The actual magnitude of the estimate is also smaller than the main effect estimates. The partial F-test here constitutes a test of the null hypothesis
\begin{align*}
  H_0: \beta_3 = 0
\end{align*}
where $\beta_3$ is the coefficient associated with `btystdave:female`. Hence, the p-value of the partial F-test that is printed in the `anova` output is equivalent to the p-value of the t-test in the output of `summary(model2)`, namely, $p = 0.0789$. 

Since we don't have enough evidence to reject $H_0$, and since the addition of the interaction term term did not increase the adjusted $R^2$, in the interest of parsimony, we may decide based on these results to exclude the interaction term.

## (2b)
We don't have a lot of formal tools available for variable selection at this stage, so the general approach here is to think in real-world terms about what variables seem likely to be related to the outcome and see what happens when they are included in the model. There's no one correct approach or outcome here. In addition to `btystdave` and `female`, let's consider the variables `age`, `minority`, `nonenglish`, `tenured`, and `onecredit`. First, let's look at a pairs plot with `courseevaluation`:
```{r 2b, fig.cap="Pairs plot for candidate predictors in 2(b)\\label{fig:2b_pairs}"}
cols <- c("minority", "nonenglish", "tenured", "onecredit")
beauty[, cols] <- lapply(beauty[, cols], as.factor)
ggpairs(beauty[, c("courseevaluation", "age", cols)])
```

The marginal distributions of each variable appear on the diagonal, while the top row and leftmost column give different visualizations of the relationship between the course evaluation scores and the other five variables. Among the candidate predictors, only `onecredit` shows a clear marginal relationship with `courseevaluation`, but since these are marginal plots, that doesn't mean the other variables won't be informative to some degree. There are very few one credit courses represented here (27 out of 463 rows), but nevertheless, let's try including it, along with the following sets of variables:
\begin{align*}
  &\text{Model m1: onecredit} \\
  &\text{Model m2: onecredit, minority, nonenglish} \\
  &\text{Model m3: onecredit, minority, nonenglish, tenured, age} 
\end{align*}
For each model, we'll look at a model summary and diagnostic plots. We can also
run an F-test comparing the model from part 2(a)ii (aka `model1`) to the expanded model.

### Model m1
`m1 <- lm(courseevaluation ~ btystdave + female + onecredit)`
```{r 2b_models1, fig.height=4, fig.cap="Diagnostic plots for model m1\\label{fig:2b_m1}"}
m1 <- lm(courseevaluation ~ btystdave + female + onecredit, data = beauty)
summary(m1)
par(mfrow=c(2, 2), mar = c(4, 4.5, 2, 1))
plot(m1)
```

#### F-Test for model m1 vs. model1
```{r 2b_model1_ftest}
anova(model1, m1)
```

### Model m2
`m2 <- lm(courseevaluation ~ btystdave + female + onecredit + minority + nonenglish)`
```{r 2b_models2, fig.height=4, fig.cap="Diagnostic plots for model m2\\label{fig:2b_m2}"}
m2 <- lm(courseevaluation ~ btystdave + female + onecredit + minority + nonenglish, 
         data = beauty)
summary(m2)
par(mfrow=c(2, 2), mar = c(4, 4.5, 2, 1))
plot(m2)
```

#### F-Test for model m2 vs. model1
```{r 2b_model2_ftest}
anova(model1, m2)
```

### Model m3
`m3 <- lm(courseevaluation ~ btystdave + female + onecredit + minority + nonenglish)`
```{r 2b_models3, fig.height=4, fig.cap="Diagnostic plots for model m3\\label{fig:2b_m3}"}
m3 <- lm(courseevaluation ~ btystdave + female + onecredit + minority + nonenglish + tenured + age,
         data = beauty)
summary(m3)
par(mfrow=c(2, 2), mar = c(4, 4.5, 2, 1))
plot(m3)
```

#### F-Test for model m3 vs. model1
```{r 2b_model3_ftest}
anova(model1, m3)
```

The diagnostics for all three models look reasonable, and they look fairly similar to the diagnostics for `model1`, the model that we retained from part 2(a)iii. The F-tests comparing each model to the baseline `model1` are all highly significant. Summaries for models `m1` and `m2` each model show highly significant coefficients. In model `m3`, however, the tests for `tenured` and `age` are non-significant, and the coefficient estimates are tiny. Given these tiny estimates, it would be reasonable to drop these variables from the model and retain model `m2`.

If we believe this model, then we can say that course evaluations bear at least a roughly linear relationship to ratings of beauty, whether a person is female or male, whether a course is one credit or not, whether the professor is a minority, and whether the professor is a native English speaker or not, with the sign and magnitude of the estimated coefficients indicating the nature and strength of the relationships. It is tempting to draw causal conclusions here, but we can't really do that without checking some additional assumptions which are beyond the scope of the problem. However, these patterns may at least be suggestive of causal relationships, which we can subsequently go about testing more rigorously.

# 3 (skipped)

