Chapter 20 Resampling

  • NOTE: This chapter is currently be re-written and will likely change considerably in the near future. It is currently lacking in a number of ways mostly narrative.

In this chapter we introduce resampling methods, in particular cross-validation. We will highlight the need for cross-validation by comparing it to our previous approach, which was to simply set aside some “test” data that we used for evaluating a model fit using “training” data. We will now refer to these held-out samples as the validation data and this approach as the validation set approach. Along the way we’ll redefine the notion of a “test” dataset.

To illustrate the use of resampling techniques, we’ll consider a regression setup with a single predictor \(x\), and a regression function \(f(x) = x^3\). Adding an additional noise parameter, we define the entire data generating process as

\[ Y \sim N(\mu = x^3, \sigma^2 = 0.25 ^ 2) \]

We write an R function that generates datasets according to this process.

gen_sim_data = function(sample_size) {
  x = runif(n = sample_size, min = -1, max = 1)
  y = rnorm(n = sample_size, mean = x ^ 3, sd = 0.25)
  data.frame(x, y)
}

We first simulate a single dataset, which we also split into a train and validation set. Here, the validation set is 20% of the data.

set.seed(42)
sim_data = gen_sim_data(sample_size = 200)
sim_idx  = sample(1:nrow(sim_data), 160)
sim_trn  = sim_data[sim_idx, ]
sim_val  = sim_data[-sim_idx, ]

We plot this training data, as well as the true regression function.

plot(y ~ x, data = sim_trn, col = "dodgerblue", pch = 20)
grid()
curve(x ^ 3, add = TRUE, col = "black", lwd = 2)

calc_rmse = function(actual, predicted) {
  sqrt(mean((actual - predicted) ^ 2))
}

Recall that we needed this validation set because the training error was far too optimistic for highly flexible models. This would lead us to always use the most flexible model.

fit = lm(y ~ poly(x, 10), data = sim_trn)

calc_rmse(actual = sim_trn$y, predicted = predict(fit, sim_trn))
## [1] 0.2297774
calc_rmse(actual = sim_val$y, predicted = predict(fit, sim_val))
## [1] 0.2770462

20.1 Validation-Set Approach

  • TODO: consider fitting polynomial models of degree k = 1:10 to data from this data generating process
  • TODO: here, we can consider k, the polynomial degree, as a tuning parameter
  • TODO: perform simulation study to evaluate how well validation set approach works
num_sims = 100
num_degrees = 10
val_rmse = matrix(0, ncol = num_degrees, nrow = num_sims)
  • TODO: each simulation we will…
set.seed(42)
for (i in 1:num_sims) {
  # simulate data
  sim_data = gen_sim_data(sample_size = 200)
  # set aside validation set
  sim_idx = sample(1:nrow(sim_data), 160)
  sim_trn = sim_data[sim_idx, ]
  sim_val = sim_data[-sim_idx, ]
  # fit models and store RMSEs
  for (j in 1:num_degrees) {
    #fit model
    fit = glm(y ~ poly(x, degree = j), data = sim_trn)
    # calculate error
    val_rmse[i, j] = calc_rmse(actual = sim_val$y, predicted = predict(fit, sim_val))
  }
}

  • TODO: issues are hard to “see” but have to do with variability

20.2 Cross-Validation

Instead of using a single test-train split, we instead look to use \(K\)-fold cross-validation.

\[ \text{RMSE-CV}_{K} = \sum_{k = 1}^{K} \frac{n_k}{n} \text{RMSE}_k \]

\[ \text{RMSE}_k = \sqrt{\frac{1}{n_k} \sum_{i \in C_k} \left( y_i - \hat{f}^{-k}(x_i) \right)^2 } \]

  • \(n_k\) is the number of observations in fold \(k\)
  • \(C_k\) are the observations in fold \(k\)
  • \(\hat{f}^{-k}()\) is the trained model using the training data without fold \(k\)

If \(n_k\) is the same in each fold, then

\[ \text{RMSE-CV}_{K} = \frac{1}{K}\sum_{k = 1}^{K} \text{RMSE}_k \]

  • TODO: create and add graphic that shows the splitting process
  • TODO: Can be used with any metric, MSE, RMSE, class-err, class-acc

There are many ways to perform cross-validation in R, depending on the statistical learning method of interest. Some methods, for example glm() through boot::cv.glm() and knn() through knn.cv() have cross-validation capabilities built-in. We’ll use glm() for illustration. First we need to convince ourselves that glm() can be used to perform the same tasks as lm().

glm_fit = glm(y ~ poly(x, 3), data = sim_trn)
coef(glm_fit)
## (Intercept) poly(x, 3)1 poly(x, 3)2 poly(x, 3)3 
##  -0.0289307   4.5471322  -0.1527198   2.3895917
lm_fit  = lm(y ~ poly(x, 3), data = sim_trn)
coef(lm_fit)
## (Intercept) poly(x, 3)1 poly(x, 3)2 poly(x, 3)3 
##  -0.0289307   4.5471322  -0.1527198   2.3895917

