---
title: "Homework 04 Solutions"
date: "2022-09-22"
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...
library(GGally) ## for ggpairs...

library(leaps) ## regsubsets(), summary(), coef()
library(car)   ## subsets(), mmps(), vif(), etc.

```


# Problem 1: Sheather, Ch 6, pp. 216--221, \#3.  

Continuing the analysis of data in "cars04.csv"\ldots{}
Notes:
\begin{itemize}
\item The data set is available as ``cars04.csv''.
\item There are some errors in the results for fitting the model (6.37) in the textbook, so I suggest you refit the models for this problem directly from the data in ``cars04.csv''.
\end{itemize}

The variables are:
\begin{description}
\item[$Y =$] Suggested Retail Price; 
\item[$x_1 =$] Engine size; 
\item[$x_2 =$] Cylinders;
\item[$x_3 =$] Horse power; 
\item[$x_4 =$] Highway mpg; 
\item[$x 5 =$] Weight;
\item[$x_6 =$] Wheel Base;
\item[$x_7 =$] Hybrid (1=Hybrid; 0=just gasoline-powered)
\end{description}

The first model proposed is
$$
Y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_3 + \beta_4x_4 + \beta_5x_5 + \beta_6x_6 + \beta_7x_7 + \epsilon
\eqno{(6.36)}
$$
After examining Box-Cox transformation proposals and "rounding" to more interpretable values, the second model proposed is
$$
\log(Y) = \beta_0 + \beta_1x_1^{0.25} + \beta_2
\log(x_2) + \beta_3\log(x_3) + \beta_4(1/x_4) + \beta_5x_5 + \beta_6\log(x_6) + \beta_7x_7 + \epsilon
\eqno{(6.37)}
$$
\em

Since there are some problems with the analysis presented in Sheather, we refit all the models and reproduce all the outputs (except that I am just going to accept the transformations suggested in (6.37)).

\rm

## Problem 1(a)

Decide whether (6.36) is a valid model. Give reasons to support your answer.

\em

We start by refitting the model and reproducing the diagnostic plots in Sheather.
```{r}
cars <- read.csv("cars04.csv")
str(cars)
names(cars)
cars.red <- cars[,c("EngineSize","Cylinders",
"Horsepower","HighwayMPG","Weight","WheelBase","Hybrid","SuggestedRetailPrice")]
dim(cars.red)
names(cars.red)
ggpairs(cars.red)
  ## Sheather did not include the response variable "SuggestedRetailPrice",
  ## but I think it is good to include since it can help us guess transformations...
summary(lm.6.36 <- lm(SuggestedRetailPrice ~ ., data=cars.red))
par(mfrow=c(2,2))
plot(lm.6.36)
```

Neither the scatter plot matrix nor the \verb|summary()| printout can tell us whether the model is valid (the scatter plot matrix can help us make guesses about transformations and/or collinearity.  The \verb|summary| output can tell us whether the fitted model will be useful for answering questions our collaborator may have).  Validity of the model depends only on
\begin{itemize}
\item Linear relationship between $Y$ and the $X$'s
\item Normal errors $\epsilon$
\item Constant error variance $\sigma^2=Var(\epsilon)$
\end{itemize}
We can use the casewise diagnostic plots for this:
\begin{description}
\item[Linear relationship between $Y$ and the $X$'s:] The ``Residuals vs. Fitted'' plot suggests a curved relationship between the residuals $\hat e$ and the fitted values $\hat y$.  This suggests that the relationship between $Y$ and the $X$'s is \ul{nonlinear}. 
\item[Normal errors $\epsilon$:] The ``Normal Q-Q'' plot of the standardized residuals $r_i$ shows quite a bit of right-skew.  This suggests that the errors $\epsilon$ are \ul{not normally distributed}.
\item[Constant error variance $\sigma=\sqrt{Var(\epsilon)}$:] The red trend line in the ``Scale-Location'' is increasing from left to right, suggesting that the variance $Var(\epsilon)$ is not constant across the data.
\end{description}
Since there is evidence of violations of all three basic assumptions of the model, we can say \ul{the model is not fully valid}.

\textbf{Note:}  The model might still be useful in some way, e.g. for prediction.  However it does not satisfy the basic statistical assumptions above, and so any tests or estimates for the $\beta$'s, for the $\hat y$'s for $\sigma$ should be viewed with some suspicion.

\rm

## Problem 1(b)

The plot of residuals against fitted values produces a curved pattern. Describe what, if anything can be learned about model (6.36) from this plot.

\em

The curved pattern in the residuals suggests we should look for transofmrations of $X$'s or $Y$'s to produce a set of raw residuals $\hat e$ with less functional dependence on the fitted values $\hat y$.  

\rm

## Problem 1(c)

Identify any bad leverage points for model (6.36).

\em

From the "Residuals vs Leverage" plot we can see: 
\begin{itemize} 
\item Observation \#223 seems to be influential, with a Cook's Distance $>$ 0.5.
\item Observation \#67 may also be influention, with a Cook's Distance near 0.5.
\item Observation \#222 does not have high leverage, but its residual is a hugh outlier---more than 6 SD's from the middle of the residual distribution.
\item Two more observations, not labelled by R, have relatively high leverage ($\approx$ 0.35) but pretty small residuals.
\end{itemize}
\textbf{Note:} The usual cutoff for "high leverage" is 
\[
        h_{high} = 2\cdot\frac{p+1}{n} = 2\cdot\frac{8}{234} \approx `r round(2*8/234,2)`
\]
so there are actually 5 or more points with leverage $h_{ii}>h_{high}$, but the points I've identified above are the ones that I might devote some worry to in practice.


\rm

## Problem 1(d)

Decide whether (6.37) is a valid model.

\em

First we refit the model and reproduce the output in Sheather (this is the part that Sheather appaers to have gotten wrong):
```{r}
attach(cars.red) ## so I don't have to type "cars.red$" all the time...
tSuggestedRetailPrice <- log(SuggestedRetailPrice)
tEngineSize <- EngineSize^(0.25)
tCylinders <- log(Cylinders)
tHorsepower <- log(Horsepower)
tHighwayMPG <- 1/HighwayMPG
tWheelBase <- log(WheelBase)
cars.t <- data.frame(tEngineSize,
                       tCylinders, tHorsepower, tHighwayMPG,
                       Weight, tWheelBase, Hybrid,tSuggestedRetailPrice)
