---
title: "Life Expectancy Example"
author: "Jeffrey B. Arnold"
date: "05/12/2015"
output: html_document
---

```{r message=FALSE}
library("car")
library("dplyr")
library("broom")
library("boot")
```


This example will use data from the Barro and Lee dataset to analyze 
life expectancy. 
This is a pedagogical example, and does not represent a sophisticated epidemiological model.
Data are 138 countries in 1985.

```{r barro_csv}
barro_raw <- read.csv("http://pols503.github.io/pols_503_sp15/data/barro.csv")
```

Subset of barro data with only the variables we need, and dropping missing observations.
```{r barro}
barro <- barro_raw %>%
  select(lifeexp, school, gdpcap85,
         civlib5, wartime) %>%
  na.omit()

glimpse(barro)
```

gdpcap85

:   Per capita GDP in 1985, thousands of international dollars

school

:   Average years of education

civlib5

:   low = 1 to high = 7 scale of civil liberties

wartime

:   Percent of recent history spent at war



We will compare two models of life expectancy.

Model 1 regresses life expectancy on `gdpcap85`, `school`, `civlib5`, and `wartime`.
```{r mod1}
mod1 <- lm(lifeexp ~ gdpcap85 + school + civlib5 + wartime,
           data = barro)
mod1
```
Model 2 regresses life expectancy on the logarithm of `gdpcap85`, the logarithm of `school`, `civlib5`, and `wartime`.
```{r}
mod2 <- lm(lifeexp ~ log(gdpcap85) + log(school) + civlib5 + wartime,
           data = barro)
mod2
```

## Comparing Models

We can compare the models from the previous section using $R^2$, adjusted-$R^2$, the standard error of the regression $\hat\sigma^2_{\epsilon}$, AIC, and BIC. 

The summary of the lm object shows the R squared values and the standard deviation of the regression ("Residual standard error"):
```{r}
mod1_summary <- summary(mod1)
mod1_summary
```
These can be extracted from the object returned by `summary`.
Explore the object to find the elements corresponding to those values:
```{r results = 'hide'}
str(mod1_summary)
```
The functions `AIC` and `BIC` do what you would expect them to do,
```{r}
AIC(mod1)
BIC(mod1)
```
The `glance` function from the **broom** package returns all of these within a model frame:
```{r mod1_glance}
glance(mod1)
```

Now let's calculate these values for Model 2:
```{r}
glance(mod2)
```

By any of those measures, Model 2 is the better fitting model.

## Cross Validation

The **boot** package contains functions to do cross validation. 

Here is an example of cross-validation using the previous models.
Although `glm` is used instead of `lm`, it produces the same results as `lm`.
This is needed for the `cv.glm` function to work:
```{r}
mod1_glm <- glm(lifeexp ~ gdpcap85 + school + civlib5 + wartime, data = barro)
```
This runs 10-fold cross validation
```{r}
cv_err_K10 <- cv.glm(barro, mod1_glm, K=10)
```

The object returned by `cv.glm` includes the original function call, the value of `K`, the predication error `delta`, and the random number seed.
For now, all we care about is the prediction error.
By default, `cv.glm` returns the mean squared error. There are two elements that it returns: the first is the average of the mean squared errors from the folds; the second corrects for bias occuring from using k-folds rather than leave-one-out cross-validation.
```{r}
cv_err_K10$delta
```
5-fold cross-validation
```{r}
cv_err_K5 <- cv.glm(barro, mod1_glm, K=5)
cv_err_K5$delta
```
And leave-one-out cross-validation (no `K` argument is given):
```{r}
mod_cv_loo <- cv.glm(barro, mod1_glm)
mod_cv_loo$delta
```

The interpretation in these cases is that the average prediction RMSE is over 30 years. 
This is not good. Not good at all.

