---
title: "Homework 04 Solutions"
date: "9/27/2021"
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 = TRUE, include = TRUE, message = FALSE)
# library(knitr)
library(GGally)
library(reshape2)
library(tidyverse)
```

# 1.
## (1a)
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*} \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_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$. \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}) &= \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*} \qed

### (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$. \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. (Sheather, Ch 5, pp. 146–147, \#1, with extra questions.)

## (2a) Develop model to predict LATE from BILL.

Let's start with a little EDA\ldots
```{r 2a scatter i,fig.height=3,fig.width=4,fig.align='center'}
overdue <- read.table("overdue.txt",header=T)
str(overdue)
plot(LATE ~ BILL, data=overdue)
```

From the scatter plot, it looks like regression is going to fail, because LATE is not a single-valued function of BILL.  Indeed, if we do a simple regression of LATE on BILL, we get
```{r 2a simple reg}
round(summary(lm(LATE ~ BILL,data=overdue))$coefficients,2)
```
and the coefficient on BILL is clearly not significantly different from zero.  

Looking back at the problem statement, we also know that the first 48 observations are for residential accounts, and the last 48 are for commercial accounts.  Maybe this helps explain the structure we see above?   

So, let's create a dummy variable for "residential", and colored the points according to this dummy variable.
```{r 2a scatter ii,fig.height=3,fig.width=4,fig.align='center'}
overdue$residential <- c(rep(1,48),rep(0,48))
plot(LATE ~ BILL, data=overdue,col=residential + 1)
legend(250,30,legend=c("Residential","Commercial"),pch=1,col=c(2,1),cex=0.5)
```
Much better: Now we clearly see that the two "branches" in the scatter plot correspond to the two types of account, with different slopes and intercepts for each type of account.  This suggests we should fit an ANCOVA model with interactions:
```{r 2a scatter iii,fig.height=3,fig.width=4,fig.align='center'}
lm.2a <- lm(LATE ~ BILL*residential, data=overdue)
print(coefs <- coef(lm.2a))
plot(LATE ~ BILL, data=overdue,col=residential + 1)
abline(a=coefs[1], b=coefs[2],col=1)
abline(a=coefs[1]+coefs[3], b=coefs[2]+coefs[4],col=2)
legend(250,30,legend=c("Residential","Commercial"),pch=1,col=c(2,1),cex=0.5)
```
The new scatter plot suggests that this may be pretty good model for the data.


## (2b) Summary \& Residual Diagnostics.

Here is the model summary and the set of four standard diagnostic plots:
```{r 2b summary}
summary(lm.2a)
par(mfrow=c(2,2))
plot(lm.2a)
```
This seems like a pretty good model.  The regression coefficients are all strongly significantly different from zero, and the $R^2=0.98$.   Looking at the residual diagnostic plots, we see that the residuals are centered at zero and show no vertical patterns; the QQ plot suggests that the residuals are approximately normal with just 3 or 4 outliers, the Scale-Location plot is consistent with constant variance, and there aer no points with Cook's distance above 0.5.  The points with high leverage are generally not outliers, and the points that are outliers have modest leverage.

## (2c) "Under 60 days or your money back?"

To think about the "Under 60 days or your money back!" slogan, let's look at the final scatter plot with the fitted regression model again, but with horizontal lines at 60, 90 and 100 days overlaid:
```{r 2c scatter iv}
plot(LATE ~ BILL, data=overdue,col=residential + 1)
abline(a=coefs[1], b=coefs[2],col=1)
abline(a=coefs[1]+coefs[3], b=coefs[2]+coefs[4],col=2)
abline(h=60,col="Blue")
abline(h=90,col="Yellow")
abline(h=100,col="Green")
legend(250,30,legend=c("Residential","Commercial"),pch=1,col=c(2,1))
```
It looks like 60 days is enough for residential accounts (even if we extrapolate a little to the right of the data using the regression line), but we need up to 100 days for commercial accounts.

If we want the slogan to work for almost all collections, perhaps we should recommend
\begin{itemize}
\item "Under 90 days or your money back!"   $~~~$or perhaps
\item "Under 100 days or your money back!"
\end{itemize}
since 90 days is enough for most collections, even commercial ones, and 100 days is enough for viturally all of them.

On the other hand, if we can restrict the advertising to Residential accounts only, then "Under 60 days" seems just fine: Based on this plot, we would expect virtually all Residential account to be collected in 60 days or less.

# 3. (Sheather, Ch 5, p. 147, \#2.)

We begin with a little exploration, and finding a good model.
```{r 3 initial eda and model}
houston <- read.csv("HoustonChronicle.csv",header=T)
str(houston)