detach()
names(cars.t)
ggpairs(cars.t)
summary(lm.6.37 <- lm(tSuggestedRetailPrice ~ ., data=cars.t))
par(mfrow=c(2,2))
plot(lm.6.37)
mmps(lm.6.37)

```

Now, following the same approach as before, we use the casewise diagnostic plots to check:
We can use the casewise diagnostic plots for this:
\begin{description}
\item[Linear relationship between $Y$ and the $X$'s:]  The red trend-line in the ``Residuals vs Fitted'' plot is nearly flat, suggesting no relationship between $\hat e$ and $\hat y$.  This suggests that \ul{we have an approximately linear} \ul{relationship} between the new $Y$ and the new $X$'s in the transformed model.
\item[Normal errors $\epsilon$:] Except for two outliers (obs. \#222 and \#67) the ``Normal Q-Q'' plot confirms that the standardized residuals $r_i$ are approximately normally distributed, suggesting that the errors $\epsilon$ are \ul{normally distributed}.
\item[Constant error variance $\sigma=\sqrt{Var(\epsilon)}$:] The red trend line in the ``Scale-Location'' plot is nearly flat, suggesting that the error variance $Var(\epsilon)$ is constant across the data set. 
\end{description}
Of course no fit will be perfectly valid, but the evidence from the casewise diagnosic plots suggests that (6.37) is a much more valid model than (6.36).

\textbf{Note:} We can also see, from the scatterplot matrix and from the marginal model plots, that we seem to have done a good job with transforms.  Especially, in the marginal model plots, the estimates of $E[Y|X]$ from the nonparametric regressions (blue) and the linear model (red) line up well, so it doesn't appear that any further transformations are needed.

\rm

## Problem 1(e)

To obtain a final model, the analyst wants to simply remove the two insignificant predictors $1/x_4$ (i.e., \texttt{tHighwayMPG}) and $\log(x_6)$ (i.e., \texttt{tWheelBase}) from (6.37). Perform a partial $F$-test to see if this is a sensible strategy.

\em

We will fit the smaller model (without the two variables that the analyst wants to remove) and then compare that with \verb|lm.6.37|:
```{r}
lm.smaller <- update(lm.6.37, . ~ . - tHighwayMPG - tWheelBase)
anova(lm.smaller, lm.6.37)
```
With a $p$-value of 0.97, there reall is no evidence in favor of keeping these two variables in the model, and so \verb|lm.smaller| seems like a sufficient model:
```{r}
summary(lm.smaller)
```

We probably should go ahead and check casewise diagnostic plots, marginal model plots, etc., to be sure that we still have a valid model and no more transformations are needed, but I will stop here.

\rm

## Problem 1(f)

The analyst’s boss has complained about model (6.37) saying that it fails to
take account of the manufacturer of the vehicle (e.g., BMW vs Toyota).
Describe how model (6.37) could be expanded in order to estimate the effect
of manufacturer on suggested retail price.

\em

This is slightly tricky to do well.  The basic idea (and, really, all you have to say to answer the question) is that you want to extract the manufacturer name from the variable \verb|Vehicle.Name| in the orignal data set \verb|cars|, and add that to the regressions above, and see what happens...

Doing the extraction is a little tricky but here's a way to proceed, if you are curious:

The \verb|Vehicle.Name| variable in the original data set is nearly unique for each car in the data set:
```{r}
dim(cars)
length(unique(cars$Vehicle.Name))
```
(There are two cars duplicated in the data set,
```{r}
tmp <- table(cars$Vehicle.Name)
tmp[tmp>1]
cars[grep(names(tmp[tmp>1])[1],cars$Vehicle.Name),]
cars[grep(names(tmp[tmp>1])[2],cars$Vehicle.Name),]

```
but this doesn't change the problem.)

With 232 unique vehicle names for 234 observation in the data set, if we just include vehicle name in a regression, it will soak up all of the variability explained by the other variables, and we won't get any useful information out.

We note that the make (i.e. the manufacturer name) of the car always comes first in the vehicle name, followed by a space.  So we can make a new variable, \verb|Make|, that just has the manufacturers' names:
```{r}
splits <- strsplit(cars$Vehicle.Name," ")
head(splits,10)
mfr.name <- sapply(splits, function(x) x[1])
head(mfr.name,10)
cars$Make <- mfr.name
length(unique(cars$Make))
table(cars$Make)
```
Now we could introduce the \verb|Make| variable into any of the regressions we fitted above, to see if the car maker affect prediction of the \verb|SuggestedRetailPrice|, separately from the other features of the car.


\rm


# Problem 2: Sheather, Ch 7, p. 261, \#3. 

Continuing with analysis of "pgatour2006.csv"\ldots{} 
Note:
\begin{itemize}
\item The data is available in ``pgatour2006.csv''.
\end{itemize}

The variables are:
\begin{description}
\item[$Y =$] PrizeMoney; 
\item[$x_1 =$] DrivingAccuracy; 
\item[$x_2 =$] GIR;
\item[$x_3 =$] PuttingAverage; 
\item[$x_4 =$] BirdieConversion; 
\item[$x 5 =$] SandSaves; 
\item[$x_6 =$] Scrambling;
\item[$x_7 =$] PuttsPerRound
\end{description}
and the model we are considering is 
$$
\log(Y) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_3 + \beta_4x_4 + \beta_5x_5 + \beta_6x_6 + \beta_7x_7 + \epsilon
\eqno{(7.10)}
$$
We want to do variable selection to choose a subset of the predictors to model $\log(Y)$.

## Problem 2(a)

Identify the optimal model or models based on $R^2_{adj}$, AIC, AIC$_C$, BIC from
the approach based on all possible subsets.

\em

Remember that "all subsets" selection divides the possible models into groups according to the number of predictors, up to the maximum number of predictors (seven, in our case):
\begin{itemize}
\item First find the 1-predictor model with the smallest RSS.
\item Then find the 2-predictor model with the smallest RSS.
\item \ldots{} and so forth up to the 7-predictor model with the smallest RSS.
\end{itemize}
These seven models can then be compared with AIC, BIC, etc.

```{r}
golf <- read.csv("pgatour2006.csv")

