---
title: "Homework 05 Solutions"
date: "2022-10-08"
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.

library(marginalmodelplots) ## because mmps doesn't work for glm's...
library(DHARMa)  ## for better glm residual plots...

library(tidyverse) ## at last...

library(foreign)  ## for "read.dta" for problem 3

```


# Problem #1

Please do ISLR, \#11, pp. 325--326. ({\textit{ Hint: You have seen me do a variation of part (a) many times in class, for demonstrations! Some examples are in the lecture notes on omitted variable bias, for example.}) Among other things, this shows that backfitting (a method used to estimate GAM's) is similar to iteratively adjusting added variable plots(!)\ldots

GAMs are generally fit using
a backfitting approach. The idea behind backfitting is actually quite
simple. We will now explore backfitting in the context of multiple
linear regression.

Suppose that we would like to perform multiple linear regression, but
we do not have software to do so. Instead, we only have software
to perform simple linear regression. Therefore, we take the following
iterative approach: we repeatedly hold all but one coefficient esti-
mate fixed at its current value, and update only that coefficient
estimate using a simple linear regression. The process is continued un-
til convergence—that is, until the coefficient estimates stop changing.

We now try this out on a toy example.

## Problem 1(a) 

Generate a response y and two predictors x1 and X2, with
n = 100.

\em

```{r}
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- 1 + 2*x1 + 3*x2 + 4*rnorm(n)
```


\rm

## Problem 1(b) 

Initialize $\beta_1$ to take on a value of your choice. It does not matter
what value you choose.

\em

```{r}
beta1 <- 10
```

\rm


## Problem 1(c) 

Keeping $\beta_1$ fixed, fit the model
\[
 Y - \beta_1 X1 = \beta_0 + \beta_2X2 + \epsilon.
\]

\em
You can do this as follows:
```{r}
a <- y - beta1 * x1
print(beta2 <- lm(a ~ x2)$coef[2])
```

\rm

## Problem 1(d) 

Keeping $\beta_2$ fixed, fit the model
\[
Y ~ \beta_2 X2 = \beta_0 + \beta_1X1 + \epsilon.
\]

\em

You can do this as follows:
```{r}
a <- y - beta2 * x2
print(beta1 <- lm(a ~ x1)$coef[2])
```

\rm

## Problem 1(e) 

Write a for loop to repeat (c) and (d) 1,000 times. Report the
estimates of $\beta_0$, $\beta_1$, and $\beta_2$ at each iteration of the for loop.
Create a plot in which each of these values is displayed, with $\beta_0$,
$\beta_1$, and $\beta_2$ each shown in a different color.

\em

```{r}
beta0 <- lm(a ~ x1)$coef[1]
b1 <- beta1
b2 <- beta2

m <- 1000
for (r in 2:m) {
  a <- y - b1 * x1
  b2 <- lm(a ~ x2)$coef[2]
  beta2 <- c(beta2,b2)

  a <- y - b2 * x2
  b1 <- lm(a ~ x1)$coef[2]
  beta1 <- c(beta1,b1)
  
  beta0 <- c(beta0,lm(a ~ x1)$coef[1])
}

plot(c(1,m),c(0,max(c(beta0,beta1,beta2))),type="n",
     xlab="iteration",ylab="value")

lines(1:m,beta0,col="blue")
lines(1:m,beta1,col="orange")
lines(1:m,beta2,col="red")

legend (800,0.5,legend=c("beta0","beta1","beta2"),col=c("blue","orange","red"),lty=1,cex=0.5)


```



\rm

## Problem 1(f) 

Compare your answer in (e) to the results of simply performing
multiple linear regression to predict Y using X1 and X2. Use
the abline() function to overlay those multiple linear regression
coefficient estimates on the plot obtained in (e).

\em

Here are the final estiamtes from backfitting:
```{r}
beta0[m]
beta1[m]
beta2[m]
```

And here's a summary of the multiple regression:
```{r}
summary(lm(y ~ x1 + x2))
```
The coefficient estiamtes are practically identical!

\rm

## Problem 1(g) 

On this data set, how many backfitting iterations were required
in order to obtain a “good” approximation to the multiple re-
gression coefficient estimates?

\em


Very few---as little as 2--3, and certainly lesss than 5!

```{r}
head(beta0)
head(beta1)
head(beta2)

m <- 10
plot(c(1,m),c(0,max(c(beta0,beta1,beta2))),type="n",
     xlab="iteration",ylab="value")

