04-Resampling.Rmd

# Resampling {#resampling}

```{r ch-4-setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r ch-4-packages, include = FALSE}
library(rsample)
library(parsnip)
library(yardstick)
library(workflows)
library(tune)
library(ggplot2)
library(dplyr)
library(broom)
library(ISLR)
library(mvtnorm)
library(rsample)
library(microbenchmark)
library(furrr)
library(doFuture)
theme_set(theme_bw())
```

## Training vs Test Data

We will explore estimating test error rates using training/testing data sets
generated by functionality in the `rsample` package. We will illustrate these 
concepts on the `Auto` dataset that is located in the `ISLR` package. The data
consist of 392 observations on nine attributes with the dependent variable being
miles per gallon, or mpg. Note that the `origin` variable should be considered
a factor, or character, variable. 

```{r}
df <- Auto %>% 
  dplyr::mutate(., origin = as.factor(origin))
```

```{r}
table(df$origin)/392
```

Why do we set a seed in the following code chunk?
```{r}
set.seed(298329)
auto_split <- rsample::initial_split(df, prop = .5)
train_data <- rsample::training(auto_split)
test_data <- rsample::testing(auto_split)
```

```{r}
table(train_data$origin)/196
```

This looks pretty even. If we want to ensure the correct proportions within
all of the strata, we can use the `strata` argument in `initial_split`.

```{r}
set.seed(298329)
auto_split <- rsample::initial_split(df, prop = .5, strata = origin)
train_data <- rsample::training(auto_split)
test_data <- rsample::testing(auto_split)
```

```{r}
table(train_data$origin)/196
```

## LM Models

What does a plot of horsepower vs. mpg look like on the training set? Add an lm 
fit, a polynomial fit, and others. 

### Linear Model

```{r}
p <- ggplot(data = train_data,
            aes(x = horsepower, y = mpg))
p + geom_point(col = "#6E0000") +
  stat_smooth(method = "lm", col = "black")
```

```{r}
# p <- ggplot(data = train_data,
#             aes(x = horsepower, y = mpg))
p + geom_point(col = "#6E0000") +
  stat_smooth(method = "lm", formula = y ~ poly(x, 15), col = "black")
```

### Test Error Rates

How can we estimate the test error rate in this problem? First, let's create and
then fit the models. 

#### Linear Model:

Declare the model first. 

```{r}
lm_model <- parsnip::linear_reg() %>% ## Class of problem
  parsnip::set_engine("lm") %>% ## The particular function that we use
  parsnip::set_mode("regression")
lm_model
```

Note that we will fit the model separately in this case because we will fit
several different types of linear models. The original model is sufficient for
fitting all of them. 

```{r}
lm_fit_1 <- lm_model %>% 
  parsnip::fit(mpg ~ horsepower, data = train_data)
broom::tidy(lm_fit_1)
```

Here is the test root mean square error (RMSE) and other metrics. 

```{r}
lm_fit_1 %>% 
  predict(test_data) %>% 
  bind_cols(test_data) %>% 
  metrics(truth = mpg, estimate = .pred)
```

#### Quadratic Model:

```{r}
lm_model %>% 
  parsnip::fit(mpg ~ poly(horsepower, 2), data = train_data) %>% 
  predict(test_data) %>% 
  bind_cols(test_data) %>% 
  metrics(truth = mpg, estimate = .pred)
```

#### Cubic Model:
```{r}
lm_model %>% 
  parsnip::fit(mpg ~ poly(horsepower, 3), data = train_data) %>% 
  predict(test_data) %>% 
  bind_cols(test_data) %>% 
  metrics(truth = mpg, estimate = .pred)
```

It looks like a quadratic fit is better than a linear fit, and a cubic
may not add much. However, should we base the results off of a single 
fit across three different models based on a single test/train set.

## Variability

You can do this for multiple splits (test/train) and you'll notice different
values for RMSE (and RSQ, MAE), i.e. there is quite a bit of variability in
the estimates. Look at Figure 5.2 in the textbook. 

Let's look at a single split and plot the results of the model fits on 
the MSE. 

```{r}
results <- tibble(d = 1:10,
                  RMSE = 1:10)

for(j in 1:10){
  results$RMSE[j] <- lm_model %>% 
    parsnip::fit(mpg ~ poly(horsepower, get("j")), data = train_data) %>% 
    predict(test_data) %>% 
    bind_cols(test_data) %>% 
    metrics(truth = mpg, estimate = .pred) %>% 
    dplyr::filter(., .metric == "rmse") %>% 
    pull(.estimate)
}
```

```{r}
p <- ggplot(data = results,
            aes(x = d, y = RMSE^2))
p + geom_line() +
  geom_point() +
  scale_x_continuous("Degree of Polynomial") +
  scale_y_continuous("Mean Squared Error", limits = c(15, 29), 
                     breaks = seq(16, 28, by = 2)) +
  theme_bw()
```