We can calculate the same for Model 2:
```{r}
mod2_glm <- glm(lifeexp ~ log(gdpcap85) + log(school) + civlib5 + wartime, data = barro)

cv.glm(barro, mod2_glm, K=10)$delta
cv.glm(barro, mod2_glm, K=5)$delta
cv.glm(barro, mod2_glm)$delta
```
In these cases the prediction RMSE is around 14 years. 
This is okay.
Maybe not as good as we would like, but a lot better than missing by 30 years.

We could also try adding a squared term to war to see whether that improves the model.
We'll use a 10-fold cross-validation for it.
```{r}
mod3 <- glm(lifeexp ~ log(gdpcap85) + log(school) + civlib5 + wartime + I(wartime ^ 2), data = barro)
cv.glm(barro, mod3, K = 10)$delta
```

It shows a small improvement in the prediction RMSE. 
Adding squared `wartime` seems to improve the model, but it is not nearly as important as getting the functional form of `gdpcap85` and `school` correct.


## Predictive Comparisons by Simulation and Bootstrap

### Simulating Coefficients from Asymptotic Normal

The first step is to draw values of $\beta$ from a multvariate normal distribution centered at $\hat\beta$ with covariance matrix $vcov(\hat\beta)$:
$$
\tilde\beta \sim N(\hat\beta, V(\hat\beta))
$$
This relies on the CLT result that $\hat\beta$ approx multivariate normal as $n \to \infty$.
We will randomly sample simulations of $\beta$ from that distribution:
```{r}
n <- 1024
simbetas <- MASS::mvrnorm(n, coef(mod2), vcov(mod2))
```

For `log(gdpcap85)` compare it at its mean versus 1 standard deviation
```{r}
barro_low <- summarize(na.omit(barro),
                      gdpcap85 = exp(mean(log(gdpcap85))),
                      school = exp(mean(log(school))),
                      civlib5 = mean(civlib5),
                      wartime = mean(wartime))

barro_high <- summarize(na.omit(barro),
                       gdpcap85 = exp(mean(log(gdpcap85)) + sd(log(gdpcap85))),
                       school = exp(mean(log(school))),
                       civlib5 = mean(civlib5),
                       wartime = mean(wartime))
```

Simulate and calculate a standard error around the difference
```{r}
xlow <- model.matrix(~ log(gdpcap85) + log(school)
                      + civlib5 + wartime, data = barro_low)
xhigh <- model.matrix(~ log(gdpcap85) + log(school)
                      + civlib5 + wartime, data = barro_high)
diffs <- rep(NA, nrow(simbetas))
for (i in 1:nrow(simbetas)) {
  betas <- simbetas[i, ]
  diffs[i] <- xhigh %*% betas - xlow %*% betas
}

```

Two packages implement this method of simulatin for comparison:

- [Zelig](http://zeligproject.org/)
- [simcf](https://github.com/chrisadolph/tile-simcf): Chris Adolph's package written while teaching this course.

## Bootstrapping

The function `Boot` from **car** can be used to bootstrap coefficients.
```{r}
Boot(mod2)
```

The **boot** package has more powerful and general functions.

Bootstrapping can be done manually with `dplyr` functions:
```{r}
simulations <- list()
for (i in 1:1024) {
  # Randomly resample data
  .data <- sample_frac(barro, replace = TRUE)
  # Run a new regression
  mod <- lm(lifeexp ~ log(gdpcap85) + log(school) + civlib5 + wartime, data = .data)
  # Predict
  lower <- predict(mod, newdata = barro_low)
  higher <- predict(mod, newdata = barro_high)
  simulations[[i]] <-
    data.frame(diff = higher - lower)
}
sims <- bind_rows(simulations)

```
From these simulations, one way to calculate a 95% cofidence interval is just the 2.5--97.5 percentiles:
```{r}
quantile(sims$diff, c(0.025, 0.975))
```

This code could have been streamlined using the `bootstrap` function in the **broom** package.

* * * 

Example adapted from Christopher Adolph (Spring 2014) "Linear Regression: Specification and Fitting" [Lecture slides]. <http://faculty.washington.edu/cadolph/503/topic5.pw.pdf>.