By default, cv.glm() will report leave-one-out cross-validation (LOOCV).

sqrt(boot::cv.glm(sim_trn, glm_fit)$delta)
## [1] 0.2598565 0.2598342

We are actually given two values. The first is exactly the LOOCV-MSE. The second is a minor correction that we will not worry about. We take a square root to obtain LOOCV-RMSE.

In practice, we often prefer 5 or 10-fold cross-validation for a number of reason, but often most importantly, for computational efficiency.

sqrt(boot::cv.glm(sim_trn, glm_fit, K = 5)$delta)
## [1] 0.2628287 0.2617056

We repeat the above simulation study, this time performing 5-fold cross-validation. With a total sample size of \(n = 200\) each validation set has 40 observations, as did the single validation set in the previous simulations.

cv_rmse = matrix(0, ncol = num_degrees, nrow = num_sims)
set.seed(42)
for (i in 1:num_sims) {
  # simulate data, use all data for training
  sim_trn = gen_sim_data(sample_size = 200)
  # fit models and store RMSE
  for (j in 1:num_degrees) {
    #fit model
    fit = glm(y ~ poly(x, degree = j), data = sim_trn)
    # calculate error
    cv_rmse[i, j] = sqrt(boot::cv.glm(sim_trn, fit, K = 5)$delta[1])
  }
}

Polynomial Degree Mean, Val SD, Val Mean, CV SD, CV
1 0.290 0.031 0.293 0.015
2 0.291 0.031 0.295 0.014
3 0.247 0.027 0.251 0.010
4 0.248 0.028 0.252 0.010
5 0.248 0.027 0.253 0.010
6 0.249 0.027 0.254 0.011
7 0.251 0.027 0.255 0.012
8 0.252 0.027 0.257 0.011
9 0.253 0.028 0.258 0.012
10 0.255 0.027 0.259 0.012

  • TODO: differences: less variance, better selections

20.3 Test Data

The following example, inspired by The Elements of Statistical Learning, will illustrate the need for a dedicated test set which is never used in model training. We do this, if for no other reason, because it gives us a quick sanity check that we have cross-validated correctly. To be specific we will always test-train split the data, then perform cross-validation within the training data.

Essentially, this example will also show how to not cross-validate properly. It will also show can example of cross-validated in a classification setting.

calc_err = function(actual, predicted) {
  mean(actual != predicted)
}

Consider a binary response \(Y\) with equal probability to take values \(0\) and \(1\).

\[ Y \sim \text{bern}(p = 0.5) \]

Also consider \(p = 10,000\) independent predictor variables, \(X_j\), each with a standard normal distribution.

\[ X_j \sim N(\mu = 0, \sigma^2 = 1) \]

We simulate \(n = 100\) observations from this data generating process. Notice that the way we’ve defined this process, none of the \(X_j\) are related to \(Y\).

set.seed(42)
n = 200
p = 10000
x = replicate(p, rnorm(n))
y = c(rbinom(n = n, size = 1, prob = 0.5))
full_data = data.frame(y, x)

Before attempting to perform cross-validation, we test-train split the data, using half of the available data for each. (In practice, with this little data, it would be hard to justify a separate test dataset, but here we do so to illustrate another point.)

trn_idx  = sample(1:nrow(full_data), trunc(nrow(full_data) * 0.5))
trn_data = full_data[trn_idx,   ]
tst_data = full_data[-trn_idx, ]

Now we would like to train a logistic regression model to predict \(Y\) using the available predictor data. However, here we have \(p > n\), which prevents us from fitting logistic regression. To overcome this issue, we will first attempt to find a subset of relevant predictors. To do so, we’ll simply find the predictors that are most correlated with the response.

# find correlation between y and each predictor variable
correlations = apply(trn_data[, -1], 2, cor, y = trn_data$y)

While many of these correlations are small, many very close to zero, some are as large as 0.40. Since our training data has 50 observations, we’ll select the 25 predictors with the largest (absolute) correlations.

selected = order(abs(correlations), decreasing = TRUE)[1:25]
correlations[selected]
##      X4942       X867      X8617      X8044       X406      X4358      X7725 
##  0.4005771  0.3847397  0.3809371  0.3692479 -0.3571329  0.3553777 -0.3459522 
##      X1986      X3784        X77      X7010      X9354      X8450      X2355 
## -0.3448612  0.3298109 -0.3252776 -0.3242813  0.3227353  0.3220087  0.3192606 
##      X4381      X2486      X5947      X5767      X1227      X1464      X8223 
##  0.3157441  0.3149892  0.3131235  0.3114936 -0.3105052 -0.3104528  0.3084551 
##       X188      X4203      X2234      X1098 
##  0.3065491  0.3039848 -0.3036512 -0.3036153