Now let's look at ten different splits and notice the variability. 

```{r}
set.seed(584)
results <- tibble(d = rep(1:10, times = 10),
                  RMSE = 1:100,
                  split = rep(1:10, each = 10))

for(i in 1:10){
  new_split <- initial_split(df, prop = .5, strata = origin)
  train_data <- rsample::training(new_split)
  test_data <- rsample::testing(new_split)
  for(j in 1:10){
    results$RMSE[(i - 1)*10 + j] <- lm_model %>% 
      parsnip::fit(mpg ~ poly(horsepower, get("j")), data = train_data) %>% 
      predict(test_data) %>% 
      bind_cols(test_data) %>% 
      metrics(truth = mpg, estimate = .pred) %>% 
      dplyr::filter(., .metric == "rmse") %>% 
      pull(.estimate)
  }
}
```

```{r}
p <- ggplot(data = results,
            aes(x = d, y = RMSE^2, col = as.factor(split)))
p + geom_line() +
  geom_point() +
  scale_x_continuous("Degree of Polynomial") +
  scale_y_continuous("Mean Squared Error", limits = c(14, 29), 
                     breaks = seq(14, 28, by = 2)) +
  scale_color_discrete("split") +
  theme_bw()
```


### Cross Validation

A better approach is use what's called cross validation. This is where we split 
the original data into a testing/training set, but then we further split our 
training set into $k$ new training/validation sets. Some people will refer to
the new sets as "analysis/assessment" sets. 

```{r}
set.seed(130498)
auto_split <- initial_split(df, prop = .5, strata = origin)
train_data <- rsample::training(new_split)
test_data <- rsample::testing(new_split)

folds <- vfold_cv(train_data, v = 5, strata = origin)
folds
```

Let's fit the simple linear regression model on the cross-validation splits. To
do this, we need to introduce the idea of a "workflow" from the `workflows()`
package. 

```{r}
## lm_model is defined above
lm_wf_1 <- workflow() %>% 
  add_model(lm_model) %>% 
  add_formula(mpg ~ poly(horsepower, 3))

lm_fit_1_rs <- 
  lm_wf_1 %>% 
  tune::fit_resamples(folds)
lm_fit_1_rs
collect_metrics(lm_fit_1_rs)
```

Let's do this for each of our ten polynomial models

```{r}
auto_recipe <- recipe(mpg ~ horsepower, train_data) %>% 
  step_poly(horsepower, degree = tune()) 
auto_recipe

tg <- tune_grid(lm_model, auto_recipe, resamples = folds, 
                grid = expand.grid(degree = 1:10))
collect_metrics(tg) %>% filter(.metric == "rmse")
p <- ggplot(data = collect_metrics(tg) %>% filter(.metric == "rmse"),
            aes(x = degree, y = mean^2))
p + geom_line()
```

How about doing this ten times? 

```{r, cache = T}
set.seed(9823)
for(i in 1:10){
  folds <- vfold_cv(train_data, v = 5, strata = origin)
  auto_recipe <- recipe(mpg ~ horsepower, train_data) %>% 
    step_poly(horsepower, degree = tune()) 
  tg <- tune_grid(lm_model, auto_recipe, resamples = folds, 
                  grid = expand.grid(degree = 1:10))
  tmp <- collect_metrics(tg) %>% 
    mutate(split = i)
  if(exists("cv_results")){
    cv_results <- bind_rows(cv_results, tmp)
  } else{
    cv_results <- tmp
  }
}
```

```{r}
p <- ggplot(data = cv_results %>% filter(.metric == "rmse"),
            aes(x = degree, y = mean^2, col = as.factor(split)))
p + geom_line() +
  geom_point() +
  scale_x_continuous("Degree of Polynomial", breaks = 1:10) +
  scale_y_continuous("Mean Squared Error", limits = c(14, 29), 
                     breaks = seq(14, 28, by = 2)) +
  scale_color_discrete("split") +
  theme_bw()
```

### Bootstrapping

- The *bootstrap* is a resampling tool (methodology?) that is used to estimate
or quantify the uncertainty associated with a given statistic or estimator.
- The term is derived from the phrase *"to pull oneself up by one's bootstraps."*
- Examples:
  - Suppose you fit a simple linear regression and compute $\hat{\beta_0}$ and
  $\hat{\beta_1}$ and you do not know how to compute the usual t-test, CI, etc.
  You can resample the data, re-fit regressions, and essentially recreate a 
  sampling distribution of those statistics. 
  - Let $X$ and $Y$ be the returns on two financial assets and $\alpha$ be the 
  percentage of our money invested in $X$ and $1-\alpha$ in $Y$. We wish to 
  choose $\alpha$ that minimizes the total risk, or variance, of the investment.
