Walk through rms::calibrate for logistic regression

The calibrate function in the rms R package allows us to compare the probability values predicted by a logistic regression model to the true probability values. This is easy enough: just plot them and make sure they are about the line $$y=x$$. However, if we knew the true probability values, there would not be any need to do statistical inference! We have to estimate what the true values are, which is a major part of what rms::calibrate does and what confuses me.

I have figured out a bit of how rms::calibrate works.

1. It estimates the true probability values using the observed $$0/1$$ labels. This makes sense. If we have a lot of $$1$$s, we should have a lot of high values of $$P(Y=1)$$.

2. It takes bootstrap samples of the predicted probability values paired with the observed outcomes.

3. It fits some kind of model with the bootstrap samples of predicted probability values paired with observed outcomes.

Step $$3$$ is where I start to get lost. What kind of model is fit? Is it, indeed, fit to the bootstrap samples of predicted probability values and observed outomes, or does it fit a new logistic regression each time? (I don't think this is how it work, but I want to make sure.) How does it go from multiple models of bootstrap data to the "apparent" and "bias-corrected" values displayed on the graph?

EXAMPLE

library(rms)
set.seed(2022)
N <- 1000
x <- runif(N, -3, 3)
z <- x
pr <- 1/(1 + exp(-z))
y <- rbinom(N, 1, pr)
plot(x, y)
L <- rms::lrm(y ~ x, x = T, y = T)
cal <- rms::calibrate(L)
plot(cal)


Both calibrate() and validate() in rms by default use the optimism bootstrap, explained for example in Chapter 7 of Elements of Statistical Learning. In your example, you do repeat the logistic-regression modeling on each of multiple bootstrapped data samples.

In general with an optimism bootstrap, you develop a new model on a bootstrap sample. You find how much better that model performs (by some measure) on its bootstrap sample than it does on the full data set. That's the "optimism bias" of that model for that measure.

By the bootstrap principle, an average "optimism" over multiple bootstrap-based models estimates the optimism of the original full model when applied to the population from which the original data were sampled. You then can correct the original model by that average optimism.

For calibrate(), the "optimism" is in the calibration error as a function of predicted probability values. You first calculate the calibration error of the full model over a range of predicted probability values. You then estimate the "optimism" in those calibration error estimates by the optimism bootstrap, and subtract the averaged "optimism" from the original calibration of the full model.

The plot then shows smoothed curves for the original calibration and the optimism-corrected calibration. In the plot you show, the original model itself ("Apparent") deviates some from the "Ideal" calibration, but there doesn't seem to be much optimism (difference between "Apparent" and "Bias-corrected") in the calibration.

You can see details in the code; for logistic regression (lrm objects) type rms:::calibrate.default at the R command prompt. That includes a cal.error() internal function for evaluating calibration, and a call to the predab.resample() function for the optimism bootstrap.

• So if I have probability values and outcomes, but I do not (and will not) have access to the model that created the probability values, I can’t replicate rms::calibrate? Among other ideas, I was hoping to check the calibration of NFL probability of spectacular plays, where it sure seems like $1\%$ chance of a spectacular play results in a successful spectacular play a lot more often than $1\%$ of the time, but I won’t get access to the model giving such probability values. All I would have are the predicted probability values and play outcomes.
– Dave
Feb 2 at 11:12
• @Dave for that application, if there are enough predictions of a "1% chance of a successful spectacular play," you could evaluate whether the observed and predicted frequencies agree to within the error expected from binomial sampling. With 120 to 130 plays typical per NFL game, you would expect at least 1 such play per game. And you'll remember those plays more than the boring ones.
– EdM
Feb 5 at 16:16
• I think I might be more interested in rms::val.prob for the NFL question, agreed? That function seems related to rms::calibrate, but I am having issues following that one, too.
– Dave
Feb 10 at 14:38