golf.red <- golf[,c("PrizeMoney","DrivingAccuracy","GIR","PuttingAverage",
                    "BirdieConversion","SandSaves","Scrambling","PuttsPerRound")]

all.subsets <- regsubsets(log(PrizeMoney) ~ ., data=golf.red)

``` 

We can get some of the results we want from the "subsets()" function:

```{r}
subsets(all.subsets,statistic="adjr2",legend=FALSE)
subsets(all.subsets,statistic="bic",legend=FALSE)
```

But to completely do the problem we need another strategy, since the "subsets()" function does not know AIC or CAIC (see the "statistic" parameter in Figure 1 below):
```{r eval=FALSE}
help(subsets)  
## opens a browser window with the help, but I've included part of the 
## contents of the browser window in Figure 1 below.
```

![help(subsets)]("help(subsets).jpg"){width=100%}

So, we have to calculate AIC and CAIC (and in the process we'll also see how to do BIC)

Note that from lecture 08, slide 7, at the bottom, we can write the log-likelihood as
\[
  \mbox{(log-likelihood)} = c_1(n) - c_2(n)\log(RSS) = -\frac{n}{2}\log(2\pi) - \frac{n}{2}\log(RSS)
\]
which means we can calculate AIC by hand as  
\begin{verbatim}  
    AIC = -2(log-liklihood) + 2*(p+2) = [n*log(2pi)] + n*log(RSS) + 2*(p+2)
\end{verbatim}
and since the \verb|[n*log(2pi)]| term will cancel when we subtract AIC's we can 
ignore it and just write
\begin{verbatim}
    AIC = n*log(RSS) + 2*(p+2)
\end{verbatim}
Similarly, 
\begin{verbatim}
    CAIC = AIC + 2*(p+2)*(p+3)/(n-p-1) = n*log(RSS) + 2*(p+2) + 2*(p+2)*(p+3)/(n-p-1)
\end{verbatim}
and
\begin{verbatim}
    BIC = n*log(RSS) + log(n)*(p+2)
\end{verbatim}
We can get the RSS for each model from 

```{r}
tmp <- summary(all.subsets)  
names(tmp)

p <- 1:7 ## number of parameters in each subset model
n <- dim(golf.red)[1]

attach(tmp)
results <- data.frame(which,rss,adjr2,bic=n*log(rss)+log(n)*(p+2),aic=n*log(rss)+2*(p+2),
                      caic=n*log(rss) + 2*(p+2) + 2*(p+2)*(p+1)/(n-p-1))
detach()
```

Note that different authors will use somewhat different definitions of AIC and BIC.
The differences are usually just in what is done with the constants $c_1(n)$ and $c_2(n)$, so the value of the criterion changes, but the model that minimizes the criterion does not change.


We can print the "results" data frame and get the minumum by hand:

```{r}
results
```

Or we can write a function to identify the row that minimizes each criteron:

```{r}
minimize <- function(res,col) {
  obj <- res[,col]
  k <- (1:length(obj))[obj==min(obj)]
  return(res[k,])
}

maximize <- function (res,col) {
  obj <- res[,col]
  k <- (1:length(obj))[obj==max(obj)]
  return(res[k,])
}

maximize(results,"adjr2")

minimize(results,"bic")

minimize(results,"aic")

minimize(results,"caic")

```

Reading the "TRUE" and "FALSE" values as when to include or not include a variable in the model, we see that for "all subsets" selection:
\begin{itemize}
\item The optimal $R^2_{adj}$ model is\\ \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + SandSaves + PuttPerRound|
\item The optimal BIC model is\\ \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling|
\item The optimal AIC model is\\ \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + SandSaves + PuttsPerRound|
\item The optimal CAIC model is\\ \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + SandSaves + PuttsPerRound|
\end{itemize}
It's a bit hard to read the graph that \verb|subsets()| gives us for $R^2_{adj}$ but clearly the best BIC model we identified is the same as the one that \verb|subsets()| gives us for BIC.

\rm

## Problem 2(b) 

Identify the optimal model or models based on AIC and BIC from the
approach based on backward selection.

\em

We can specify a version of "backward selection" with the "method" parameter in the
\verb|regsubsets()| function.  Instead of looking at all models at each subset size, this 
\begin{itemize}
\item First calculates RSS for the largest model, 7 predictors, in our case 
\item Then takes the 6-predictor model that increases RSS the least from the 7-predictor model
\item Then takes the 5-predictor model that increases RSS the least from the 6-predictor model
\item \ldots{} and so forth down to the 1-predictor model
\end{itemize}
We can then compare these seven models with AIC, BIC, etc.
```{r}

backward <- regsubsets(log(PrizeMoney) ~ ., data=golf.red, method ="backward")

tmp <- summary(backward)
names(tmp)

attach(tmp)
results <- data.frame(which,bic=n*log(rss)+log(n)*(p+2),aic=n*log(rss)+2*(p+2))
detach()

results

minimize(results,"bic")

minimize(results,"aic")

```

We see that for "backward selection"
\begin{itemize}
\item The optimal BIC model is\\ \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling|
\item The optimal AIC model is\\ \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + SandSaves + PuttsPerRound|
\end{itemize}

\rm

## Problem 2(c) 

Identify the optimal model or models based on AIC and BIC from the
approach based on forward selection.

\em
We can specify a version of "forward selection" with the "method" parameter in the
\verb|regsubsets()| function.  Instead of looking at all models at each subset size, this 
\begin{itemize}
\item First finds the 1-predictor model with the smallest RSS
\item Then takes the 2-predictor model that decreases RSS the most from the 1-predictor model
\item Then takes the 3-predictor model that decreases RSS the most from the 2-predictor model
\item \ldots{} and so forth up to the 7-predictor model
\end{itemize}
We can then compare these seven models with AIC, BIC, etc.

```{r}