cat("\n")
for (i in names(houston)) {
  cat(i,"\n")
  print(summary(houston[,i]))
  cat("Number of unique values:",length(unique(houston[,i])),"\n\n")
}
```
Apparently, there are only two "year" values: 1994 and 2004; so we will convert this to a factor variable, just to make interpretation of the coefficients easier; and we make a scatterplot matrix ("pairs" plot) of all of the variables except district.  I'm also going to rename the two continuous variables to something more suggestive, and re-order the variables, so that the "pairs" plot makes a little more sense.
```{r 3 ggpairs}
houston$Year <- as.factor(houston$Year)
names(houston)[c(2,3)] <- c("Pct.Repeating.1st.Grade", "Pct.Low.income.students")
houston <- houston[,c(1,3,2,4,5)]
ggpairs(houston[,-1])  # there are too many districts for the plot, so we omit it
```
There is some evidence in the "pairs" plots for all three of the hypotheses in parts (a), (b) and (c), but we now fit regression and ANCOVA models to see which effects are "significant" (i.e. large enough that they are not likely just due to noise in the data).

## (3a) Low income associated with repeating first grade?

Let's start with a simple regression of repeating first grade on low income
```{r 3a model,fig.height=3,fig.width=4,fig.align='center'}
summary(lm.3a <- lm(Pct.Repeating.1st.Grade ~ Pct.Low.income.students,data=houston))
par(mfrow=c(1,1))
plot(Pct.Repeating.1st.Grade ~ Pct.Low.income.students,data=houston)
abline(lm.3a)
```
```{r 3a model diagnostics}
par(mfrow=c(2,2))
plot(lm.3a)
```
This actually looks like a good regression; the summary shows a very significant increasing association between percent of students repeating first grade and percent of students with low (family) incomes: the coefficient on percent low income students is $\hat\beta_1 = 0.07550$, with $SE(\hat\beta_1)=0.01823$.  The effect seems rather small, however---we expect an increase of only 0.08\% of kids repeating first grade, for every 1\% increase in kids in poverty.

## (3b) More students repeating first grade in 2004-2005 than in 1994-1995?

For this question, we will just regress repeating first grade on year
```{r 3b model,fig.height=3,fig.width=4,fig.align='center'}
summary(lm.3b <- lm(Pct.Repeating.1st.Grade ~ Year,data=houston))
par(mfrow=c(1,1))
jy <- jitter(as.numeric(houston$Year)-1)
plot(Pct.Repeating.1st.Grade ~ jy,data=houston,xlab="Year (jittered)")
abline(lm.3b)
```
```{r 3b model diagnostics}
par(mfrow=c(2,2))
plot(lm.3b)
```
Again, this looks like a good regression model, in the sense that the assumptions underlying linear regression are approximately satisfied\footnote{Remember, the clustering around different $x$ values is just due to the fact that $x$ is categorical.  The important thing here is that we don't see any interesting patterns within or between those clusters.}.  However, the estimated coefficient on Year, $\hat\beta_1=1.2$ with $SE(\hat\beta_1)=0.73$ is not statistically significantly different from zero.  From this we would conclude that there really is not enough evidence to say that there are more students repeating first grade in the 2004 school year than in 1994 school year.


## (3c) Difference in the relationship of income with repeating between the two school years?

Now let's fit the interactive ANCOVA model.  This will both help answer question (3c), and help us to understand whether our answers to questions (3a) and (3b) need any additional elaboration.
```{r 3c model fit}
summary(lm.3c <- lm(Pct.Repeating.1st.Grade ~ Pct.Low.income.students*Year,data=houston))
```
The summary shows that neither the main effect for Year, nor the interaction between low income and Year, are significant in the model.  This is illustrated in the scatterplot below: the regression lines for the two categories (year=1994-1995 and year=2004-2005) are nearly identical.  The residual diagnostic plots that follow suggest that the modeling assumptions hold up fairly well here.
```{r 3c model graph,fig.height=3,fig.width=4,fig.align='center'}
par(mfrow=c(1,1))
plot(Pct.Repeating.1st.Grade ~ Pct.Low.income.students,data=houston,
     col=as.numeric(houston$Year))