lines(1:m,beta0[1:m],col="blue")
lines(1:m,beta1[1:m],col="orange")
lines(1:m,beta2[1:m],col="red")

legend (8,0.5,legend=c("beta0","beta1","beta2"),col=c("blue","orange","red"),lty=1,cex=0.5)

```


\rm

# Problem #2

The \texttt{warpbreaks} data set included in R gives the results of an experiment to determine the the effect of wool type (A or B) and tension (low, medium or high) on the number of warp breaks per loom. Data was collected for nine looms for each combination of settings. You can get more information on this data set from \texttt{help(warpbreaks)}.  Use \texttt{View(warpbreaks)} and \verb|xtabs(breaks ~ wool + tension, data=warpbreaks)| to familiarize yourself with the data.  (You don't have to turn in any of this initial exploration.)

## Familiarizing with data (not to turn in!)

\em

The below EDA shows us few things, but the plot visualization is clearest relative to future exploration we will be doing. We can observe the neither the rows nor the columns inform the viewer of the full relationship between such class structure and the number of breaks of the wool. Specifically note that Wool A appears to break more with loose (L) tension (compared to Wool B), but that this structure (that Wool A breaks more than Wool B) isn't observed / as strong across all levels of tension.

```{r message = F, warning = F}

xtabs(breaks ~ wool + tension, data = warpbreaks)

warpbreaks %>% head

warpbreaks %>% 
  ggplot() +
  geom_histogram(aes(x = breaks)) +
  facet_wrap(wool ~ tension)
```

\rm

## Problem 2(a)

Fit the two Poisson regression models \verb|breaks ~ wool + tension| and \verb|breaks ~ wool * tension| (dont't forget \texttt{family=poisson}!).  Consider \texttt{summary()}'s, \texttt{plot()}'s of the fitted glm's, \texttt{mmplot()}'s and also \texttt{binnedplot()}'s (from \texttt{library(arm)}) and/or plots from \texttt{DHARMa}, and any other plots you think are useful.  Comment on the fits of each model, using the graphical evidence you have obtained \textit{ (if you think some plot(s) are \ul{\textbf{not}} useful for this, please tell which ones \& why)}.

\em

### 2(a) No interaction term

Our anova model without an interaction term between wool and tension in predicting number of breaks doesn't fit the data to well. Although this anova model only provides us with (at most) 6 different fitted values for 54 observations, we can see that the lack of interaction terms impacts the our residual vs fitted plot, where we find that the true average value for each of the Wool + tension combinations does not exactly equal the predicted average (this can be seen by the smoothed regression line of the residual vs fitted not being 0 for all fitted values). For this anova, it's probably not the best to *over-interpret* the standard summary plots (due to the 6 fitted values). A similar observation relative to the residuals vs fitted can be seen in marginal plot, where the fitted curve doesn't match the raw average of the number of breaks verse predicted number of breaks. Note that the x-axis values are of the log(breaks) scale, whereas the y-axis is on the breaks scale.

Our "DHARMa" residual plots present a different way to see that our residuals still relationships with our predicted values (aka the 6 classes between Wool and tension), suggesting that the model that we fit doesn't create / simulate residuals similar to the truely observed value. The over-dispersion seen in the QQ plot of the ranks of the residuals makes sense and it suggests that the true residuals are more variable than the simulation model would expect (which does make sense relative to our observation that the model doesn't well capture the mean of each group).

It should be noted that I did not use the binned residual plot as it would show almost the same information as the regular residual vs fitted plot (relative to 6 bins). The binned residual plot may have captured the information in the `summary`'s residual vs fitted and absolute residual vs fitted plots. We include the code to do so below just in case that would be useful - just commented out.

```{r}
p_glm1 <- glm(breaks~wool+tension, data = warpbreaks,
              family = poisson)