forward <- regsubsets(log(PrizeMoney) ~ ., data=golf.red, method ="forward")

tmp <- summary(forward)
names(tmp)

attach(tmp)
results <- data.frame(which,bic=n*log(rss)+log(n)*(p+2),aic=n*log(rss)+2*(p+2))
detach()

results

minimize(results,"bic")

minimize(results,"aic")



```

We see that for "forward selection"
\begin{itemize}
\item The optimal BIC model is\\ \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + PuttsPerRound|
\item The optimal AIC model is\\ \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + SandSaves + PuttsPerRound|
\end{itemize}


\rm

## Problem 2(d) 

Carefully explain why the models chosen in (a) & (c) are not the same while
those in (a) and (b) are the same.

\em

Here are the results that I got from model selection:
\begin{description}

\item[Best Models from 2(a) ``all subsets'':] ~
\begin{itemize}
\item BIC: \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling|
\item AIC: \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + SandSaves + PuttsPerRound|
\end{itemize}

\item[Best Models from 2(b)``backward'':] ~
\begin{itemize}
\item BIC: \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling|
\item AIC: \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + SandSaves + PuttsPerRound|
\end{itemize}

\item[Best Models from 2(c) ``forward'':] ~
\begin{itemize}
\item BIC: \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + PuttsPerRound|
\item AIC: \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + SandSaves + PuttsPerRound|
\end{itemize}

\end{description}
We know that the "all subsets" models are the best of all possible models on the seven variables we have to work with, whereas the "backward" and "forward" models come from heuristics that may or may not find the best of all possible models on these seven variables.
\begin{itemize}

\item The BIC model for "backward" is the same as for "all subsets" because the "backward" selection method stumbled on the best overall model for BIC (look at line 3 of the ``results'' data frame in parts (a) and (b) to verify this).

\item The BIC model for "forward" selection is different because the "forward" heuristic didn't identify the best overall model for BIC (look at line 3 of the ``results'' output in part (a), vs line 4 of the ``results'' output in part (c) to verify this. Note also that the ``line 3'' model in part (c) is \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + PuttsPerRound| instead of the BIC-optimal model \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling|.

\item The AIC models are the same in each case, because both the "forward" and the "backward" heuristics stumbled on the best overall model for AIC (look at line 5 of the ``results'' data frame in parts (a), (b) and (c) to verify this).
\end{itemize}

\rm

## Problem 2(e) 

Recommend a final model. Give detailed reasons to support your choice.

\em

There are really only two models to consider, from the model selection above.  They are:
```{r}
lm.AIC <- lm(log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling + 
               SandSaves + PuttsPerRound, data=golf.red)

lm.BIC <- lm(log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling,data=golf.red)
```
Looking at the casewise diagnostic plots,
```{r fig.height=2.5,fig.width=8}
par(mfrow=c(1,4))
plot(lm.AIC)
plot(lm.BIC)
```
we see that both models seem to satisfy the assumptions of linear regression equally well.
\begin{description}
\item[Linearity:] There is no real trend in the Residuals vs Fitted plot for either model, so no transforms seem to be needed for either model.
\item[Normality:] The Normal Q-Q plots for both models look about the same.
\item[Constant Variance:] Except for ``edge effects'' the trend lines in both Scale-Location plots are pretty horizontal, supporting the idea of constant variance for both models.
\end{description}
In addition the Residuals vs Leverage plots look about the same for both models.

Now let's look at summaries and, \textit{since the models are nested}, a partial F-test to see whether the additional two variables in \verb|lm.AIC| are provide an improvement over \verb|lm.BIC|.

```{r}
summary(lm.AIC)
summary(lm.BIC)

