11  Resampling in rtemis

Each resampling design of Chapter 10 has a setup_*() function in rtemis. resample() applies a configuration to a dataset and returns the training indices of every resample; train() applies it through its outer_resampling_config argument. This chapter creates resamples of the heart failure data, measures how the variability of a performance estimate depends on the design, and shows the effect of grouped resampling on the Parkinson’s data.

11.1 Resampling configurations

Table 11.1: Resampling configurations
Design Function Main arguments
k-fold cross-validation setup_KFold() n_resamples (the number of folds), stratify_var, id_strat, seed
Repeated stratified subsampling setup_StratSub() n_resamples, train_p, stratify_var, id_strat, seed
Stratified bootstrap setup_StratBoot() n_resamples, stratify_var, id_strat, seed
Bootstrap setup_Bootstrap() n_resamples, id_strat, seed
Leave-one-out setup_LOOCV() —

k-fold cross-validation, subsampling and the stratified bootstrap stratify on the outcome unless stratify_var names another column; a numeric outcome is stratified within strat_n_bins quantile bins. id_strat names a column of group identifiers, such as patient IDs, and keeps all cases of a group in the same side of every split.

11.2 Creating resamples

hf <- read_heart_failure()
folds <- resample(hf, setup_KFold(n_resamples = 10L, seed = 2026L))
folds
<Resampler>
     type: KFold
resamples: 
            Fold_1: 1, 2, 3, 4, ...
            Fold_2: 1, 2, 3, 5, ...
            Fold_3: 1, 2, 3, 4, ...
            Fold_4: 2, 3, 4, 5, ...
            Fold_5: 1, 3, 4, 5, ...
            Fold_6: 1, 2, 3, 4, ...
            Fold_7: 1, 2, 4, 5, ...
            Fold_8: 1, 2, 3, 4, ...
            Fold_9: 1, 2, 3, 4, ...
           Fold_10: 1, 2, 3, 4, ...
   config:
           <KFoldConfig>
            n_resamples: 10
           stratify_var: NULL
           strat_n_bins: 4
               id_strat: NULL
                   seed: 2026

The Resampler object holds the training indices of each resample in its resamples property. The test cases of a resample are the cases not in its training set:

fold_summary <- rbindlist(lapply(names(folds@resamples), function(r) {
  test_rows <- setdiff(seq_len(nrow(hf)), folds@resamples[[r]])
  data.table(
    resample = r,
    n_training = length(folds@resamples[[r]]),
    n_test = length(test_rows),
    deaths_test = sum(hf[test_rows, death] == "Yes")
  )
}))
fold_summary
    resample n_training n_test deaths_test
      <char>      <int>  <int>       <int>
 1:   Fold_1        260     30           6
 2:   Fold_2        261     29           5
 3:   Fold_3        262     28           5
 4:   Fold_4        260     30           6
 5:   Fold_5        262     28           5
 6:   Fold_6        261     29           5
 7:   Fold_7        261     29           6
 8:   Fold_8        261     29           5
 9:   Fold_9        262     28           5
10:  Fold_10        260     30           6

Stratification places between 5 and 6 of the 54 deaths in each test fold. Every case is a test case exactly once:

test_count <- tabulate(
  unlist(lapply(folds@resamples, function(idx) setdiff(seq_len(nrow(hf)), idx))),
  nbins = nrow(hf)
)
table(test_count)
test_count
  1 
290 

A bootstrap resample has as many training indices as cases, with repetitions:

boot <- resample(hf, setup_StratBoot(n_resamples = 10L, seed = 2026L))
idx <- boot@resamples[[1]]
c(
  n_training = length(idx),
  n_distinct = length(unique(idx)),
  n_out_of_bag = nrow(hf) - length(unique(idx))
)
  n_training   n_distinct n_out_of_bag 
         290          217           73 

11.3 The variability of the estimate

A resampling configuration passed to train() as outer_resampling_config determines how the model is assessed. Two designs are compared here for the logistic regression of the heart failure data. The first is a single 75/25 holdout split; with setup_StratSub(n_resamples = 50L), one call fits 50 models on 50 independent random splits, and the AUC of each split is the estimate that a single holdout split would have produced:

mod_holdout <- train(
  hf,
  hyperparameters = setup_GLM(),
  outer_resampling_config = setup_StratSub(
    n_resamples = 50L, train_p = 0.75, seed = 2026L
  ),
  execution_config = setup_SerialExecution(seed = 2026L),
  verbosity = 0L
)
auc_holdout <- vapply(
  mod_holdout@metrics_test@res_metrics,
  function(m) m@metrics[["overall"]][["auc"]],
  numeric(1)
)

The second is 10-fold cross-validation, repeated with 50 different random fold assignments; each repetition yields one cross-validation estimate, the mean of its 10 test-fold AUCs:

auc_cv <- vapply(seq_len(50L), function(s) {
  mod <- train(
    hf,
    hyperparameters = setup_GLM(),
    outer_resampling_config = setup_KFold(n_resamples = 10L, seed = s),
    execution_config = setup_SerialExecution(seed = 2026L),
    verbosity = 0L
  )
  mod@metrics_test@mean_metrics[["auc"]]
}, numeric(1))
draw_boxplot(
  list(
    `Single 75/25 split` = auc_holdout,
    `10-fold cross-validation` = auc_cv
  ),
  ylab = "Estimated AUC",
  boxpoints = "all"
)

