43 Calibration in rtemis

calibrate() fits a calibration map to a classification model’s predicted probabilities, using setup_Isotonic() by default; available_calibration() lists the algorithms that produce non-decreasing maps. draw_calibration() plots calibration curves. This chapter assesses the calibration of models from earlier chapters, and recalibrates the weighted elastic net of Chapter 41 with calibration cases from cross-validation.

hf <- read_heart_failure()
cv <- setup_KFold(n_resamples = 10L, seed = 2026L)
ex <- setup_SerialExecution(seed = 2026L)
available_calibration()

      Isotonic: Isotonic regression
  MonotonicHAL: Monotonic Highly Adaptive Lasso (hal9001)

43.1 Calibration of several models

The calibration intercept and slope are estimated by a logistic regression of the outcome on the logit of the out-of-fold predicted probabilities (Section 42.1):

calibration_summary <- function(observed, prob) {
  y <- as.integer(observed == "Yes")
  eps <- 1e-6
  lp <- qlogis(pmin(pmax(prob, eps), 1 - eps))
  fit <- glm(y ~ lp, family = binomial)
  c(
    mean_prob = mean(prob),
    observed_rate = mean(y),
    slope = coef(fit)[["lp"]],
    brier = mean((y - prob)^2)
  )
}
models_hf <- list(
  GLM = setup_GLM(),
  CART = setup_CART(maxdepth = 4L),
  Ranger = setup_Ranger(min_node_size = 20L, seed = 2026L),
  `GLMNET, ifw` = setup_GLMNET(alpha = 0, which_lambda_cv = "lambda.min", ifw = TRUE)
)
oof_hf <- lapply(models_hf, function(hp) {
  mod <- train(
    hf,
    hyperparameters = hp,
    outer_resampling_config = cv,
    execution_config = ex,
    verbosity = 0L
  )
  oof_predictions(mod)
})
cal_table <- rbindlist(lapply(oof_hf, function(o) {
  as.list(calibration_summary(o[["observed"]], o[["prob"]]))
}), idcol = "model")
cal_table

         model mean_prob observed_rate    slope  brier
        <char>     <num>         <num>    <num>  <num>
1:         GLM    0.1850        0.1862  0.65056 0.1396
2:        CART    0.1703        0.1862 -0.02688 0.1981
3:      Ranger    0.1892        0.1862  0.90697 0.1362
4: GLMNET, ifw    0.4540        0.1862  1.39350 0.2092

draw_calibration(
  lapply(oof_hf, `[[`, "observed"),
  lapply(oof_hf, `[[`, "prob")
)

The logistic regression’s mean predicted probability matches the observed rate, and its slope, 0.65, is below 1: with 11 coefficients and about 49 deaths per training fold, its predictions are somewhat too extreme. The depth-4 tree has the lowest slope, -0.03, since its leaves contain few patients and some predict probabilities of 0 or 1. The ridge model trained with inverse frequency weights has a mean predicted probability of 0.454 against an observed rate of 0.186.

43.2 Recalibration

The procedure of Section 42.4 requires three sets of cases: those that fit the model, those that fit the calibration map, and those that assess the result. The data are divided once into a development set of 75% of the patients and a test set of 25%:

split <- resample(hf, setup_StratSub(n_resamples = 1L, train_p = 0.75, seed = 2026L))
hf_dev <- hf[split[[1]]]
hf_test <- hf[-split[[1]]]
hp <- setup_GLMNET(alpha = 0, which_lambda_cv = "lambda.min", ifw = TRUE)

A 10-fold cross-validation on the development set provides an out-of-fold prediction for every development case, made by a model that did not see it. These are the calibration cases:

mod_dev_cv <- train(
  hf_dev,
  hyperparameters = hp,
  outer_resampling_config = setup_KFold(n_resamples = 10L, seed = 2027L),
  execution_config = ex,
  verbosity = 0L
)
oof_dev <- oof_predictions(mod_dev_cv)

The final model is fitted to the whole development set and evaluated on the test set. calibrate() fits the isotonic map to the out-of-fold predictions and returns a calibrated model, whose predict() applies the model and then the map:

mod_final <- train(
  hf_dev,
  dat_test = hf_test,
  hyperparameters = hp,
  execution_config = ex,
  question = hf_question,
  verbosity = 0L
)
mod_calibrated <- calibrate(
  mod_final,
  predicted_probabilities = oof_dev[["prob"]],
  true_labels = oof_dev[["observed"]],
  verbosity = 0L
)
prob_raw <- as.numeric(predict(mod_final, features(hf_test)))
prob_cal <- as.numeric(predict(mod_calibrated, features(hf_test)))
rbind(
  uncalibrated = calibration_summary(hf_test[["death"]], prob_raw),
  calibrated = calibration_summary(hf_test[["death"]], prob_cal)
)

             mean_prob observed_rate  slope  brier
uncalibrated    0.4563        0.1918 1.1992 0.2173
calibrated      0.1813        0.1918 0.7548 0.1406

On the 73 test patients, the calibrated probabilities have a mean close to the observed death rate and a lower Brier score. The AUC changes from 0.642 to 0.657: the isotonic map preserves the ranking but merges nearby probabilities into ties. With 14 deaths in the test set, the slope is estimated imprecisely, and a calibration curve on these patients would rest on very few events per bin.

setup_MonotonicHAL() fits a smooth, strictly increasing map instead, which leaves the AUC unchanged:

mod_calibrated_hal <- calibrate(
  mod_final,
  hyperparameters = setup_MonotonicHAL(seed = 2026L),
  predicted_probabilities = oof_dev[["prob"]],
  true_labels = oof_dev[["observed"]],
  verbosity = 0L
)
prob_hal <- as.numeric(predict(mod_calibrated_hal, features(hf_test)))
calibration_summary(hf_test[["death"]], prob_hal)

    mean_prob observed_rate         slope         brier 
       0.1825        0.1918        1.2645        0.1441

The calibration map is fitted to 217 out-of-fold predictions with 40 deaths. A larger calibration set estimates the map more precisely, and each recalibration spends calibration cases on a correction; the unweighted logistic regression and random forest of Section 43.1 were calibrated in the mean without one.