anova(lm.BIC,lm.AIC)
```
Both the $t$-statistics in \verb|summary(lm.AIC)| and the partial $F$-test suggest that neither \verb|SandSaves| nor \verb|PuttsPerRound| is needed.

\textbf{Conclusion:} The smaller model, \verb|lm.BIC|, \verb|log(PrizeMoney) ~ 1 + GIR + BirdieConversion + Scrambling|, seems to be the preferable model.

\textbf{Important Note:} Because we are repeatedly using the data, for variable selection, and to conduct $t$-tests and an $F$-test, the results of the tests may be misleading because of capitalization on chance.  The usual result of testing ``significance'' after model selection on the same data is that the results are \textit{more significant} than they should be.  

In this case, we are seeing non-significance for the two variables \verb|SandSaves| and \verb|PuttsPerRound|, despite this tendency toward inflated significance.  So these two variables probably really are non-significant predictors.

\rm

## Problem 2(f) 

Interpret the regression coefficients in the final model. Is it necessary to be
cautious about taking these results to literally?

\em

\textbf{The final model is} \verb|lm.BIC|.

Since we are considering \verb|log(PrizeMoney)| we can interpret coefficient estimates in \verb|summary(lm.BIC)| in terms of percent change in prize money:
\begin{itemize}
\item An increase of 1 unit in \verb|GIR| (1 extra percent of the time the player was able to hit the greens in regulation) is associated with a $100\times0.15658\%\approx 16\%$ increase in \verb|PrizeMoney|.
\item An increase of 1 unit in \verb|BirdieConversion| (1 extra percent of the time the player makes a ``birdie'' or better after hitting the green in regulation) is associated with a $100\times0.20625\%\approx 21\%$ increase in \verb|PrizeMoney|.
\item An increase of 1 unit in \verb|Scrambling| (1 extra percent of the time the player misses the green in regulation but still makes ``par'') is associated with a $100\times0.09178\%\approx 9\%$ increase in \verb|PrizeMoney|.
\end{itemize}
So it looks like all three of \verb|GIR|, \verb|BirdieConversion| and \verb|Scrambling| are importantly associated with a golfer's prize money.  We don't know if increases in \verb|BirdieConversion| and \verb|Scrambling| \textit{cause} a player to earn more prize money, or if there is a third lurking variable (like natural talent) that influences all three variables, \verb|BirdieConversion|, \verb|Scrambling| and \verb|PrizeMoney|.

\textbf{Important Note:}  All three variables have small SE's and are highly significant, but just like with testing after variable selection, calculating SE's after variable selection usually results in SE's that are too small, and so predictors look more significant than they actually are.

In this case, all the $p$-values are very very small, so we might guess that these variable probably really are significant predictors.

\rm


# Problem 3: [Based on Gelman \& Hill (2009), p. 51, \#5]  

The subfolder \texttt{beauty} in the hw04 folder in the "Files" area for our course on canvas contains data from Hamermesh and Parker (2005) on student
evaluations of instructors’ beauty and teaching quality for several courses at the
University of Texas. The teaching evaluations were conducted at the end of the
semester, and the beauty judgments were made later, by six students who had
not attended the classes and were not aware of the course evaluations.  Various documents in the folder give background and some variable definitions (some variables are defined in the ".log" file there, others' definitions you will have to deduce from pdf's in the subfolder).

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

## Problem 3(a)

Fit a regression model predicting \texttt{courseevaluation} (average student evaluations) from \texttt{btystdave} (the average of 6 standardized beauty ratings for each instructor) and \texttt{female}.  Then fit the same model with the interaction between \texttt{btystdave} and \texttt{female} added in.
\begin{enumerate}
\item Graph each fitted model on a scatter plot of \texttt{courseevaluation} vs \texttt{btystdave}.  Indicate clearly in the graph what the various parameters in the model represent geometrically.
\item Display the four standard diagnostic plots in R and comment on their features, for each model.  Comment on whether the fit seems adequate from the evidence in these plots, for either model.  In case there are problems with the fit, indicate what they are and how you might improve things.
\item Produce summaries of the two fitted models; comment on the coefficient estimates and their standard errors, and on $R^2$, for each model  Use a partial $F$ test to determine whether the interaction should be kept. \textit{ Your comments should include not only technical points (``$B$'' in the ``$ABA^{-1}$'' metaphor for applied statistics from the course syllabus), but also what it means for understanding how factors may influence course evaluations (``$A^{-1}$'').}   
\end{enumerate}

\em
\ul{\bf Part (i):}
```{r fig.width=8}
beauty <- read.csv("ProfEvaltnsBeautyPublic.csv")
## str(beauty) ## too long to include here...
lm.1 <- lm(courseevaluation ~ btystdave + female,data=beauty)
print(beta.lm.1 <- coef(lm.1))
lm.2 <- lm(courseevaluation ~ btystdave * female,data=beauty)
print(beta.lm.2 <- coef(lm.2))
par(mfrow=c(1,2))
plot(courseevaluation ~ btystdave, data=beauty,main="lm.1",
     col=c("blue","red")[beauty$female+1])
abline(beta.lm.1[1],beta.lm.1[2],col="blue")
abline(beta.lm.1[1]+beta.lm.1[3],beta.lm.1[2],col="red")
legend(0.5,2.5,lty=1,col=c("blue","red"),
       legend=c("male instructor","female instructor"),cex=0.55)
plot(courseevaluation ~ btystdave, data=beauty,main="lm.2",
     col=c("blue","red")[beauty$female+1])
abline(beta.lm.2[1],beta.lm.2[2],col="blue")
abline(beta.lm.2[1]+beta.lm.2[3],beta.lm.2[2]+beta.lm.2[4],col="red")
legend(0.5,2.5,lty=1,col=c("blue","red"),
       legend=c("male instructor","female instructor"),cex=0.55)

```
For the additive model (\verb|lm.1|, the plot on the left), the slope of the line relating \verb|courseevaluation| to \verb|btystdave| is `r round(beta.lm.1[2],2)`, for both male and female instructors. The intercept for male instructors is  `r round(beta.lm.1[1],2)`, and for female instructors it's `r round(beta.lm.1[1]+beta.lm.1[3],2)`.  In this model, male instructors get about a `r round(-beta.lm.1[3],2)` boost in course evaluations, vs. female instructors.

For the interactive model (\verb|lm.2|, the plot on the right), the slope of the line relating  \verb|courseevaluation| to \verb|btystdave| is `r round(beta.lm.2[2],2)` for male instructors and `r round(beta.lm.2[2]+ beta.lm.2[4],2)` for female instructors.  The intercept for male instructors is `r round(beta.lm.2[1],2)` whereas for female instructors it is `r round(beta.lm.2[1]+beta.lm.2[3],2)`.  Thus according to this model, evaluations for courses with male instructors benefit by about `r round(-beta.lm.2[3],2)` to start with, and the slope as a function of \verb|btystave| is about `r round(-beta.lm.2[4],2)` higher for male instructors, than for female instructors.

\ul{\bf Part (ii):}
```{r fig.width=8,fig.height=2.5}
par(mfrow=c(1,4))
plot(lm.1)
plot(lm.2)
```
There is very little difference between the two sets of casewise diagnostic plots. In both models, there is no functional pattern in the raw residual plot, slightly short tails in the residuals in the Normal Q-Q plots, approximately constant variance in the Scale-Location plots, and very similar, and unconcerning residuals vs. leverage plots.

So, the fits of both models seem adequate (except perhaps for short tails in the Normal Q-Q plots, which don't concern me very much), and I would not do any transformations, etc. to improve either fit.

About the only difference I see is that the dsitribution of fitted values appear to be slightly higher for \verb|lm.2| than \verb|lm.1| (e.g. see the x-axis of the residuals vs fitted plots).

\ul{\bf Part (iii):}
```{r}
summary(lm.1)

summary(lm.2)

