Cross-Validation, Tuning, and Thresholds

1. Why one validation split can be fragile

Performance measured on a small set depends on which observations happened to be sampled. K-fold cross-validation divides the training portion into k folds, repeatedly fits on k−1 folds, and predicts the remaining fold. This uses training data more efficiently for model comparison.

Example: Split ten typed ids into 5 folds and list which ids each fold would hold out.

library(tidymodels)
set.seed(1)
tiny_folds <- vfold_cv(tibble(id = 1:10), v = 5)
purrr::map(tiny_folds$splits, ~ assessment(.x)$id)

Output

[[1]]
[1] 4 9

[[2]]
[1] 3 5

[[3]]
[1]  7 10

[[4]]
[1] 1 2

[[5]]
[1] 6 8

Each fold holds out two ids, and every id appears in exactly one fold. In fold 1 the model would be fitted on the other eight ids and scored on the two held-out ones; the five scores are then averaged.

This chapter uses only the training rows. The test rows were already scored once in Chapter 15, so they are no longer untouched.

2. Set up a cross-validation workflow

library(tidymodels)
ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE) |>
  mutate(diabetes = factor(diabetes, levels = c("Yes", "No")),
         clinic_site = factor(clinic_site), sex = factor(sex),
         family_history = factor(family_history))
set.seed(500)
split <- initial_split(ehr, prop = 0.8, strata = diabetes)
train <- training(split)
rec <- recipe(diabetes ~ age + bmi + sbp + pulse + family_history + sex + clinic_site,
              data = train) |>
  step_impute_median(all_numeric_predictors()) |>
  step_impute_mode(all_nominal_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_zv(all_predictors())
wf <- workflow() |>
  add_recipe(rec) |>
  add_model(logistic_reg() |> set_engine("glm") |> set_mode("classification"))
set.seed(500)
folds <- vfold_cv(train, v = 5, strata = diabetes)

Example: Print the folds and count the held-out rows in each.

folds

Output

#  5-fold cross-validation using stratification 
# A tibble: 5 x 2
  splits           id   
  <list>           <chr>
1 <split [384/96]> Fold1
2 <split [384/96]> Fold2
3 <split [384/96]> Fold3
4 <split [384/96]> Fold4
5 <split [384/96]> Fold5

purrr::map_int(folds$splits, ~ nrow(assessment(.x)))

Output

[1] 96 96 96 96 96

Each of the five splits shows 384/96: 384 rows to fit on and 96 rows to assess. All five assessment folds hold 96 rows, and 5 x 96 is the 480 training rows.

Every analysis fold refits the recipe on that fold’s analysis rows. That prevents the fold’s assessment rows from informing its imputation medians.

3. Collect out-of-fold predictions

set.seed(500)
cv_fit <- fit_resamples(
  wf,
  resamples = folds,
  metrics = metric_set(roc_auc, pr_auc, accuracy, sens, spec),
  control = control_resamples(save_pred = TRUE)
)
collect_metrics(cv_fit)

Output

# A tibble: 5 x 6
  .metric  .estimator   mean     n std_err .config        
  <chr>    <chr>       <dbl> <int>   <dbl> <chr>          
1 accuracy binary     0.904      5 0.0166  pre0_mod0_post0
2 pr_auc   binary     0.322      5 0.0500  pre0_mod0_post0
3 roc_auc  binary     0.829      5 0.0447  pre0_mod0_post0
4 sens     binary     0.0987     5 0.0487  pre0_mod0_post0
5 spec     binary     0.989      5 0.00596 pre0_mod0_post0

cv_predictions <- collect_predictions(cv_fit)
head(cv_predictions)

Output

# A tibble: 6 x 7
  .pred_class .pred_Yes .pred_No id    diabetes  .row .config        
  <fct>           <dbl>    <dbl> <chr> <fct>    <int> <chr>          
1 No            0.0117     0.988 Fold1 No           5 pre0_mod0_post0
2 No            0.00544    0.995 Fold1 No          10 pre0_mod0_post0
3 No            0.0270     0.973 Fold1 No          18 pre0_mod0_post0
4 No            0.00505    0.995 Fold1 No          23 pre0_mod0_post0
5 No            0.0938     0.906 Fold1 No          26 pre0_mod0_post0
6 No            0.0721     0.928 Fold1 No          32 pre0_mod0_post0

Example: Read the table that collect_metrics() printed. The mean ROC AUC over the five folds is 0.829 with a standard error of 0.045, so the true value could plausibly be a good bit higher or lower. Accuracy is 0.904 but sensitivity is only 0.099 and specificity is 0.989 (these are the means of the five fold-level values): at the default 0.5 cutoff the model says No almost every time, which is why accuracy looks fine. head(cv_predictions) shows one row per training patient, scored by a model that did not train on that patient.

The fold estimates describe model performance during development. They do not replace a truly independent test. If you try many models and keep improving their CV score, you can adaptively overfit the CV process; record the search process and avoid treating the best observed CV mean as unbiased.

4. Choose a threshold on training folds only

Let’s illustrate a synthetic target of at least 80% sensitivity on pooled out-of-fold predictions. This target is an instructional choice, not a medical recommendation.

candidate_thresholds <- seq(0.05, 0.95, by = 0.05)
positive <- cv_predictions$diabetes == "Yes"
scores <- purrr::map_dfr(candidate_thresholds, function(t) {
  predicted_positive <- cv_predictions$.pred_Yes >= t
  tibble(
    threshold = t,
    sensitivity = sum(positive & predicted_positive) / sum(positive),
    specificity = sum(!positive & !predicted_positive) / sum(!positive)
  )
})
scores

Output

# A tibble: 19 x 3
   threshold sensitivity specificity
       <dbl>       <dbl>       <dbl>
 1      0.05      0.756        0.572
 2      0.1       0.644        0.754
 3      0.15      0.578        0.832
 4      0.2       0.511        0.887
 5      0.25      0.356        0.920
 6      0.3       0.289        0.949
 7      0.35      0.244        0.968
 8      0.4       0.222        0.984
 9      0.45      0.178        0.984
10      0.5       0.0889       0.989
11      0.55      0.0444       0.993
12      0.6       0.0222       1    
13      0.65      0.0222       1    
14      0.7       0.0222       1    
15      0.75      0.0222       1    
16      0.8       0            1    
17      0.85      0            1    
18      0.9       0            1    
19      0.95      0            1

eligible <- scores |>
  filter(sensitivity >= 0.80) |>
  arrange(desc(specificity), desc(threshold))
locked_threshold <- if (nrow(eligible) > 0) eligible$threshold[1] else 0.5
locked_threshold

Output

[1] 0.5

If no cutoff meets the goal, report that failure rather than invent a compliant result. Our simple fallback of 0.5 is an implementation convenience for this exercise, not a justified clinical threshold. This code does not use the test rows.

Example: Read the scores table above. Sensitivity falls as the cutoff rises: 0.756 at 0.05, 0.089 at 0.50, and 0 at 0.80. These are pooled over all 480 out-of-fold rows, so 0.089 differs slightly from the fold-average 0.099 printed by collect_metrics(). No row reaches the 0.80 target, so eligible is empty and locked_threshold is the fallback 0.5. The honest report is that the target is not met on this 0.05 grid; Try it 3 searches a finer 0.01 grid.

5. What is hyperparameter tuning?

A parameter is learned from data (for example, a regression coefficient). A hyperparameter is a setting chosen by the analyst (for example, how deep a decision tree may grow or how many observations a node must hold before it can be split further). We compare settings via folds and choose one according to a metric chosen in advance.

Example: Fit the logistic workflow on train and look at a few coefficients (parameters the data decided), then write down a forest setting we choose ourselves.

fit(wf, data = train) |> tidy() |> select(term, estimate) |> slice_head(n = 3)

Output

# A tibble: 3 x 2
  term        estimate
  <chr>          <dbl>
1 (Intercept)  15.7   
2 age          -0.0321
3 bmi          -0.0608

rand_forest(min_n = 5) |> set_mode("classification")

Output

Random Forest Model Specification (classification)

Main Arguments:
  min_n = 5

Computational engine: ranger

The three coefficients came out of the fitting. min_n = 5 was typed by us before any fitting; the next chapter tries other values with the folds.

When selecting a final model, decide on CV performance, complexity, interpretability, and safety before consulting test outcomes. After a single held-out test, poor performance is a finding to report; further iteration needs a newly independent evaluation set or clear disclosure that the old set has been used for development.

6. Calibration and confidence intervals

Out-of-fold predictions can also inform calibration checks and interval estimation. For small positive counts, simple point estimates can be highly uncertain. Do not choose a model solely on a tiny apparent difference in CV AUC. Later studies should consider bootstrap confidence intervals, subgroup metrics, temporal and external validation.

Example: Put a 95% confidence interval around the out-of-fold sensitivity at a 0.5 cutoff.

n_pos <- sum(positive)
hits <- sum(positive & cv_predictions$.pred_Yes >= 0.5)
c(hits = hits, positives = n_pos)

Output

     hits positives 
        4        45

round(binom.test(hits, n_pos)$conf.int, 3)

Output

[1] 0.025 0.212
attr(,"conf.level")
[1] 0.95

Only 4 of the 45 positive rows were caught (pooled sensitivity 0.089, not the fold average of 0.099), and the interval runs from about 0.025 to 0.212. With so few positives, any single sensitivity number is very uncertain.

Try it

Change the number of folds to 4. How does each assessment fold’s size change?

Show solution

Answer: Each assessment fold grows from 96 rows (5 folds) to 120 rows (4 folds), because 480 training rows are split into fewer, larger pieces. Each fitting set shrinks from 384 to 360 rows.

Step by step

Build 4 folds with the same seed and stratification.
Count the assessment rows in each fold.

set.seed(500)
folds4 <- vfold_cv(train, v = 4, strata = diabetes)
folds4

Output

#  4-fold cross-validation using stratification 
# A tibble: 4 x 2
  splits            id   
  <list>            <chr>
1 <split [360/120]> Fold1
2 <split [360/120]> Fold2
3 <split [360/120]> Fold3
4 <split [360/120]> Fold4

purrr::map_int(folds4$splits, ~ nrow(assessment(.x)))

Output

[1] 120 120 120 120

3. The splits read 360/120: 360 rows to fit on and 120 to assess, four times. That is 480 / 4 = 120 rows per fold. Fewer folds means fewer fits (faster), but each model sees less of the training data.

Calculate the number of out-of-fold observations.

Show solution

Answer: 480. Each training row is held out exactly once, so there is one out-of-fold prediction per training row.

Step by step

Count the rows in the collected predictions.
Compare with the training rows, and check that no row appears twice.

nrow(cv_predictions)

Output

[1] 480

nrow(train)

Output

[1] 480

n_distinct(cv_predictions$.row)

Output

[1] 480

3. All three numbers are 480, so each training patient has exactly one prediction made by a model that never saw that patient. (With repeated CV you would get one per repeat.)

Change the illustrative sensitivity target to 0.90; describe what typically happens to the selected cutoff and specificity.

Show solution

Answer: A higher sensitivity target forces a lower cutoff, and a lower cutoff gives lower specificity. In this chapter no cutoff on the 0.05 grid reaches 0.90 (or even 0.80), so the rule falls back to 0.5; on a finer grid the 0.90 target is met only at a cutoff of 0.02 or lower, and at 0.02 specificity is only 0.356.

Step by step

Apply the chapter’s rule with a 0.90 target to the scores table.

scores |>
  filter(sensitivity >= 0.90) |>
  nrow()

Output

[1] 0

2. Zero rows qualify, so locked_threshold would be the fallback 0.5. Now build a finer grid and compare the two targets.

fine <- purrr::map_dfr(seq(0.01, 0.50, by = 0.01), function(t) {
  predicted_positive <- cv_predictions$.pred_Yes >= t
  tibble(
    threshold = t,
    sensitivity = sum(positive & predicted_positive) / sum(positive),
    specificity = sum(!positive & !predicted_positive) / sum(!positive)
  )
})
pick <- function(target) {
  fine |>
    filter(sensitivity >= target) |>
    arrange(desc(specificity), desc(threshold)) |>
    slice(1)
}
bind_rows(`0.80` = pick(0.80), `0.90` = pick(0.90), .id = "target")

Output

# A tibble: 2 x 4
  target threshold sensitivity specificity
  <chr>      <dbl>       <dbl>       <dbl>
1 0.80        0.04       0.8         0.520
2 0.90        0.02       0.911       0.356

3. For a target of 0.80 the best cutoff is 0.04 (sensitivity 0.800, specificity 0.520). For 0.90 it drops to 0.02 (sensitivity 0.911, specificity 0.356). Catching more of the positives costs many more false alarms, and these cutoffs are tiny because the model’s scores for Yes patients are low.

Explain why imputation must be learned separately within each fold.

Show solution

Answer: The median used to fill missing values must come only from the rows the model is fitted on. If it were computed on all rows, the held-out fold would influence the fill value, so the held-out fold would no longer be unseen (data leakage) and the CV score would be slightly too optimistic.

Step by step

Take fold 1 and compute the median bmi from three row sets: all training rows, the fold’s fitting (analysis) rows, and the held-out (assessment) rows.

fold1 <- folds$splits[[1]]
c(
  all_train = median(train$bmi, na.rm = TRUE),
  analysis = median(analysis(fold1)$bmi, na.rm = TRUE),
  assessment = median(assessment(fold1)$bmi, na.rm = TRUE)
)

Output

 all_train   analysis assessment 
      30.5       30.2       31.8

2. The three medians differ. The all-training median already contains information from the held-out rows; the analysis median does not. 3. fit_resamples() re-learns the recipe on each fold’s analysis rows only, so each fold uses its own analysis median. The same holds for dummy columns and any scaling.

Explain the difference between CV tuning and final locked-test assessment.

Show solution

Answer: CV tuning is for choosing: it compares models and settings using only training rows, and you may repeat it. The locked test assessment is for reporting: it scores one already-chosen model and cutoff once on rows that played no part in the choices.

Step by step

In CV, every row is used for both fitting and scoring (in different folds), and the score is used to pick between options. Look at the CV ROC AUC and how uncertain it is.

collect_metrics(cv_fit) |> filter(.metric == "roc_auc")

Output

# A tibble: 1 x 6
  .metric .estimator  mean     n std_err .config        
  <chr>   <chr>      <dbl> <int>   <dbl> <chr>          
1 roc_auc binary     0.829     5  0.0447 pre0_mod0_post0

2. The mean is 0.829 with a standard error of 0.045. This is a development estimate: the more options we compare, the more it flatters the winner. 3. The test set is not used in any of those choices, so scoring it once gives an estimate that has not been tuned to. Once you have tried several things on it, it becomes a second tuning set, and you need fresh data for an honest final score.

Checkpoint

You can run CV on training data, interpret out-of-fold predictions, choose an instructional threshold without touching test data, and identify adaptive overfitting risk.