---
title: "Homework 05 Solutions"
date: "10/4/2021"
output: 
  pdf_document:
    keep_tex: true
fig_caption: true
header-includes:
- \usepackage{amsmath}
- \usepackage{amsthm}
- \usepackage{xcolor}
- \usepackage{txfonts}
- \DeclareMathOperator{\Cov}{Cov}
- \DeclareMathOperator{\Corr}{Corr}
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, include = TRUE, message = FALSE)
library(car)
library(GGally)
library(knitr)
library(pander)
library(tidyverse)
```

# 1. (Sheather 6.7.5)
```{r 1}
golf <- read_csv("pgatour2006.csv")[, c(3, 5:10, 12)]
# plot(golf)
# golf %>% mutate(PrizeMoney = log(PrizeMoney)) %>% plot
```

```{r}
m1 <- lm(log(PrizeMoney) ~ ., golf)
```

(a) Agree. Since the $Y$ variable is heavily right-skewed and ranges over several orders of magnitude, a log transformation probably makes sense. To confirm this, you can regress `PrizeMoney` against the other predictors and then regress `log(PrizeMoney)` against the other predictors. The residuals are clearly non-normal in the first case and are approximately normal in the second case, which suggests that the log transformation of the outcome variable is appropriate. There is no apparent need to transform the predictors.

(b) I simply regressed the log-transformed `PrizeMoney` variable (`LogPrizeMoney`) against the seven predictors. The diagnostic plots (Figure \ref{f:diagnostics1}) look about as good as one could ask. Since I know nothing about golf, I have no reason to expect any particular interactions, so for the sake of simplicity, and since the diagnostics look good, I will leave the model as is. A summary of the model is below, and Figure \ref{f:golf_pairs} contains a pairs plot of the data. Some collinearity is evident in the pairs plots, in particular between `PuttingAverage` and `BirdieConversion` and between `PuttingAverage` and `PuttsPerRound`.

> It is tempting to remove some of the predictors with high p-values, but until we are working with a golf expert as a collaborator we do not know which ones to remove.  If we wanted to remove some predictors without further information, we might consider removing ones for which the 95% confidence interval for $\hat\beta$ (a) is narrow and (b) contains zero.  `DrivingAccuracy`, `SandSaves` and maybe `Scrambling` could be considered for removal.  On the other hand `PuttingAverage` has a very wide confidence interval, which suggests that we may not have enough information yet to have a good estimate of its $\hat\beta$.  (See also the comments after part (e) below.)


## Summary of the golf model:
```{r 1_summary}
pander(summary(m1), caption = "")
```

```{r 1_pairs, fig.cap="Pairs plot of golf data, with log transformation of PrizeMoney variable.\\label{f:golf_pairs}"}
# golf <- golf %>%
#   mutate(LogPrizeMoney = log(PrizeMoney)) %>%
#   select(-PrizeMoney)
ggpairs(golf %>% mutate(PrizeMoney = log(PrizeMoney)))
```

```{r 1_diagnostics, fig.cap="Problem 1: Diagnostic plots.\\label{f:diagnostics1}"}
par(mfrow=c(2, 2), mar = c(4, 4.5, 2, 1))
plot(m1)
```

(c) There are no points that clearly require investigation. Point 185 has a large residual, but in a data set of size 196 it is not particularly surprising to see a point with a standardized residual of 3 or so. There are no bad leverage points, and since the diagnostics look good in general, there's nothing obvious to investigate.  

(d) As noted above, there's some collinearity in the predictors, so some of the predictors may be redundant; and the standard errors of the estimates may be large because of this. As always, the assumption that the outcome is linearly related to the predictors is a strong one, and it means the model will not capture nonlinearities.

(e) There are two problems with relying only on insignificant $t$-values to eliminate predictors in a model.  The first is that a $t$-value can be small because the predictor is (nearly) collinear with another predictor or group of predictors.  In that case, it's not obvious which predictor (if any!) to eliminate, without a conversation with a subject matter expert.  The second is that hypothesis tests are asymmetrical, in the sense that when we reject the test, we are making an assertion that $H_0$ is false (with the attendant possibility of Type I error), whereas when we fail to reject the test, it could be because we have strong evidence that $H_0$ is true or because we simply don't have enough evidence to make a judgment either way. 

>> A TA from a previous iteration of the class comments further on the second reason:

>> The value of the $t$-statistic associated with a coefficient $\beta_k$ in a regression is a function of the magnitude of $\beta_k$, the variance of the predictor $X_k$, the amount of noise around the regression line, and the sample size. We may be able to reject the null for coefficients with very small magnitudes if, say, the variance of the associated predictor is high; conversely, we may fail to reject the null for coefficients with large magnitudes because the range of the associated predictor is small. Just because a test is insignificant does not mean that the predictor is insignificant in a real world sense. (If, however, the confidence interval for the predictor were very small and contained 0, then we might be able to confidently say that the coefficient was insignificant in a real world sense, as long as the model were generally credible.) For more details, see Professor Shalizi's lecture notes from 36-401:
\begin{itemize}
  {\color{blue} \item \href{http://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/08/lecture-08.pdf}{Fall 2015, Lecture 8. Statistical Signifcance: Uses and Abuses}} (esp. Sections 5, 6, \& 7)
  {\color{blue}\item \href{http://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/lecture-15.pdf}{Fall 2015, Lecture 15. What, Exactly, Is R Testing?}}
  {\color{blue}\item \href{http://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/26/lecture-26.pdf}{Why Variable Selection Using p-Values Is a Bad Idea}}
\end{itemize}


# 2.  The Beauty Data Again

## (a)
Why might we want to transform a predictor? There are at least three reasons:
\begin{itemize}
  \item For interpretability, for example if we want the intercept to have a particular meaning, or if we want to be able to think about percentage changes in $y$ relative to $x$, rather than linear changes.
  \item To reduce the leverage of particular points.
  \item To satisfy modeling assumptions, for example if $y$ is not linearly related to $x$ but is linearly related to some transformation $g(x)$.
\end{itemize}

Now let's think about the predictors we have in this problem. The distributions of all the non-indicator predictors are given in Figure \ref{f:hist}.

Regarding the first item 1, I see no compelling reason to transform the predictors for the sake of interpretability. True, a y-value with respect to an "age" of 0, for example, will not be meaningful, but we won't ever be predicting in this range anyway.

With respect to the second item above, there are some predictors with very skewed distributions, so we could consider transforming them to reduce the skewness. In particular, `btystdvariance`, `didevaluation`, and `students` are strongly right-skewed and strictly positive, so we could consider log-transforming them.

We also have `btystdavepos`, which is strongly right-skewed, and `btystdaveneg`, which is strongly left-skewed. However, both these variables have a large number of values which are 0, so obviously no deterministic transformation will change the fact that there are a large number of instances with the same value. Since the rest of the values are fairly spread out within their ranges, transformations won't necessarily be useful.

With respect to the third item above, we'd need to fit the model and then check whether the assumptions appear reasonable.

```{r 2a, fig.cap="Problem 2: Histograms of all non-indicator variables in the beauty dataset.\\label{f:hist}", fig.height=7}
beauty <- read_csv("ProfEvaltnsBeautyPublic.csv") %>%
  select(-profnumber, -multipleclass, -starts_with("class"))