anova(lm.1,lm.2)
```
The most noticable thing about both models is that their $R^2$'s are pitiful: 0.066 for \verb|lm.1| and 0.073 for \verb|lm.2|.  So, even though we have some significant predictors in both models, they are not doing a very good job of accounting for the variation in \verb|courseevaluation|.

In both models, \verb|btystave| and \verb|female| are highly significant predictors, and in \verb|lm.2|, the interaction term \verb|btystdave:female| is just barely non-significant ($p$-value $=$ 0.0789).

Because we are only adding one more variable, the interaction term, the $F$-test statistic is just the square of the $t$-statistic for the interaction in the \verb|summary()| output, and hence there is no new information: the $p$-value is again 0.0789.  So, we would probably prefer the simpler model \verb|lm.1|.

In practical terms, model \verb|lm.1| says that a course with a male instructor starts off with a course evaluation of about `r round(beta.lm.1[1],2)`, but this drops to about `r round(beta.lm.1[1]+ beta.lm.1[3],2)` if the course has a female instructor.  The increase in course evaluation per unit increase in beauty rating is the same for both male and female instructors, about `r round(beta.lm.1[2],2)`.

The standard errors for the coefficients are all quite small (but perhaps a bit smaller than they should be since we used the data twice: once to select a model and once to estimate coefficients), and generally they are reassuring that the effects quoted in the last paragraph are "real".

\rm

## Problem 3(b)

Now let's look at \textit{ all} of the variables in the data set.  Should any of the variables in the data set be transformed before being used in a regression model?  List each variable that is not a dummy variable, and for each of these, 
\begin{itemize}
\item Say whether the variable should be transformed (yes or no)
\item If yes, indicate what transformation you would make
\item Justify these two answers, using both evidence from the data and other considerations
\end{itemize}
Note: being able to communicate with a client or collaborator matters, so there may be instances where either (a) a transformation might help, but you decide against it since it would be difficult to explain to a client/collaborator, or (b) an automatic method like Box-Cox might suggest one power, but you pick a simpler power ``nearby'' because it is easier to explain to a collaborator/client.

\em

Since we're not going to consider \texttt{profnumber}, \texttt{multipleclass} and \texttt{class1} through \texttt{class30} in part (c), I will eliminate them now before we consider transformations.

I will use Box-Cox to suggest transformations, and then choose more interpretable transformations based on Box-Cox.  
```{r}
names(beauty)
prof.loc <- grep("profnumber",names(beauty))
multiclass.loc <- grep("multipleclass",names(beauty))
class.locs <- grep("class",names(beauty))
beauty.red <- beauty[,-c(prof.loc,multiclass.loc,class.locs)]
```
From the \verb|str()| command, our background knowledge of college, and some examination of the papers and variable glossary in the \verb|beauty| subdirectory, we find
\small
\begin{description}\itemsep=0pt\parsep=0pt\parskip=0pt
\item[tenured          ]  dummy: instructor tenured? 0=no, 1=yes
\item[minority         ]  dummy: instructor minority? 0=no, 1=yes
\item[age              ]  continuous: 36 59 51 40 31 62 33 51 33 47 ...
\item[beautyf2upper    ]  continuous: 6 2 5 4 9 5 5 6 5 6 ...
\item[beautyflowerdiv  ]  continuous: 5 4 5 2 7 6 4 4 3 5 ...
\item[beautyfupperdiv  ]  continuous: 7 4 2 5 9 6 4 6 7 7 ...
\item[beautym2upper    ]  continuous: 6 3 3 2 6 6 4 3 5 6 ...
\item[beautymlowerdiv  ]  continuous: 2 2 2 3 7 5 4 2 5 3 ...
\item[beautymupperdiv  ]  continuous: 4 3 3 3 6 5 4 3 3 6 ...
\item[btystdave        ]  continuous: 0.202 -0.826 -0.66 -0.766 1.421 ...
\item[btystdf2u        ]  continuous: 0.289 -1.619 -0.188 -0.665 1.721 ...
\item[btystdfl         ]  continuous: 0.458 -0.0735 0.458 -1.1365 1.521 ...
\item[btystdfu         ]  continuous: 0.8758 -0.577 -1.5456 -0.0927 1.8444 ...
\item[btystdm2u        ]  continuous: 0.682 -1.132 -1.132 -1.736 0.682 ...
\item[btystdml         ]  continuous: -0.9 -0.9 -0.9 -0.313 2.038 ...
\item[btystdmu         ]  continuous: -0.195 -0.655 -0.655 -0.655 0.723 ...
\item[courseevaluation ]  continuous: 4.3 4.5 3.7 4.3 4.4 4.2 4 3.4 4.5 3.9 ...
\item[didevaluation    ]  continuous: 24 17 55 40 42 182 33 25 48 16 ...
\item[female           ]  dummy: instructor female? 0=no, 1=yes
\item[formal           ]  dummy: web pic of instructor wears tie/jacket (dress)? 0=no, 1=yes
\item[fulldept         ]  dummy: everyone in dept has a web pic? 0=no, 1=yes
\item[lower            ]  dummy: lower division (freshman/sophomre)? 0=no, 1=yes
\item[nonenglish       ]  dummy: instructor non-native English speaker? 0=no, 1=yes
\item[onecredit        ]  dummy: a one-credit course? 0=no, 1=yes
\item[percentevaluating]  continuous: 55.8 85 100 87 87.5 ...
\item[profevaluation   ]  continuous: 4.7 4.6 4.1 4.5 4.8 4.4 4.4 3.4 4.8 4 ...
\item[students         ]  continuous: 43 20 55 46 48 282 41 41 60 19 ...
\item[tenuretrack      ]  dummy: is instructor in tenure-track? 0=no, 1=yes
\item[blkandwhite      ]  dummy: (don't know what this is!)
\item[btystdvariance   ]  continuous: 2.13 1.39 2.54 1.76 1.69 ...
\item[btystdavepos     ]  continuous: 0.202 0 0 0 1.421 ...
\item[btystdaveneg     ]  continuous: 0 -0.826 -0.66 -0.766 0 ...
\end{description}
\normalsize

There are a few relationships we can see or guess immediately:
```{r}
attach(beauty.red)
sum(abs(btystdave - (btystdf2u+btystdfl+btystdfu+btystdm2u+btystdml+btystdmu)/6))
sum(abs(btystdave - (btystdavepos+btystdaveneg)))
detach()
```
So,up to rounding error, we see that 
\begin{itemize}
\item \verb|btystdave| is in fact \verb|(btystdf2u+btystdfl+btystdfu+btystdm2u+btystdml+btystdmu)/6|
\item \verb|btystdave| is also \verb|btystdavepos+btystdaveneg|
\end{itemize}
So we will keep \verb|btystdave| and eliminate the variables that sum up to it.  (There is no loss in eliminating \verb|btystdavepos| and \verb|btystdaveneg|, since these are just the "positive part" and "negative part" of \verb|btystdave|. There is a little loss in eliminating \verb|btystdf2u, btystdfl, btystdfu, btystdm2u, btystdml| and \verb|btystdmu|, but hopefully not too much since we will keep the unstandardized versions.)
```{r}
btystd.locs <- grep("btystd",names(beauty.red))[-c(1,8)] 
## we are keeping btystdave, which comes first in this list,
## and btystdvariance, which is eighth...
beauty.red <- beauty.red[,-btystd.locs]
names(beauty.red)
```
The variables we want to consider for transformation are all of the non-dummies that are left:
```{r fig.width=8,fig.height=10}
trans.names <- names(beauty.red)[c(3:12,19:21,24)]