summary(p_glm1)
```

```{r}
par(mfrow = c(2,2))
plot(p_glm1)
```
 
```{r warning = F}
mmplot(p_glm1, locfit.control=list(nn = "gcv")) 
```

```{r}
sim_resid1 <- simulateResiduals(p_glm1)
plot(sim_resid1)
```

```{r}
# library(arm)
binnedplot(predict(p_glm1), resid(p_glm1))
```


### 2(a) Interaction term

For the model with an interaction term (who's diagnostic plots are below), we find the model is able to capture all the different groupings of wool an tensions mean (but design). This is seen in the residual vs fitted and the fact that both smoothed fitted and true lines are very close to each other in overall marginal plot.  Although this TA would discourage from interpreting the standard summary plots too much, we can observe that distribution of errors for the wool + tension group with average log(breaks) of 3.8 has a log-based residuals that are larger than other wool + tension groups (as seen in the residual vs fitted and absolute residual vs fitted plots). The difference between the variability of the residuals per each wool + tension group also seems to lead to overdispersion of the residuals relative to simulations draw from out fitted models (visualized in the "DHARMa" plots).

```{r}
p_glm2 <- glm(breaks~wool*tension, data = warpbreaks,
              family = poisson)
```

```{r}
summary(p_glm2)
anova(p_glm1, p_glm2)
```

```{r}
par(mfrow = c(2,2))
plot(p_glm2)
```

```{r message = F, warning = F}
mmplot(p_glm2)
```
```{r message = F, warning = F}
#library(arm)
binnedplot(predict(p_glm2), resid(p_glm2))
```

```{r message = F, warning = F}
plot(simulateResiduals(p_glm2))
```


\rm

## Problem 2(b)

We would like to know whether the interaction between \texttt{wool} and \texttt{tension} is needed.  
\begin{enumerate}\def\theenumi{\roman{enumi}}
\item Compare AIC's for the two models, and then compare BIC's for the two models.  Do AIC and BIC choose the same model? 
\item Conduct a likelihood ratio test for including the interaction. Does the LR test choose the same model as either AIC or BIC? \textit{ (Hint/recipe: For the likelihood ratio test, use the \texttt{logLik()} function to obtain the log-likelihood of the two models,  compute $-2$ times the difference between the null and alternative model log-liklihoods, and obtain a $p$-value from the tail of an appropriate $\chi^2$ distribution.)} 
\end{enumerate}

\em

### (i)

As can be seen in the below table, the second poisson model with interaction terms has a lower AIC and BIC than the first model. If we remember that both of these score functions look like
\[
-2 \cdot \text{(log likelihood)} + \text{constant}\cdot \text{(number of parameters)}
\]
then we can conclude the both AIC and BIC recommend going with the poisson model **with** interaction terms (aka lower AIC/BIC is better).

```{r}
df_aic_info <- data.frame(AIC = c(AIC(p_glm1), AIC(p_glm2)),
                          BIC = c(BIC(p_glm1), BIC(p_glm2))) %>% t()
colnames(df_aic_info) <- c("poisson model, no interaction",
                           "poisson model, interaction")

knitr::kable(df_aic_info)
```

### (ii)

The LRT also suggests we should reject the null hypothesis that the model with an interaction only explains as much information as the the model without the interaction term. This would lead us to select the second model like AIC and BIC did.

```{r}
diff_log_like <- -2 * (logLik(p_glm1) - logLik(p_glm2))
1-pchisq(diff_log_like, df = 4)
```


\rm

## Problem 2(c)

For the best model from part (b), we would like to know if overdispersion is a problem. Use the Pearson residuals to conduct a test of overdispersion (hint: follow the recipe from lecture \#17.1).  Does this agree with the conclusion of \texttt{testDispersion()} from \texttt{library(DHARMa)}?   

\em

The below p-values suggest that the residuals in model 2 are over-dispersed and match the simulation test from DHARMa.

```{r}
n <- nrow(warpbreaks)

pp2 <- length(coef(p_glm2))

e_y_2 <- predict(p_glm2, type = "response")
z_2 <- (warpbreaks$breaks - e_y_2)/ sqrt(e_y_2)
test_stat_2 <- sum(z_2^2)

pchisq(test_stat_2, n-pp2, lower.tail = F)


```


```{r}
testDispersion(p_glm2, alternative = "greater")
```


\rm

## Problem 2(d)

Refit the model from part (c) using \texttt{family=quasipoisson}, and note any differences in overdispersion parameter, coefficient estimates, SE's, etc.  Are there any differences in which predictors are significant?

\em


```{r}
p_glm2d <- glm(breaks~wool*tension, data = warpbreaks,
              family = quasipoisson)