is_ind <- function(var) {all(na.omit(var) %in% c(0, 1))}
indicator_vars <- names(beauty)[sapply(beauty, is_ind)]

beauty_long <- gather(select(beauty, -indicator_vars))
ggplot(beauty_long, aes(x = value)) +
  geom_histogram() +
  facet_wrap(~ key, scales = "free")
```

## (b)
I tried fitting the model with no transformations as well as with log transformations of `btystdvariance`, `didevaluation`, and `students`. The diagnostics for the former (Figure \ref{f:diagnostics2}) looked slightly better than the latter, so I will stick with the model with no transformations. The t-values and VIF values are in Table \ref{t:t_vif}.

```{r 2b_summary}
m2 <- lm(courseevaluation ~ ., data = beauty)
summary(pander(coef(m2), caption = ""))
```

```{r 2b_table}
out <- cbind(coef(summary(m2))[, 3], c(0, vif(m2)))
colnames(out) <- c("t-value", "VIF")
out <- out[order(rownames(out)), ]
kable(out, digits = 3, caption = "Problem 2b: Table of t-values and Variance Inflation Factors for the model coefficients.\\label{t:t_vif}")
```

```{r 2b_diagnostics, fig.cap="Problem 2b: Diagnostics.\\label{f:diagnostics2}"}
par(mfrow=c(2, 2), mar = c(4, 4.5, 2, 1))
plot(m2)
```

(c) Using the rule of thumb that a VIF above 5 is high, many of the predictors have extremely high VIFs. If the purpose of the model is to perform inference on the coefficients, then you may wish to eliminate some of these predictors. If the purpose of the model is to generate accurate predictions, then you may wish to leave them in, though we would need to evaluate predictive performance on various subsets of predictors in order to make that decision.

(d) The same issues considered in part (e) of the previous problem apply here.  The two reasons to not consider $t$-statistics and VIF's alone are:
\begin{itemize}
\item $t$-statistics tell you about whether a variable has a coefficient significantly different from zero, {\em after} accounting for all of the other predictors.  If the variable with a non-significant $t$-statistic is collinear with another variable or group of variables, it is no obvious which (if any!) predictors to eliminate, without help from a subject matter expert.
\item The $t$-statistic is being used in a hypothesis test of the hypothesis $H_0:~\beta=0$.  If the $t$-statistic is large, we have sufficient evidence to reject $H_0$. But if the $t$-statistic is small, it might be {\em either} because the $\beta$ really is zero, {\em or} because we don't have enough evidence (data) yet to rule out $\beta=0$. 
\end{itemize}

<!---
(d) Table \ref{t:t_vif2} contains the t-values and VIF values for a model with `profevaluation` as the dependent variable. The VIFs are virtually unchanged, which is unsurprising: the VIF for a given predictor is calculated by regressing that predictor on all the other predictors, so we wouldn't expect the results to change much when we swap out one predictor among a large set of variables.

The t-values, on the other hand, are quite different, which is also not surprising. However, I would only rarely use the t-values to inform variable selection. Since the VIFs haven't changed much, the decision about which variables to include or exclude wouldn't change either. 

```{r 2d_table}
m3 <- lm(profevaluation ~ ., data = beauty)
out <- cbind(coef(summary(m3))[, 3], c(0, vif(m3)))
colnames(out) <- c("t-value", "VIF")
out <- out[order(rownames(out)), ]
kable(out, digits = 3, caption = "Problem 2d: Table of t-values and Variance Inflation Factors for the model coefficients.\\label{t:t_vif2}")
```

-->