powers <- NULL
for (i in trans.names) {
  x <- beauty.red[,i]
  if (min(x)==0) { 
    x <- x + 0.01*max(x)
    } else if (min(x)<0) {
      x <- x - min(x)*1.01
      }
  ## forces each variable to be strictly positive, for Box-Cox
  powers <- c(powers,powerTransform(x)$lambda)
}
names(powers) <- trans.names
powers <- round(powers,2)
powers

par(mfrow=c(5,3))
for (i in trans.names) {
  plot(density(beauty.red[,i]),main=i,xlab=paste("Suggested power =",powers[i]))
}
    
```

I'm going to round these suggested powers as follows:

```{r}
bjpowers <- powers
bjpowers["age"] <- 1 ## already fairly symmetric
bjpowers["beautyf2upper"]   <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautyflowerdiv"] <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautyfupperdiv"] <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautym2upper"]   <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautymlowerdiv"] <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["beautymupperdiv"] <- 0.5  ## right skewed, sqrt() easier to talk about 
bjpowers["btystdave"] <- 1 ## some right skew but reasonably bounded between -2 and +2
bjpowers["courseevaluation"] <- 1   ## even though 2.2 is suggested, I'm leaving this 
                                    ## alone, to make the linear regression coefficients
                                    ## more interpretable
bjpowers["didevaluation"]   <- -0.5 ## I may try a log here too.
bjpowers["percentevaluating"] <- 2  ## close, perhaps more interpretable
bjpowers["profevaluation"]   <- 3   ## I may try 2 later..
bjpowers["students"]        <- -0.5 ## I may try a log here too.
bjpowers["btystdvariance"]  <- 0.5  ## 0.29 is called for but hard to talk about
```


```{r fig.width=8,fig.height=10}
par(mfrow=c(5,3))
for (i in trans.names) {
  plot(density(beauty.red[,i]^bjpowers[i]),main=paste(i,bjpowers[i],sep="^"),xlab="")
}

```
We'll go with these transformations, though you can see in my comments above that some others might be better\ldots

\rm

## Problem 3(c)

Fit the model that regresses \texttt{courseevaluation} onto all other variables, except for \texttt{profnumber}, \texttt{multipleclass}, and the 30 \texttt{class} variables (\texttt{class1} through \texttt{class30}).  Use the transformations you recommended in part (b).  Make a table indicating 
\begin{itemize}
\item The t-statistics for each variable
\item The VIFs for each variable
\end{itemize}
in your model.

\em

```{r}
beauty.trans <- beauty.red
for (i in trans.names) {
  beauty.trans[,i] <- beauty.red[,i]^bjpowers[i]
  names(beauty.trans)[grep(i,names(beauty.trans))] <- paste("t",i,sep=".")
}
names(beauty.trans)[grep("t.btystdave",names(beauty.trans))] <- "btystdave" 
## didn't transform btystdave...
names(beauty.trans)[grep("t.courseevaluation",names(beauty.trans))] <- "courseevaluation" 
## didn't transform courseevaluation...
names(beauty.trans)[grep("t.age",names(beauty.trans))] <- "age" 
## didn't transform age...
names(beauty.trans)

lm.1 <- lm(courseevaluation ~ ., data=beauty.trans)
tab <- cbind(summary(lm.1)$coef, vif=c(NA,vif(lm.1)))

round(tab,4)
```


\rm

## Problem 3(d)

On the basis of this table, and what you know about the definitions of the variables, would you eliminate any variables in your model?  Why or why not?

\em

\begin{itemize}

\item Since none of them are individually significant, all have moderately large vif's, and they are approximately summarized by \verb|btystdave|, I would eliminate the six students' ratings \verb|t.beautyf2upper, t.beautyflowerdiv, t.beautyfupperdiv, t.beautym2upper, t.beautymlowerdiv| and \verb|t.beautymupperdiv|.
\item Even though \verb|t.profevaluation| is highly significant with a low vif, I would probably remove it.  In my experience, evalation of the instructor and evaluation of the course are highly correlated, and so this variable is probably soaking up variation in course evaluation that I'd like to see explained by other variables in the model.
\item The variables \verb|t.didevaluation,  t.percentevaluating|, and  \verb|t.students| are interesting: In untransformed form there is  clear (nonlinear) relationship among these variables: \texttt{percentevaluating = 100*didevaluation/students}.  Since \verb|t.percentevaluating| has almost no practical effect on course evaluation, I will eliminate that one, and see what the effect is.
\item I would not get rid of \verb|btystave| even though it has a sky-high vif and is only marginally significant: I'd hope that removing the other variables above would allow us to see the effect of \verb|btystdave| on course evaluation.
\end{itemize}
In a "first pass", this is all the farther I'd go; removing these variables is going to strongly affect the $t$ statistics and the vifs.  If I refit the model without these variables, I get
```{r}
lm.2 <- update(lm.1, . ~ . - t.beautyf2upper - t.beautyflowerdiv - t.beautyfupperdiv -
                 t.beautym2upper - t.beautymlowerdiv - t.beautymupperdiv - 
                 t.profevaluation - t.percentevaluating)
