**Question**  Find all variables associated with an excessively large VIF. How should the VIFs be interpreted? (Should you remove *all* of these variables from your model? Discuss possible strategies to decide on the variables to exclude from your model due to the large VIF values.)

<!--

Answer Q5. x4 and x5 have VIFs over 10 (the threshold). Values over 10 mean than if we fit the linear regression like x4 = a0 + a1*x1 + a2*x2 + a3*x3 + a4*x5, we'll get an R^2 of over 90%. (Since $VIF=\frac{1}{1-R^2_4}$ where the R^2_4 is for this special regression for x4. Solve for R^2_4 when VIF=10 and you'll get R^2=90%.)

It makes sense to remove at least one variable since at least one must be close to redundant given the other regressors. (i.e. there is multi-collinearity present). Given there are two values with nearly the same large VIF, it makes sense to take only one out at a time. (We also know this makes sense since we simulated the data...but we can always check the VIFs after one regressor is removed and if the VIF for the other variable is still large we can then take out the second variable.)

-->

The code chunk below fits a *modified* version of Model 1 and saves the result in an object named *MM1*.

```{r}
MM1 <- lm(formula=y~x1+x2+x3+x4, data=df) # "Modified" Model 1
tidy.MM1 <- tidy(MM1)
glance.MM1 <- glance(MM1)
```

Look at the output from the Modified Model 1. 

```{r echo=FALSE, eval=FALSE}
tidy.MM1  # not included for students
glance.MM1  # not included for students
```

**Question**  What is the value of $p$ associated with Modified Model 1?

<!--

Answer Q6. $p$ is the number of regression coefficients (number of regressors), including the intercept (including a constant regressor). So here $p=5$.

-->

**Question**  According to the Modified Model 1 output, which regression coefficients are significant at the 5% level? Do the estimated coefficients seem "close" to the true values used to simulate the data? 

<!--

Answer Q7. The estimated intercept and coefficient on x2 are not (individually) significant. The estimated values are not spot on, but they seem closer to the corresponding true values than under Model 1. 

-->

**Question**  Note the size of the standard errors associated with each estimated coefficient. How do these compare with the coefficients that were not significant in Model 1? (See Q2.) Given the way the data were simulated, do you have any ideas about why there might still be insignificant coefficients, even for the Modified Model 1? 

<!--

See the *glance.MM1* output. The standard errors for the remaining coefficients haven't changed all that much, but we have removed the coefficient that had the very large standard error (x5). Multi-collinearity will tend to inflate the standard errors.

Why still relatively large standard errors? Well, there is still a small sample size here! It seems that the "signal" for the intercept, and for x2, are too small to be detected amid all of the noise in the data. 

-->

**Question**  What are the values of: $R^2$ and the Adjusted $R^2$ values for the fitted Modified Model 1? Compare these values to those from the original Model 1. (See Q4.) Do you think removing regressor $x_5$ has improved or diminished the model?

<!-- 

Answer to Q9. Need to look at *glance.MM1* to see the result: r.squared = 0.5255 or 52.55% and adjusted r.squared = 0.4833 or 48.33%. We really can only compare the adjusted r.squared since the ordinary r.squared doesn't account for the number of regressors. We can see that the adjusted r.squared actually improves with the removal of x5, so the models seems better without x5 in it (on this measure, at least). 

-->
    
Next, calculate the variance inflation factors (VIF) associated with the regressors in Modified Model 1.

```{r}
vif(MM1)
```



Change the value of $p$ in the code chunk below to match the correct value for Modified Model 1. Then run the code chunk to extract the *leverage* and *Cook's D* measures associated with the fit of Modified Model 1 that exceed the "rule of thumb" threshold of $2*p/n$. 

```{r fig.height=5, fig.width=8, eval=FALSE}
df <- df %>% mutate(rowid=1:n)
aug.MM1 <- augment(MM1)
aug.MM1 <- aug.MM1 %>% mutate(rowid=1:n)

p <- 1 # change this once you know the correct value of p

threshold <- 2*p/n # threshold for .hat and .cooksd

exceed_hat <- aug.MM1 %>% 
  arrange(desc(.hat)) %>% 
  select(rowid, .hat) %>% 
  filter(.hat > threshold) 

exceed_cooksd <- aug.MM1 %>% 
  arrange(desc(.cooksd)) %>% 
  select(rowid, .cooksd) %>% 
  filter(.cooksd > threshold) 

df_longer <- df %>% 
  pivot_longer(-c(y, rowid), 
               names_to="regressor", 
               values_to="factor")  

df_longer %>% 
  ggplot(aes(x=factor, y=y, colour=regressor)) + 
  geom_point() + 
  facet_wrap(~regressor, nrow=2) + 
  theme_bw() + 
  geom_text(label=df_longer$rowid, 
            nudge_x = 3, 
            colour="black", 
            size=2,
            check_overlap = T) + 
  ggtitle("y against each regressor") 
```

<!-- code chunk below has the correct value of p (not shown to students) -->

```{r fig.height=5, fig.width=8, eval=FALSE, echo=FALSE}
df <- df %>% mutate(rowid=1:n)
aug.MM1 <- augment(MM1)
aug.MM1 <- aug.MM1 %>% mutate(rowid=1:n)

p <- 5 # correct value here

threshold <- 2*p/n 

exceed_hat <- aug.MM1 %>% arrange(desc(.hat)) %>% dplyr::select(rowid, .hat) %>% filter(.hat > threshold) 

exceed_cooksd <- aug.MM1 %>% arrange(desc(.cooksd)) %>% dplyr::select(rowid, .cooksd) %>% filter(.cooksd > threshold) 

df_longer <- df %>% pivot_longer(-c(y, rowid), names_to="regressor", values_to="factor")  

df_longer %>% ggplot(aes(x=factor, y=y, colour=regressor))+ geom_point() + facet_wrap(~regressor, nrow=2) + theme_bw() + geom_text(label=df_longer$rowid, nudge_x = 3, colour="black", size=2, check_overlap = T) + ggtitle("y against each regressor") 
```

**Question**  How many points exceed the leverage threshold? Find these points using their *rowid* value which is also shown on the faceted plot. (Some will be easier to find than others! Try filtering the data points to see the values of the regressors to make it easier to find the points.)

<!-- 

Answer Q11. Need to view the *exceed_hat* object. Then look at the plot produced above. 

-- code chunk below not shown to students -->

```{r echo=FALSE, eval=FALSE}
exceed_hat

df %>% dplyr::filter(rowid == 21 | rowid == 16 | rowid == 44)
```

**Question**  How many points exceed the Cook's D threshold? 

<!-- 

Answer Q11. Need to view the *exceed_cooksd* object. No points are flagged as influential according to this measure. 

-- code chunk below not shown to students -->

```{r eval=FALSE, echo=FALSE}
exceed_cooksd
```