```

```{r}
summary(p_glm2d)
```

Fitting the second model using quasipoisson allows us to have a disperison parameter that is not 1. The new model's dispersion paramter is now `r round(summary(p_glm2d)$dispersion,4)`, which is much larger than 1 - which does suggest that we had a problem with residuals being overdispersed in the original model.



\rm


# Problem #3

[Based on some exercises in Gelman \& Hill.] The data set \texttt{nes5200\_{}processed\_{}voters\_{}realideo.dta} has data from the National Election Survey in every congressional election year in the US from 1948 to 2002.  You can read it into R using the function \texttt{read.dta()} from \texttt{library(foreign)}.  We are interested in learning whether some of the survey data is useful in predicting respondents' vote for president.

## Problem 3(a)

Extract from the full NES data set above a data frame with just data from the presidential election year 1976, and the following variables: \texttt{ideo7, black, female, rep\_{}pres\_{}intent, presapprov, age, educ1, urban, income, union} and \texttt{perfin1}.  The variable \texttt{rep\_{}pres\_{}intent} is 1 if the respondent intended to vote for the Republican candidate (Gerald Ford, who assumed the presidency when Richard Nixon resigned amid an impeachment investigation), and 0 if the respondent intended to vote for the Democratic candidate (Jimmy Carter).  Familiarize yourself with all of the variables by using \texttt{summary()} on your 1976 data frame\footnote{A command like \texttt{attr(nes\_{}data,"var.labels")} on the full data set will give you brief definitions of all of the variables.} to look at the categories. 
You don't have to turn anything in for this, but you should do it.

\em

Nothing is required to be turned in here.

```{r}
nes_all <- read.dta("nes5200_processed_voters_realideo.dta")
#attr(nes_all, "var.labels")

nes_1976 <- nes_all %>% filter(year == "1976") %>%
  dplyr::select(ideo7, black, female,
                          rep_pres_intent,
                          presapprov, age,
                          educ1, urban, income, union,
                          perfin1)
#summary(nes)
```


\rm

## Problem 3(b)

Fit a logistic regression predicting \texttt{rep\_{}pres\_{}intent} from the other variables, for 1976.  Assess the fit using appropriate graphical methods.  Interpret the results: which variables seem to matter for predicting voters' preference for president?  How do they affect the odds of voting for the Republican candidate?

\em



```{r}
l_glm_1976 <- glm(rep_pres_intent ~ ., data = nes_1976, family = "binomial")
```

### 3(b) Fit assessment

For standard logistic regression, where the data are assumed to be **Bernoulli** random variables (aka 0/1), the standard summary plots are not very helpful in letting us understand the fit. The summary does suggest that the model captures some of the variability that we see. Interestly, the binned residual plot we see that the average logit predicted value is less extreme than it should be at either end the spectrum (extremely likely to vote Democrat or Republican). At the same time, the DHARMa residual diagnostics suggest that the ranking and ranking conditional on model predictions fit simulation expectations well.

```{r}
binnedplot(predict(l_glm_1976), resid(l_glm_1976))
```
```{r warning = F, message = F}
library(DHARMa)
plot(simulateResiduals(l_glm_1976))
```
 


### 3(b) model interpretation


```{r}
oldwidth <- options()$width
options(width=200)
summary(l_glm_1976)
options(width=oldwidth)
```

To predict if a voter with vote for a Republican (or Democrat) we see some interesting trends. First, holding all other features equal, if an individual is **black** (not white), we should expect $e^{\beta_\text{black}}$ (`r exp(coef(summary(l_glm_1976))["black","Estimate"])`) multiplicative decrease in the odds ratio of voting for a Republican over a Democrat. Similarly, but in reverse, **not having any household member being in a union** (instead of the reverse) is observed as increasing the odds ratio of voting for a Republican  over a Democrat by `r exp(coef(summary(l_glm_1976))["union2. no, no one in the household belongs t","Estimate"])`. Not surprisingly, given that the current president is a Republican, if **one disapproves of the current president** (as opposed to approving of the current president), would lead to a multiplicative decrease in the odds ratio of voting for a Republican  over a Democrat by multiplication by `r exp(coef(summary(l_glm_1976))["presapprov2. disapprove","Estimate"])`. 

There are also some non-binary variable that were significant, but they require more nuanced interpretations. For different **ideological** views we see impacts on whether the individual will vote Republican  or Democrat. In this case, the base case is “extremely liberal” and interestedly we see that if one is moves from “extremely liberal to “liberal or “slightly liberal” we’d expect a multiplicative decrease in the odds ratio of voting for a Republican  over a Democrat by `r exp(coef(summary(l_glm_1976))["ideo72. liberal","Estimate"])` or `r exp(coef(summary(l_glm_1976))["ideo73. slightly liberal","Estimate"])` respectively. We don’t observe significant differences (holding all other variables equal) between the odds ratio between “extremely liberal” vs any time of conservative. In context this is surprising, but one could hypothesis that this may be related to the fact that certain differences between “extremely liberal” and “conservatives” are already captured in other variables. Another significant categorical variable is **maximum education level attained**, where the base case is only some grade school (0-8 grades). Compared with this base level we do see that as we increase education levels the estimated coefficient for voting Republican  over Democrat increases. We can’t directly use the p-values of the higher education levels to test if a high school education vs college education significantly impacts a voter’s likelihood of voting for a Republican, but looking at the difference of the coefficients and the standard error we might guess (without directly doing the test) that they are not given the standard error for each $\beta$ value is larger than the difference between the $\beta$ values.



\rm


## Problem 3(c)

Make a new data frame by converting all of the variables in the 1976 data frame to numeric variables, using \texttt{as.numeric()}.  Repeat part (b), using this new data frame.

\em

```{r}
nes_1976_numeric <- nes_1976 %>% sapply(as.numeric) %>% data.frame
l_glm_1976_numeric <- glm(rep_pres_intent ~ ideo7 + black + female + presapprov + age + 
                educ1 + urban + income + union + perfin1, 
              data = nes_1976_numeric, family=binomial)