- **Warning**: The bootstrap is not an oracle. You don't necessarily get 
something for free. 

### Example 

Suppose that we wish to invest a fixed sum of money in two financial assets that
yield returns of $X$ and $Y$ (random variables). We will invest a fraction 
$\alpha$ of our money in $X$ and $1-\alpha$ in $Y$. We need to find the value of 
$\alpha$ in order to minimize the total risk, or the variance, of the investment,
*i.e.*, minimize $\textrm{var}(\alpha X + (1-\alpha) Y)$. Note that you can show
(using Calculus) that the value of $\alpha$ that minimizes the risk is given by
$$\alpha = \frac{\sigma^2_Y - \sigma_{XY}}{\sigma^2_X+\sigma^2_Y-2\sigma_{XY}},$$
where $\sigma^2_X = \textrm{X}$, $\sigma^2_Y = \textrm{Y}$, and $\sigma_{XY} = \textrm{cov}(XY)$. 

**Goal**: Collect data from $X$ and $Y$ and estimate $\alpha$. How would we 
estimate the standard deviation of $\widehat{\alpha}$? 

**Quick Aside**: This is an unrealistic exercise, but it illustrates the idea 
of a sampling distribution. We will tackle the more realistic scenario by 
implementing the bootstrap after this.

**Strategy**: Let's sample 100 pairs observations $(X, Y)$ from a normal 
population and compute $\widehat\alpha$. Replicate the experiment 1000 times in
order to compute $\widehat\alpha_1, \widehat{\alpha}_2, \dots, \widehat{\alpha}_{1000}$. This is essentially constructing a sampling 
distribution from which a estimate of the statistic's uncertainty, or standard
error (deviation), can be obtained. *Why*?

Here is the implementation in R. Note that we will use the following 
distributions $X \sim N(0, 1)$ and $Y \sim N(0, 1.25)$ with $\sigma_{XY} = .5$.
Therefore, the optimal value of $\alpha$ is 0.6.

```{r}
set.seed(98)
data <- rmvnorm(n = 1000, mean = c(0, 0), 
                sigma = matrix(c(1, .5, .5, 1.25), nc = 2))
alpha <- (var(data[, 2]) - cov(data)[1, 2]) / 
  (var(data[, 1]) + var(data[, 2]) - 2 * cov(data)[1, 2])
alpha
```

```{r}
results <- vector("numeric", 1000)
for (i in seq_along(results)) {
  df_tmp <- rmvnorm(n = 100, mean = c(0, 0), 
                    sigma = matrix(c(1, .5, .5, 1.25), nc = 2))
  alpha <- (var(df_tmp[, 2]) - cov(df_tmp)[1, 2]) / 
    (var(df_tmp[, 1]) + var(df_tmp[, 2]) - 2 * cov(df_tmp)[1, 2])
  results[i] <- alpha
}
df <- tibble(alpha = results,
             method = "true")
```

```{r}
p <- ggplot(data = df,
            aes(x = alpha))
p + geom_density(fill = "#6E0000", alpha = .3)
```

How could we use this information to compute a CI?
```{r}
summarize(df, m = mean(alpha), s = sd(alpha)) %>%
  mutate(lower =  0.6294525 - 2*s, upper = 0.6294525 + 2*s)
```

### Bootstrap Problem

This is the more realistic scenario in that we will have one set of data that we
have observed and need to estimate a standard error of some statistic. How can 
we do this in R? We essentially need to build up a sampling distribution. 

```{r}
colnames(data) <- c("x", "y")
data <- data.frame(data)
cov_tmp <- cov(data)
# alpha_1 <- (cov_tmp[2, 2] - cov_tmp[1, 2]) / 
#   (cov_tmp[1, 1] + cov_tmp[2, 2] - 2 * cov_tmp[1, 2])
alpha <- (var(data[, 2]) - cov(data)[1, 2]) / 
  (var(data[, 1]) + var(data[, 2]) - 2 * cov(data)[1, 2])
alpha
```

I am going to re-sample observations (w/ replacement) from my data set a large
number of times and recompute alpha on each "bootstrap" sample. 

```{r}
set.seed(2938)
B <- 1000 # Number of bootstrap samples
results <- vector("numeric", B)
for (i in seq_along(results)) {
  boots <- data %>%
    dplyr::sample_n(size = 100, replace = TRUE)
  alpha <- (var(boots[, 2]) - cov(boots)[1, 2]) / 
    (var(boots[, 1]) + var(boots[, 2]) - 2 * cov(boots)[1, 2])
  results[i] <- alpha
}
df_boots <- tibble(alpha = results,
                   method = "boots")
```