coefs <- coefficients(lm.3c)
abline(a=coefs[1],b=coefs[2],col=1)
abline(a=coefs[1]+coefs[3],b=coefs[2]+coefs[4],col=2)
legend(78,18,legend=c("1994-1995","2004-2005"),pch=1,col=c(1,2),cex=0.5)
```
```{r 3c model diagnostics}
par(mfrow=c(2,2))
plot(lm.3c)
```
We can also apply an F test using the ANOVA function in R to compare the full ANCOVA model in part (c) with the models in parts (a) and (b).
```{r}
anova(lm.3a,lm.3c)
```
Since the $F$ statistic testing $H_0$: ${\tt Pct.Repeating.1st.Grade ~\sim~ Pct.Low.income.students}$ vs $H_A$: ${\tt Pct.Repeating.1st.Grade ~\sim~ Pct.Low.income.students * Year}$ is nonsignificant ($F=0.2541$, $p=0.7761$), we cannot reject the simpler model in part (a): If we are starting with  ${\tt Pct.Repeating.1st.Grade ~\sim~ }$ ${\tt Pct.Low.income.students}$, there is really no need to add ${\tt Year}$ (let alone an interaction with ${\tt Year}$) to the model.

```{r}
anova(lm.3b,lm.3c)

```
Here the $F$ statistic testing $H_0$: ${\tt Pct.Repeating.1st.Grade ~\sim~ Year}$ vs $H_A$: ${\tt Pct.Repeating.1st.Grade ~\sim~ }$ ${\tt Pct.Low.income.students * Year}$ is highly significant ($F=7.2221$, $p=0.001099$), so we would reject the simpler model in part (b): starting with ${\tt Pct.Repeating.1st.Grade ~\sim~ Year}$, we really do get a better model by adding ${\tt Pct.Low.income.students}$ to the model.

Note that we cannot perform an $F$ test like ${\tt anova(lm.3a, lm.3b)}$ directly\footnote{There are ways to compare non-nested models, e.g. cross-validation, AIC, BIC, etc., but it doesn't work with $F$ tests or likelihood ratio tests.}, since these models are not nested.  We can infer, informally, from the above tests though, that of all three models,  ${\tt Pct.Repeating.1st.Grade ~\sim~ Pct.Low.income.students}$ seems to be the best.  Thus, there is no difference between years in the way that low income is related to repeating first grade.

# 4. (course evaluations question from Gelman  \& Hill)
## (4a)
```{r 4a}
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 4a_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 4a_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 4a_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 4a_iii_model1}
summary(model1)
```

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

#### Here's a partial F-test, computed with `anova(model1, model2)`, where `model1` is the model with no interaction:
```{r 4a_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.

## (4b)
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 4b, 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 4b_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 4b_model1_ftest}
anova(model1, m1)
```

### Model m2
`m2 <- lm(courseevaluation ~ btystdave + female + onecredit + minority + nonenglish)`
```{r 4b_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 4b_model2_ftest}
anova(model1, m2)
```

### Model m3
`m3 <- lm(courseevaluation ~ btystdave + female + onecredit + minority + nonenglish)`
```{r 4b_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 4b_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.

## (end of solutions)