We subset the training and test sets to contain only the response as well as these 25 predictors.

trn_screen = trn_data[c(1, selected)]
tst_screen = tst_data[c(1, selected)]

Then we finally fit an additive logistic regression using this subset of predictors. We perform 10-fold cross-validation to obtain an estimate of the classification error.

add_log_mod = glm(y ~ ., data = trn_screen, family = "binomial")
boot::cv.glm(trn_screen, add_log_mod, K = 10)$delta[1]
## [1] 0.3742339

The 10-fold cross-validation is suggesting a classification error estimate of almost 30%.

add_log_pred = (predict(add_log_mod, newdata = tst_screen, type = "response") > 0.5) * 1
calc_err(predicted = add_log_pred, actual = tst_screen$y)
## [1] 0.48

However, if we obtain an estimate of the error using the set, we see an error rate of 50%. No better than guessing! But since \(Y\) has no relationship with the predictors, this is actually what we would expect. This incorrect method we’ll call screen-then-validate.

Now, we will correctly screen-while-validating. Essentially, instead of simply cross-validating the logistic regression, we also need to cross validate the screening process. That is, we won’t simply use the same variables for each fold, we get the “best” predictors for each fold.

For methods that do not have a built-in ability to perform cross-validation, or for methods that have limited cross-validation capability, we will need to write our own code for cross-validation. (Spoiler: This is not completely true, but let’s pretend it is, so we can see how to perform cross-validation from scratch.)

This essentially amounts to randomly splitting the data, then looping over the splits. The createFolds() function from the caret() package will make this much easier.

caret::createFolds(trn_data$y, k = 10)
## $Fold01
##  [1] 17 23 27 44 45 76 85 87 93 97
## 
## $Fold02
##  [1]  6 14 15 26 37 38 55 68 69 71
## 
## $Fold03
##  [1]  3  4  7 29 39 52 54 57 59 82
## 
## $Fold04
##  [1] 19 21 40 46 48 56 73 78 91 96
## 
## $Fold05
##  [1] 25 34 36 58 61 65 66 75 83 89
## 
## $Fold06
##  [1]  2  9 10 62 74 79 80 90 92 98
## 
## $Fold07
##  [1]  8 31 32 41 43 53 60 67 88 95
## 
## $Fold08
##  [1] 12 18 33 35 42 49 51 64 84 94
## 
## $Fold09
##  [1]  11  13  16  20  28  47  50  77  99 100
## 
## $Fold10
##  [1]  1  5 22 24 30 63 70 72 81 86
# use the caret package to obtain 10 "folds"
folds = caret::createFolds(trn_data$y, k = 10)

# for each fold
# - pre-screen variables on the 9 training folds
# - fit model to these variables
# - get error on validation fold
fold_err = rep(0, length(folds))

for (i in seq_along(folds)) {

  # split for fold i  
  trn_fold = trn_data[-folds[[i]], ]
  val_fold = trn_data[folds[[i]], ]

  # screening for fold i  
  correlations = apply(trn_fold[, -1], 2, cor, y = trn_fold[,1])
  selected = order(abs(correlations), decreasing = TRUE)[1:25]
  trn_fold_screen = trn_fold[ , c(1, selected)]
  val_fold_screen = val_fold[ , c(1, selected)]

  # error for fold i  
  add_log_mod = glm(y ~ ., data = trn_fold_screen, family = "binomial")
  add_log_prob = predict(add_log_mod, newdata = val_fold_screen, type = "response")
  add_log_pred = ifelse(add_log_prob > 0.5, yes = 1, no = 0)
  fold_err[i] = mean(add_log_pred != val_fold_screen$y)
  
}

# report all 10 validation fold errors
fold_err
##  [1] 0.4 0.9 0.6 0.4 0.6 0.3 0.7 0.5 0.6 0.6
# properly cross-validated error
# this roughly matches what we expect in the test set
mean(fold_err)
## [1] 0.56
  • TODO: note that, even cross-validated correctly, this isn’t a brilliant variable selection procedure. (it completely ignores interactions and correlations among the predictors. however, if it works, it works.) next chapters…

20.4 Bootstrap

ISL discusses the bootstrap, which is another resampling method. However, it is less relevant to the statistical learning tasks we will encounter. It could be used to replace cross-validation, but encounters significantly more computation.

It could be more useful if we were to attempt to calculate the bias and variance of a prediction (estimate) without access to the data generating process. Return to the bias-variance tradeoff chapter and think about how the bootstrap could be used to obtain estimates of bias and variance with a single dataset, instead of repeated simulated datasets.

20.5 Which \(K\)?

  • TODO: LOO vs 5 vs 10
  • TODO: bias and variance

20.6 Summary

  • TODO: using cross validation for: tuning, error estimation

20.8 rmarkdown

The rmarkdown file for this chapter can be found here. The file was created using R version 4.0.2.