```{r}
p <- ggplot(data = df_boots,
            aes(x = alpha))
p + geom_density(fill = "#6e0000", alpha = .25)
```

How can we compare the two figures?

```{r}
df <- bind_rows(df, df_boots)
p <- ggplot(data = df,
            aes(x = alpha, fill = method))
p + geom_density(alpha = .25) +
  scale_fill_brewer(palette = "Dark2")
```


```{r}
df %>%
  group_by(method) %>%
  summarize(m = mean(alpha), s = sd(alpha)) %>%
  mutate(lower = 0.629 - 2*s, upper = 0.629 + 2*s)
```

### Easier Bootstrapping?

How about using the `rsample` and `purrr` package?


```{r}
set.seed(13)
resample <- rsample::bootstraps(data, times = 1000)
rsample_boots <- purrr::map_dbl(resample$splits,
                     function(x){
                       df <- analysis(x)
                       alpha <- (var(df[, 2]) - cov(df)[1, 2]) /
                         (var(df[, 1]) + var(df[, 2]) - 2 * cov(df)[1, 2])
                       return(alpha)
                       }
                     )
```

```{r}
df <- bind_rows(df,
                tibble(alpha = rsample_boots,
                       method = "rsample"))
p <- ggplot(data = df,
            aes(x = alpha, fill = method))
p + geom_density(alpha = .25) +
  scale_fill_brewer(palette = "Dark2")
```

### Another Method

```{r}
df %>%
  group_by(method) %>%
  summarize(m = mean(alpha), s = sd(alpha)) %>%
  mutate(lower = 0.629 - 2*s, upper = 0.629 + 2*s)
```

### Profiling
Which bootstrap implementation is preferred? 

```{r}
system.time({
  B <- 1000 # Number of bootstrap samples
  results <- vector("numeric", B)
  for (i in seq_along(results)) {
    boots <- data %>%
      sample_n(size = 100, replace = TRUE)
    alpha <- (var(boots[, 2]) - cov(boots)[1, 2]) / 
      (var(boots[, 1]) + var(boots[, 2]) - 2 * cov(boots)[1, 2])
    results[i] <- alpha
  }
})
```

```{r}
system.time({
  resample <- rsample::bootstraps(data, times = 1000)
  rsample_boots <- purrr::map_dbl(resample$splits,
                                  function(x){
                                    df <- as.data.frame(x)
                                    alpha <- (var(df[, 2]) - 
                                                cov(df)[1, 2]) /
                                      (var(df[, 1]) + var(df[, 2]) - 2 * cov(df)[1, 2])
                                    return(alpha)
                                  })
})
```

### Microbenchmark

```{r}
## I say don't use all the cores
all_cores <- parallel::detectCores(logical = FALSE) - 1
registerDoFuture()
cl <- parallel::makeCluster(all_cores, setup_strategy = "sequential")
plan(future::cluster, workers = cl)
```

```{r}
mbm = microbenchmark(
    purrr = {
    resample <- rsample::bootstraps(data, times = 1000)
    rsample_boots <- purrr::map_dbl(resample$splits,
                                    function(x){
                                      df <- analysis(x)
                                      alpha <- (var(df[, 2]) - 
                                                  cov(df)[1, 2]) /
                                        (var(df[, 1]) + var(df[, 2]) - 2 * cov(df)[1, 2])
                                      return(alpha)
                                    })
  },
  furrr = {
    resample <- rsample::bootstraps(data, times = 1000)
    rsample_boots <- furrr::future_map_dbl(resample$splits,
                                    function(x){
                                      df <- analysis(x)
                                      alpha <- (var(df[, 2]) - 
                                                  cov(df)[1, 2]) /
                                        (var(df[, 1]) + var(df[, 2]) - 2 * cov(df)[1, 2])
                                      return(alpha)
                                    })
  },
  old = {
    B <- 1000 # Number of bootstrap samples
    results <- vector("numeric", B)
    for (i in seq_along(results)) {
      boots <- data %>%
        sample_n(size = 100, replace = TRUE)
      alpha <- (var(boots[, 2]) - cov(boots)[1, 2]) / 
        (var(boots[, 1]) + var(boots[, 2]) - 2 * cov(boots)[1, 2])
      results[i] <- alpha
    }
  },
  times = 20,
  unit = "s"
)
```

```{r}
autoplot(mbm)
```


```{r}
p <- ggplot(data = mbm,
            aes(x = expr, y = time/1000000000, fill = expr))
p + geom_violin() +
  labs(y = "runtime (s)",
       x = "method") +
  coord_flip() +
  scale_fill_brewer(palette = "Set1") +
  theme_bw()
```


Yowza, what if we use more cores? There is a `future` implementation of `purrr`,
called `furrr`. In some cases, you can speed up the work. There is always a 
tradeoff in moving to multiple cores though.