---
title: "Homework 07 Solutions"
date: "10/30/2020"
output: 
  pdf_document:
    keep_tex: true
fig_caption: true
header-includes:
  - \usepackage{amsmath}
  - \usepackage{amsthm}
  - \DeclareMathOperator{\Cov}{Cov}
  - \DeclareMathOperator{\Corr}{Corr}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(cache = FALSE, echo = FALSE, include = TRUE, message = FALSE)
options(scipen = 1, digits = 2)

# library(car)
library(GGally)
# library(knitr)
library(marginalmodelplots)
library(pander)
library(tidyverse)
library(arm)
library(DHARMa)
```

# Sheather 8.3.2
First of all, note that the outcome values are not independent, because exactly 10 jurisdictions produce a top 10 finalist in a given year. Since the data ranges over 9 years, that means that the sum of the `Top10` values must be $9\times 10 = 90$. That's actually a pretty strong constraint; it means, for example, that if we observe 10 jurisdictions with the value 9 for `Top10`, we automatically know the value for the other 41 jurisdictions (i.e., 0). Of course, we never expect regression assumptions to hold exactly in practice, but I'd take this model with a big grain of salt since the samples are so heavily dependent.

That said, let's go ahead and visualize the data (Figure \ref{f:pairs1}).

Three predictors, `LogPopulation`, `LogContestants`, and `Latitude` have high marginal correlations with `Top10`. `LogTotalArea` and `Longitude` have pretty small correlations, but of course these are marginal correlations, so they could still be strongly related to the outcome.

There are some unsurprising correlations between some pairs of predictors. `LogPopulation` and `LogContestants` are positively correlated, which makes sense: larger populations produce more contestants. Longitude is positively correlated with `LogTotalArea`: western states are larger on average than eastern ones.

The predictors have reasonably symmetric and unimodal distributions, presumably thanks in part to the log transformations.

```{r scatter, fig.cap="Pairs plot for Miss America data, excluding abbreviation column.\\label{f:pairs1}"}
dat <- read_table2("MissAmerica2008-binomial.txt")
dat_original <- dat
# attach(dat)
dat %>% dplyr::select(-abbreviation) %>% ggpairs
```

## (a): Full model
Here's the full model, with a summary below. Standardized deviance residuals are plotted against each of the predictors in Figure \ref{f:resid1}. Marginal model plots are in Figure \ref{f:marginal1}.

The residuals look essentially patternless, and the nonparametric curves for the fitted and the observed values in the marginal model plots look roughly the same. Since `Longitude` is not significant, and the problem requests a model in which all the predictors are significant, let's drop `Longitude` and refit (see below).

```{r model1}
attach(dat)
m1 <- glm(cbind(Top10, 9 - Top10) ~ LogPopulation + LogContestants + 
            LogTotalArea + Latitude + Longitude, family = binomial)

pander(summary(m1), caption = "")
```

```{r residuals1, fig.cap="Standardized deviance residuals (left) against fitted values for the full model.\\label{f:resid1}"}
H <- influence(m1)$hat  # Hat matrix
resid_deviance_std <- residuals(m1)/sqrt(1 - H)