```
### 3(c) Fit assessment

The fit of the above model is worse (in terms of shrinking the residual deviance). From the marginal plot, this may be related to the less optimal fit relative to ideology and age. Other than that - it does look like a similar fit to the factor based model.

```{r}
binnedplot(predict(l_glm_1976_numeric), resid(l_glm_1976_numeric))

plot(simulateResiduals(l_glm_1976_numeric))

mmplot(l_glm_1976_numeric,exclude = c(2,3,4,7,9,10))
```


### 3(c) model interpretation

This numerical model suggests similar parameters are useful. A similar interpretation of **black** (vs white), **not having any household member being in a union** (instead of the reverse) and presidental approval exists for this model as the last. Specifically, holding all other features equal, if an individual is **black** (not white), we should expect $e^{\beta_\text{black}}$ (`r exp(coef(summary(l_glm_1976_numeric))["black","Estimate"])`) multiplicative decrease in the odds ratio of voting for a Republican  over a Democrat, **not having any household member being in a union** (instead of the reverse) is observed as increasing the odds ratio of voting for a Republican over a Democrat by a multiplication of `r exp(coef(summary(l_glm_1976_numeric))["union","Estimate"])`, and if **one disapproves of the current president** (as opposed to approving of the current president), would lead to a decrease in the odds ratio of voting for a Republican over a Democrat by a multiplication of`r exp(coef(summary(l_glm_1976_numeric))["presapprov","Estimate"])`. 

We now can only fit a linear trend (that ranges from "extremely liberal" to "extremely conservative"), but it is significant, a one "unit" increase (between the levels) of ideology towards conservatism, we should expect $e^{\beta_\text{ideo7}}$ (`r exp(coef(summary(l_glm_1976_numeric))["ideo7","Estimate"])`) multiplicative decrease in the odds ratio of voting for a Republican over a Democrat. Education attainment is no longer a significant variable.


```{r}
summary(l_glm_1976_numeric)
```



\rm


## Problem 3(d)

Compare the results from parts (b) and (c}).  How are they consistent? Inconsistent?  In this case, would you recommend converting all of the variables to numeric?

\em

Both models contain similar information, and, in terms of significant variables, all variables that are significant in the first model (except education) are significant in the second model as well. Naturally, the "linearizing" of the ordered factor variables reduced the interpretability of the model. Additionally, although not significant, the interpretation of the linearization of the levels of income is probably a worse assumption. Requiring any of these ordered factors to be linear requires (1) a strictly increasing or decreasing trend, and (2) that the change between levels should be the same. In the interpretation of the first model we saw (1) violated with ideological levels and (2) violated with education levels. As such I wouldn't recommend converting all the variables to numeric in this case.



\rm


## Problem 3(e)

NOT TO TURN IN.  If you have time, it is interesting to repeat parts (a) through (d) for the presidential election year 2000.  In this year, the democratic candidate was Al Gore, who served as vice-president to Bill Clinton, who was impeached but not removed from office.  The Republican candidate was George Bush.  How are the results from 1976 and 2000 similar?  How are they different?

\em

```{r}
nes_2000 <- nes_all %>% filter(year == "2000") %>%
  dplyr::select(ideo7, black, female,
                          rep_pres_intent,
                          presapprov, age,
                          educ1, urban, income, union,
                          perfin1)
