# How to estimate a calibration curve with bootstrap (R)

Question: I have fitted a probabilistic model (bayesian network) for modeling a binary outcome variable. I would like to create a high-resolution calibration plot (e.g. spline) corrected for overfitting with bootstrapping. Is there a standard procedure for calculating such a curve?

Considerations: I could do this easily with train/test splitting, but I would rather not throw away any data as I have less than 20,000 samples. So I naturally thought about bootstrapping. I know that one such function (calibrate) is implemented in Frank Harrell's rms package, but unfortunately the model I use is not supported by the package.

Bonus question: is it possible to recalibrate a miscalibrated model with bootstrapping? The reason I ask this is that I tried to recalibrate a model by

1. split data in train/test
2. fitting model to train set
3. recalibrate model to train set (with a cubic spline)
4. evaluate calibration on test set

The models recalibrated in the fashion above were perfectly calibrated on the train set but not so much on the test set, which probably indicates mild overfitting. I also tried further splitting the test set, calibrating on one split and evaluating the calibration on the second split. I got better results (still not perfectly calibrated though), but the sets became quite small (~1000 samples) and thus the calibration unreliable

• Not sure what's available in R, but Python's sklearn has an excellent probability calibration module which supports cross-validation and isotonic/monotonic regression, which are key for high quality probability calibration.Might be a good place to get some ideas. scikit-learn.org/stable/modules/calibration.html – olooney May 17 '18 at 16:23
• Thank you for your answer! I am bound to R and not particularly keen on starting with reticulate. I could implement the procedure myself if I knew it, but I haven't found it anywhere.. I guess I am just hoping that prof Harrell will see this question :D – Gino_JrDataScientist May 17 '18 at 20:42
• Please include a (small) data example to illustrate what you're talking about. Thank you. – Jim May 24 '18 at 21:57
• You should consider taking a look at this introduction to the rms package (and its calibrate function) in R: r-bloggers.com/introduction-to-the-rms-package – rpatel May 25 '18 at 15:57
• Hi rpatel, thanks for the suggestion. I had mentioned the rms::calibrate function in my original question, noting that it does not support the model class I am using. I also own Harrell's Regression Modelling Strategies book, but I am unable to find any detailed description on how the calibrate function works. – Gino_JrDataScientist May 28 '18 at 7:58

After discussing with prof Frank Harrell by email, I devised the following procedure for estimating the optimism-corrected calibration curve, partially based on his Tutorial in Biostatistics (STATISTICS IN MEDICINE, VOL. 15,361-387 (1996)):

1. fit a risk prediction model on all data
2. fit a flexible model (gam with spline and logit link) to the model's predicted probabilities vs outcome, and query the gam at a grid of predicted probabilities $$p=(0.01,0.02,...,0.99)$$. This is the apparent calibration curve and we call it $$cal_{app}$$
3. draw bootstrap sample with replacement, same size of original data
4. fit risk prediction model on bootstrap sample
5. use the bootstrap model to predict probabilities from the bootstrap sample, fit a gam between the predicted probabilities and the outcome, and query the gam at a grid of predicted probabilities (let us call these points $$cal_{boot}$$)
6. use the bootstrap model to predict probabilities from the original sample, fit a gam between the predicted probabilities and the outcome, and query the gam at a grid of predicted probabilities obtaining a calibration curve ($$cal_{orig}$$)
7. compute the optimism at every point $$p$$ of the grid like so $$Optimism(p)=cal_{boot}(p) - cal_{orig}(p)$$
8. repeat steps 3-7 some 100 times, average the optimism at each point $$p$$
9. compute the optimism corrected calibration like so $$cal_{corr}(p)=cal_{app}(p)-$$

Important note: The procedure above is inspired by Harrell's work and my discussion with him, but all errors are mine alone.