attach(dat)
par(mfrow = c(2, 3),  mar = c(4, 4.5, 2, 1))
plot(LogPopulation, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(LogContestants, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(LogTotalArea, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(Latitude, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(Longitude, resid_deviance_std, ylab = "Standardized deviance residuals")
```

```{r marginal1, fig.cap="Marginal model plots for the full model\\label{f:marginal1}"}
mmplot(m1)
```

## (a): Reduced model without `Longitude`

```{r model2}
attach(dat)
m2 <- glm(cbind(Top10, 9 - Top10) ~ LogPopulation + LogContestants + 
            LogTotalArea + Latitude, family = binomial)
```

The model summary is below. Residuals are plotted in Figure \ref{f:resid1} and the marginal model plots are in Figure \ref{f:marginal2}.

Once again, the residuals look essentially patternless, and the curves in each marginal model plot follow each other closely. All the predictors are now statistically significant. The AIC for this model (`r m2$aic`) is also lower than for the full model (`r m1$aic`). This model seems reasonable.

```{r summary2}
pander(summary(m2), caption = "")
```

```{r residuals2, fig.cap="Standardized deviance residuals (left) against fitted values for the reduced model.\\label{f:resid2}"}
H <- influence(m2)$hat  # Hat matrix
resid_deviance_std <- residuals(m2)/sqrt(1 - H)

attach(dat)
par(mfrow = c(2, 3),  mar = c(4, 4.5, 2, 1))
plot(LogPopulation, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(LogContestants, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(LogTotalArea, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(Latitude, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(Longitude, resid_deviance_std, ylab = "Standardized deviance residuals")
```

```{r marginal2, fig.cap="Marginal model plots for the reduced model\\label{f:marginal2}"}
mmplot(m2)
```

## (b)
Recall that a "bad" leverage point was defined as a leverage point that is also an outlier. To detect this, let's examine the diagnostic plot that R generates with standardized Pearson residuals on the y-axis and leverage on the x-axis, in Figure \ref{f:leverage}.

Point 12 has both high leverage and a high standardized residual value, so it could be considered a bad leverage point. Let's refit the model without that point. 

Once again, the residual plots (Figure \ref{f:resid3}) and marginal model plots (Figure \ref{f:marginal3}) look good. From the summary below, we see that the coefficient for `Latitude` is no longer significant. Point 12 corresponded to Hawaii, which has a very low latitude, so it makes sense that that point was  responsible for making `Latitude` statistically significant. The model without Hawaii is probably a better model of the continental US. (Alaska is still in the picture, but evidently it's not a bad leverage point, since it didn't show up as such on the plot.)

```{r leverage, fig.cap="Residuals vs. leverage for the reduced model\\label{f:leverage}"}
plot(m2, which = 5)  # Residuals vs. leverage
```

```{r model3}
dat <- dat[-12, ]
attach(dat)
m3 <- glm(cbind(Top10, 9 - Top10) ~ LogPopulation + LogContestants + 
            LogTotalArea + Latitude, family = binomial)

pander(summary(m3), caption = "")
```

```{r residuals3, fig.cap="Standardized deviance residuals (left) against fitted values for the reduced model without point 12.\\label{f:resid3}"}
H <- influence(m3)$hat  # Hat matrix
resid_deviance_std <- residuals(m3)/sqrt(1 - H)

attach(dat)
par(mfrow = c(2, 3),  mar = c(4, 4.5, 2, 1))
plot(LogPopulation, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(LogContestants, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(LogTotalArea, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(Latitude, resid_deviance_std, ylab = "Standardized deviance residuals")
plot(Longitude, resid_deviance_std, ylab = "Standardized deviance residuals")
```

```{r marginal3, fig.cap="Marginal model plots for the reduced model without point 12.\\label{f:marginal3}"}
mmplot(m3)
```

As described above, a "bad" point is related to high residuals and high leverage (making it influential). In order to assess this, we need the residuals to meet our expectations (leverage is only related to the dependent variables, not the model). In Figure \ref{f:binned2} we can observe that, only once does a binned grouping have a average residual beyond the 2 standard deviations, this suggests that the relationship between residuals and fitted values isn't too concerning. Moreover, we can't observe any nonconstant pattern with the residual averages. Additionally, from the marginal plots we know that only 2 observations have expected values below -4, which might give us pause to overinterpret the plot. Under the model assumptions, we can simulate potential outcomes and compare these to the actual values. Figure \ref{f:sim2} suggests that the residuals observed compared to those from the simulations do not suggest significant differences, specifically the ranks for the residuals vs what ranking would be expected seem to have similar distributions. 

```{r binned-resid, fig.cap="Binned residual plot for second model.\\label{f:binned2}"}
binnedplot(x = predict(m2), y = resid(m2))
```

```{r  sim-assessment, fig.cap="Simulation based assessment of residuals.\\label{f:sim2}"}
simulateResiduals(m2,n =10000, seed = 1) %>% plot()
```




\pagebreak

## (c)
Let's interpet the coefficients from the model without Hawaii.

Mathematically, the intercept means that when the covariates are all 0, the estimated log odds of producing a top 10 finalist is -9.741, which corresponds to a probability of 0.000059. Of course, a jurisdiction with 0 for all covariate values is not realistic or of interest in the model.

Note that the odds are
\begin{align*}
  \frac{\theta(X)}{1 - \theta(X)}
\end{align*}
where $\theta(X)$ represents the probability of producing a top 10 finalist (given covariates $X$) for a single Bernoulli trial. Since we're pretending that $Y_i \sim Bin(10, \theta(X_i))$, we have that $\theta(X)$ corresponds to the probability of producing a top 10 finalist in a single year. It does not represent the probability of producing at least one top 10 finalist over 10 ten years or anything like that.

The coefficient on `LogPopulation` indicates that, if we observe two jurisdictions that have the same covariate values except that they differ by 1 in `LogPopulation`, then the log odds of producing a top 10 finalist for the more populous jurisdiction are 0.6171 higher than for the less populous jurisdiction. Equivalently, the odds are `exp(0.671) = 1.9` higher.

The interpretation is analogous for the other coefficients.

# Problem 2

```{r bernoulli-data, message = F}
dat2 <- read.csv("MissAmerica2008-bernoulli.txt",sep = " ") 
dat2 <- dat2[-10*(1:51),] ## delete these rows
## and then continue with the problem using this modified data frame.
```

##  (a)

It appears that the "bernoulli" dataset is just yearly for each state if their contestant was in the top 10. This means that the "binimial" dataset was a summary of the "bernoulli" for each state. 

```{r message = F}
testthat::expect_equal(nrow(dat2), 9 * nrow(dat_original))

info_per_state <- dat2 %>% group_by(abbreviation) %>%
  summarize(t = n(),  count = sum(Top10)) 
merger <- dat_original %>% full_join(info_per_state, by ="abbreviation")

testthat::expect_equal(merger$Top10, merger$count)

```

## (b)

Using this Bernoulli data, we show the model without `Longitude` below. Because both this model and the second model in part 1a are trying to to model $\theta[X_i]$ - the probability in being in the top 10 contestants related to state information and they both use the same information we expect the model's parameters to be very similar. 

```{r model2-revisited}
m2_revisited <- glm(Top10 ~ LogPopulation + LogContestants + 
            LogTotalArea + Latitude, family = binomial,
            data = dat2)

pander(summary(m2_revisited), caption = "")
```

# (c)


```{r}
binomial_m2_aic <- AIC(m2) 
bernoulli_m2_aic <- AIC(m2_revisited) 
```

We should not expect the AICs of the two models to be the same, nor are they (with the Bernoulli `m2` model having an AIC of `r bernoulli_m2_aic` and the binomial `m2` model having an AIC of `r binomial_m2_aic`). This directly related to the definition of AIC, which is 

\[-2 \cdot \text{Log likelihood} + 2 \cdot p\;,
\]
where $p$ is the number of parameters of the model (5 in our models). 

The first model has an assumption of the observations being $\text{Binomial}(9, \theta[X_i])$ and the second uses the assumption that the observations are $\text{Bernoulli}(\theta[X_i])$. The likelihoods for a state therefore look slightly different, specifically:

\[
\begin{pmatrix}9 \\ x_i \end{pmatrix} \theta[X_i]^{y_i}(1-\theta[X_i])^{9-y_i}
\]
for the Binomial assumption and 
\begin{align*} 
& \prod_{z=2000}^{2009}  \theta[X_i]^{y_{iz}}(1-\theta[X_i])^{1-y_{iz}} \\
\text{or } & \theta[X_i]^{y_{i}}(1-\theta[X_i])^{9-y_{i}}
\;,
\end{align*}
 for the Bernoulli assumption^[This is because we can see $y_i = \sum_{z=2000}^{2009}  y_{iz}$.]. As the log likelihood is over all $i$, this means that the Log likelihoods for the model from part 1 (the Binomial assumption) has an additional
 
 \[
 \sum_{i=1}^{51} \log\left(\begin{pmatrix}9 \\ y_i \end{pmatrix}\right)
 \]
 term. THe difference between the two AICs should then be by $2 * \sum_{i=1}^{51} \log\left(\begin{pmatrix}9 \\ y_i \end{pmatrix}\right)$ (from the definition of the AIC). The final 2 lines of code below show that to be true.
 

```{r, echo = T}
binomial_m2_aic - bernoulli_m2_aic  
```

```{r, echo = T}
sapply(dat_original$Top10, function(x) choose(9, x)) %>%
  log() %>% sum() %>% "*"(.,-2)
```

