# Proper way to get an averaged calibration curve from nested resampling

**URL:** <https://discourse.datamethods.org/t/proper-way-to-get-an-averaged-calibration-curve-from-nested-resampling/6207>\
**Category:** model validation\
**Tags:** calibration, machine-learning, prediction\
**Created:** [November 12, 2022, 2:53pm UTC](https://discourse.datamethods.org/t/proper-way-to-get-an-averaged-calibration-curve-from-nested-resampling/6207 "2022-11-12T14:53:23Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![carldelfin](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/carldelfin/32/3380_2.png) [@carldelfin](https://discourse.datamethods.org/u/carldelfin)\
**Post date:** [November 12, 2022, 2:53pm UTC](https://discourse.datamethods.org/t/proper-way-to-get-an-averaged-calibration-curve-from-nested-resampling/6207/1 "2022-11-12T14:53:23Z")

</div>

Dear all,

I am working on some projects where the overarching goal is to assess whether a selection of machine learning algorithms perform better than “traditional” statistical methods such as logistic regression (with and without regularization) in terms of (binary) predictive performance, similar to [this 2021 study](https://journals.sagepub.com/doi/full/10.1177/09622802211002867) by Austin, Harrell, and Steyerberg.

In my case, most relevant datasets range in sample size from ~500 to ~10 000 and the number of predictors from ~10 to ~50. I work with the [Tidymodels](https://www.tidymodels.org/) meta-package and tune the models via [Bayesian optimization](https://tune.tidymodels.org/reference/tune_bayes.html). I use a nested resampling approach — _k_-fold cross-validation repeated _n_ times for both the inner and outer loops (although possibly with different _k_ and _n_ for inner and outer) — to obtain estimates of predictive performance. For single value estimates, I can just populate a matrix and calculate an average, like so:

| log\_loss | outer\_repeat | outer\_fold |
| --- | --- | --- |
| 0.43 | 1 | 1 |
| 0.40 | 1 | 2 |
| 0.39 | 1 | 3 |
| … | _n_ | _k_ |

And then:

```auto
mean_log <- mean(c(0.43, 0.40, 0.39 ...))

```

But — and herein lies my main concern — I also want to plot some kind of averaged calibration curve. My current solution is that for each prediction on held-out data in the outer loop, I first get model predictions:

```auto
test_pred <- predict(final_mod_fit, data_test, type = "prob") %>%
        bind_cols(data_test) %>%
        rename(prob_yes = .pred_yes) %>%
        select(prob_yes, outcome) %>%
        mutate(outcome_num = ifelse(outcome == "yes", 1, 0))

```

`test_pred` looks like this:

| prob\_yes | outcome | outcome\_num |
| --- | --- | --- |
| 0.03 | no | 0 |
| 0.45 | yes | 1 |
| 0.67 | yes | 1 |
| … | … | … |

Then I use [stats::lowess](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/lowess) to get the x and y coordinates of the smoothed calibration curve, with the predicted probabilities and observed outcomes in `test_pred` as input:

```auto
cal <- data.frame(with(test_pred, lowess(prob_yes, outcome_num, iter = 0)))

```

I thus get one `cal` dataframe for each _k_ in each _n_. After some processing to merge them all, I get a dataframe structured like this:

| x | y | outer\_repeat | outer\_fold |
| --- | --- | --- | --- |
| 0.01 | 0.03 | 1 | 1 |
| 0.01 | 0.03 | 1 | 2 |
| 0.02 | 0.05 | 1 | 3 |
| … | … | _n_ | _k_ |

I can plot this with an interaction term of _n_ and _k_ for color, like so:

```auto
p1 <- ggplot(results, aes(x, y)) +
    geom_abline(slope = 1, linetype = "dashed", color = "#606060") +
    geom_line(aes(color = interaction(outer_fold:outer_repeat)), size = 1) 

```

Now, finally, for the _averaged_ calibration I just add a loess-smoothed curve:

```auto
p2 <- ggplot(results, aes(x, y)) +
    geom_abline(slope = 1, linetype = "dashed", color = "#606060") +
    geom_line(aes(color = interaction(outer_fold:outer_repeat)), size = 1) +
    geom_smooth(method = "loess")

```

The resulting plots (using a small example dataset) look like this:

 ![calplot](https://discourse.datamethods.org/uploads/default/original/2X/4/492c1729cd586dd1a0824a1e58abf04e5bca3695.png)

My question boils down to whether this is a suitable approach?

My thinking is that with a large enough number of _n_ repeats for the outer loop, just showing the averaged/smoothed curve for each model would suffice. Maybe also include the averaged GLM curve as a reference in order to compare it to ML algorithms. Still, I am very much open to suggestions if there are better approaches!

Finally, a bonus question:

The hyperparameter optimization stage takes a _long_ time (~10+ hours for some models/datasets). What would you deem as an _absolute_ minimum number of repeats of the outer loop? This will help determine if I can run this on a local workstation or if I need HPC resources.

If you’ve read this far, thank you!

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [November 14, 2022, 12:14am UTC](https://discourse.datamethods.org/t/proper-way-to-get-an-averaged-calibration-curve-from-nested-resampling/6207/2 "2022-11-14T00:14:15Z")

</div>

This doesn’t address your questions but I like your approach. Besides averaging _loess_ calibration curves it’s always good to plot all the raw ones going into the average. This is to see if there are any kinds of systematic calibration errors that the averages cover up. But just thinking out loud.