tab <- cbind(summary(lm.2)$coef, vif=c(NA,vif(lm.2)))

round(tab,4)
```
Things are starting to come into focus: 
\begin{itemize}\renewcommand{\labelitemii}{$\circ$}
\item \verb|bytstdave| is a significant predictor now
\item Most of the predictors are significant or very nearly significant; the exceptions are \verb|tenured,  lower| and \verb| t.btystdvariance|.  (I would probably try to remove these in a further round of variable selection\ldots)
\item The variables \verb|t.didevaluation| and \verb|t.students| still have somewhat high vifs, but both are significant or nearly so.  I might try to take one of these out (I think they are just collinear with each other: their correlation is  `r round(cor(beauty.trans$t.didevaluation,beauty.trans$t.students),2)`), or I might try to leave them in and see if the rather large effects they have are ``real''.
\item The other variables have coefficients with signs that we might expect:
\begin{itemize}
\item minority instructors have somewhat lower course evaluations ($\hat\beta_{minority}=-0.1584$)
\item each year older an instructor is, the course evaluation lowers a bit ($\hat\beta_{age}=-0.0073$)
\item increasing beauty rating has a somewhat positive effect on course rating ($\hat\beta_{btystdave}=-0.0915$)
\item keeping in mind that \verb|t.didevaluation = 1/sqrt(didevaluation)|, the \ul{fewer} students that evaluate the class, the lower the course rating ($\hat\beta_{t.didevaluation}=-1.5430$)
\item similarly, since \verb|t.students = 1/sqrt(students)|, the \ul{fewer} students in the class, the higher the course rating ($\hat\beta_{t.students}=2.6575$)
\item female instructors take a hit a bit larger than minority instructors in course evaluation ($\hat\beta_{female}=-0.1825$)
\item If the instructor dresses formally for their web photo, or if the department faculty all have web photos, course evaluation increases enough to essentially offset the effect of being female or minority ($\hat\beta_{formal}=0.1426$, $\hat\beta_{fulldept}=0.1997$)
\item If the instructor is not a native speaker of English, the course evaluation takes a hit that's about twice as big as the minority instructor hit ($\hat\beta_{nonenglish}=-0.3116$)
\item One credit course get over a half-point advantage in course evaluation ($\hat\beta_{onecredit}=0.5339$)
\item We have no idea what \verb|blkandwhite| is, but it has a positive effect on course evaluation! ($\hat\beta_{onecredit}=0.2023$)
\end{itemize}
Of course, all the interpretations above sound like we are saying each predictor ``causes'' a change in course evaluation, but we really don't know about lurking variables, common causes, etc., so we can't be sure of ``causes'' here.
\end{itemize}


\rm

## Problem 3(e)

Why might the methods used in parts (c) and (d) not be adequate for deciding which variables to keep, and which ones to eliminate, in a regression model?

\em

Essentially, the question is, "what could go wrong with selecting variables on the basis of $t$-statistics and vif's?"

There are a few potential difficulties with this approach:
\begin{enumerate}\renewcommand{\labelenumi}{\arabic{enumi}.}
\item Removing variables one-at-a-time can change other coefficients and $t$-statistics, so if you remove the variables in one order based on $t$-statistics, and I remove them in another order, we could end up with very different models.
\item High vif's do not mean that the model is invalid; they just mean that the estimated coefficients will be hard to interpret.  Removing variables on the basis of vif's can remove important predictor variables from a model
\item One-variable-at-a-time approaches are examples of ``greedy'' algorithms: since they do not consider all variables at each step of variable selection, they can miss the very best models (forward selection with BIC in Problem \#2 is an example of this).
\end{enumerate}
In the case of this problem, these dificulties are offset by "subject matter knowledge": since we know a lot about college or university life, we can make better guesses about what order to remove variables in, what groups of variables to remove all at once, etc.

Just for kicks, I tried "all subsets" variable selection on the model fitted in part (c) here (except that I took \verb|profevaluation| out, because of its close relationship with \verb|courseevaluation|).  This approach uses no subject matter knowledge (none of our knowledge about university life!), and replaces it with just mathematical optimality:

```{r}
all.subsets <- regsubsets(courseevaluation ~ . - t.profevaluation, data=beauty.trans,
                          nvmax=dim(beauty.trans)-1)  # -1 for t.profevaluation
tmp <- summary(all.subsets)
attach(tmp)
p <- 1:dim(tmp$which)[1]
n <- dim(beauty.trans)[1]
results <- data.frame(which, bic=n*log(rss)+log(n)*(p+2), aic=n*log(rss)+2*(p+2))
detach()
```


\tiny

```{r}
oldwidth <- options()$width
options(width=200)
results
minimize(results,"bic")
minimize(results,"aic")
options(width=oldwidth)
```

\normalsize

The best BIC model here is

\begin{verbatim}
courseevaluation ~ 1 +  t.beautyfupperdiv +  female + nonenglish + onecredit +
    t.percentevaluating + blkandwhite
\end{verbatim}

and the best AIC model is

\begin{verbatim}
courseevaluation ~ 1 + minority + t.beautyf2upper + t.beautyfupperdiv + t.beautym2upper +
     t.didevaluation + female + nonenglish + onecredit + t.percentevaluating +
     blkandwhite
\end{verbatim}

This turns out to be a nice illustration of how mathematical optimality may not be what makes the most sense substantively:
\begin{itemize}
\item Each model tries to make a new ``beauty'' variable instead of taking \verb|btystdave|.  There isn't really much substantive sense to the beauty variables that the models pick out.
\item The other variables in these models are a subset of the model from part (d).  They are fine, but they also seem to be missing significant predictors with intepretable coefficients.  
\end{itemize}
If I were advising the university about what (besides teaching quality) affects course evaluation, I would not want to miss some of the variables that the optimal AIC and BIC models miss!
\rm