```

```{r}
l_glm_2000 <- glm(rep_pres_intent ~ ., data = nes_2000, family = "binomial")
```

### 3(e) Fitted assessment

Interestingly, this model (relative to the model for 1976) seems a bit better fitted. Most of the binned residuals still fit in the expected variablity range, with much fewer bins (still at the extremes) that are under-estimated in the extremes. Some of this may actually just be related to the fact that the 2000 sample has `r nrow(nes_2000)` individuals whereas the 1976 sample had more (`r nrow(nes_1976)`).

```{r}
binnedplot(predict(l_glm_2000), resid(l_glm_2000))
```

```{r warning = F, message = F}
plot(simulateResiduals(l_glm_2000))
```

### 3(e) model interpretation

Some similarities exist between the two models (in terms of which variables are important), yet we see new trends emerging and some naturally being reverse based on intuition. 

For binary variables, **black** and **previous president disapproval** still is significant (though the impact of previous president disapproval is not the reverse - favoring a vote for a Republican over a Democrat)^[This makes sense become the previous president is now a Democrat, while in 1976 it was a Republican.]. Specifically, holding all other features equal, if an individual is **black** (not white), we should expect $e^{\beta_\text{black}}$ (`r exp(coef(summary(l_glm_2000))["black","Estimate"])`) multiplicative decrease in the odds ratio of voting for a Republican over a Democrat, and if **one disapproves of the current president** (as opposed to approving of the current president), would lead to an increase in the odds ratio of voting for a Republican over a Democrat by a multiplication of `r exp(coef(summary(l_glm_2000))["presapprov2. disapprove","Estimate"])`. 

We see that ideology and education level is not longer a significant variable (and unlikely and comparison due to the very high standard errors estimated for each $\beta$ value). On the other hand, **suburban** (vs urban) looks to cause a significant to a `r exp(coef(summary(l_glm_2000))["urban2. suburban areas","Estimate"])` multiplicative decline in the odds ratio of voting for a Republican over a Democrat. We also find the age is now a significant variable (with a 1 year increase in age) related to a `r exp(coef(summary(l_glm_2000))["age","Estimate"])` multiplicative decrease in the odds ratio of voting for a Republican over a Democrat.

```{r}
oldwidth=options()$width
options(width=200)
summary(l_glm_2000)
options(width=oldwidth)
```

### 3(e) numeric

```{r}
nes_2000_numeric <- nes_2000 %>% sapply(as.numeric) %>% data.frame
l_glm_2000_numeric <- glm(rep_pres_intent ~ ideo7 + black + female + presapprov + age + 
                educ1 + urban + income + union + perfin1, 
              data = nes_2000_numeric, family=binomial)
```

#### 3(e) numeric, fitted assessment

Again, this model looks a bit worse than the non-linear transformation of factor variables. We see at differences in education levels are not well captured, and maybe a bit suboptial fit in the linearized ideological variable (via the marginal plot). Interestly, the binned residuals look like a better fit (with only the most extreme bin for Republican voting to be highly underestimated). 

```{r}
binnedplot(predict(l_glm_2000_numeric), resid(l_glm_2000_numeric))

plot(simulateResiduals(l_glm_2000_numeric))

mmplot(l_glm_2000_numeric,exclude = c(2,3,4,7,9,10))
```

#### 3(e) numeric, model interpretation

The TA leaves this up to the student to interpret, but we note that ideology of the individual is now significant, and urban is no longer significant. Additionally similar comments as in (d) relate to the linearized and factor based models.

```{r}
summary(l_glm_2000_numeric)
```

### 3(e) Overall...

This exercise suggests that different voter characteristics affected voter decisions across the years (although some like previous president approval^[With the party in power flip.] and race seem pretty constant). Additionally, ideology seems to have some impact (but maybe not linear in all cases). Interestingly, the model fit on the 2000 data seems to better conform with the linear model assumptions. The differences between the models fit on the different years relate to the impact of education level (significant for 1976 but not 2000) and urban vs suburban (significant for 2000 but not 1976). As captured in the discussion about interested 1976 trend for ideology, we should make sure **not** to claim causal structure related for voter preferences.

\rm