The 50 holdout estimates range from 0.567 to 0.812, with a standard deviation of 0.058. The 50 cross-validation estimates range from 0.668 to 0.739, with a standard deviation of 0.015. The AUC reported by a single holdout split of these data depends strongly on which patients fall in the test set; a 10-fold cross-validation estimate depends much less on the random division. Both distributions describe variation due to the division of one dataset, and neither includes the variation that a different sample of 290 patients would produce.

11.4 Resampling in train()

The training and test metrics of every resample are stored in the fitted object (Section 9.1). The outer resampler, with its indices, is stored in outer_resampler, so the assessment can be reproduced and inspected:

mod_cv <- train(
  hf,
  hyperparameters = setup_GLM(),
  outer_resampling_config = setup_KFold(n_resamples = 10L, seed = 2026L),
  execution_config = setup_SerialExecution(seed = 2026L),
  question = hf_question
)
<Resampled Classification Model>
GLM (GLM (generalized linear model))
⟳ Tested using 10 independent folds.

  <Resampled Classification Training Metrics>
  Aggregate Confusion Matrix across resamples.
                     Predicted
          Reference  Yes  No    
                Yes  111   375
                 No   62  2062

  Showing mean (sd) across resamples.
        Sensitivity: 0.228 (0.025)
        Specificity: 0.971 (0.006)
  Balanced Accuracy: 0.600 (0.011)
                PPV: 0.644 (0.042)
                NPV: 0.846 (3.8e-03)
                 F1: 0.336 (0.027)
           Accuracy: 0.833 (4.7e-03)
                AUC: 0.780 (0.010)
        Brier Score: 0.123 (1.6e-03)

  <Resampled Classification Test Metrics>
  Aggregate Confusion Matrix across resamples.
                     Predicted
          Reference  Yes  No   
                Yes    9   45
                 No    9  227

  Showing mean (sd) across resamples.
        Sensitivity: 0.170 (0.147)
        Specificity: 0.962 (0.047)
  Balanced Accuracy: 0.566 (0.065)
                PPV: NA (NA)
                NPV: 0.836 (0.029)
                 F1: NA (NA)
           Accuracy: 0.814 (0.040)
                AUC: 0.709 (0.117)
        Brier Score: 0.140 (0.015)
identical(mod_cv@outer_resampler@resamples, folds@resamples)
[1] TRUE

The seed of the resampling configuration determines the division into folds; the seed of the execution configuration determines any randomness in model fitting. Both are recorded with the model.

11.5 Leave-one-out cross-validation

Leave-one-out cross-validation fits one model per case. For the diabetes data, it fits 442 linear models, each tested on one patient:

db <- read_diabetes()
mod_loocv <- train(
  db,
  hyperparameters = setup_GLM(),
  outer_resampling_config = setup_LOOCV(),
  execution_config = setup_SerialExecution(seed = 2026L),
  verbosity = 0L
)
oof_loocv <- oof_predictions(mod_loocv)
regression_metrics(oof_loocv[["observed"]], oof_loocv[["predicted"]])
<Regression Metrics>
 MAE: 44.36
 MSE: 3001.75
RMSE: 54.79
  R²: 0.49

Per-resample \(R^2\) is undefined for a single test case, so the metrics are computed on the pooled out-of-fold predictions.

11.6 Grouped resampling

The Parkinson’s data contain between 101 and 168 recordings per patient (Section 7.4). Two linear models of total UPDRS are assessed by 5-fold cross-validation. The first assigns recordings to folds at random, ignoring the patient; the column subject is removed so that it is not used as a feature:

pk <- read_parkinsons()
mod_pk_random <- train(
  pk[, !"subject"],
  hyperparameters = setup_GLM(),
  outer_resampling_config = setup_KFold(n_resamples = 5L, seed = 2026L),
  execution_config = setup_SerialExecution(seed = 2026L),
  question = parkinsons_question,
  verbosity = 0L
)

The second assigns patients to folds with id_strat = "subject". train() uses the named column to group the resamples and excludes it from the features:

mod_pk_grouped <- train(
  pk,
  hyperparameters = setup_GLM(),
  outer_resampling_config = setup_KFold(
    n_resamples = 5L, id_strat = "subject", seed = 2026L
  ),
  execution_config = setup_SerialExecution(seed = 2026L),
  question = parkinsons_question,
  verbosity = 0L
)

Every patient’s recordings are in a single test fold:

test_subjects <- lapply(mod_pk_grouped@outer_resampler@resamples, function(idx) {
  unique(pk[-idx, subject])
})
any(duplicated(unlist(test_subjects)))
[1] FALSE
lengths(test_subjects)
Fold_1 Fold_2 Fold_3 Fold_4 Fold_5 
    10      7      8      7     10 

The mean test \(R^2\) is 0.163 with random assignment of recordings and -0.127 with assignment by patient:

c(
  random_recordings = rsq_random,
  grouped_by_patient = rsq_grouped
)
 random_recordings grouped_by_patient 
            0.1631            -0.1265 

With random assignment, every patient in a test fold also contributes about four fifths of their recordings to the training set, and the model is assessed on its ability to predict other recordings of patients it has seen. With assignment by patient, it is assessed on patients it has not seen, which is the situation of the question. The negative \(R^2\) shows that the linear model, applied to a new patient, predicts total UPDRS less accurately than the mean of that test fold. The difference between the two estimates is a measure of how much the random design overstates performance on new patients.