Create calibration plot in R with vectors of predicted and observed values

Question

I am currently working on a project regarding the external validation of a logistic regression model for binary classification. I would like to create a calibration plot and compute the calibration-in-the-large (intercept) and calibration slope, like the following figure.

I noted that most online tutorials involved using the lrm object in R to compute the calibration-in-the-large and calibration slope. However, no lrm object is available for my project as the logistic regression results are calculated elsewhere. Instead, I have two vectors, one contains all predicted values (e.g. 0.1, 0.5, 0.8, 0.1 ...), the other contains observed values (e.g. 0, 1, 1, 0 ...).

I am currently plotting the calibration curve manually by breaking all observations into quantiles and calculating the mean observed risk and the mean predicted risk in each quantile. I wonder if there is a way to get the calibration-in-the-large and calibration slope by using the two vectors only. It would be best if there is a way to create a sophisticated calibration plot (like the attached one with zone of confidence interval ) directly from two vectors without any manual calculation.

Thanks a lot!

Answer 1

It sounds like what you want to do is check if, when an event probability of $p$ is predicted by the model, the event indeed happens with probability $p$. After all, there is a sense in which the model is lying if it keeps predicting events to have probability of $0.5$ that almost never (or almost always) happen. Likewise, low and high predicted probabilities should correspond to events that rarely and frequently occur, respectively.

An R function to do this is rms::val.prob from Frank Harrell's rms package. If the claimed probabilities are reflective of the reality of how often the event occurs, that will correspond to a 45-degree line.

library(rms)
set.seed(2023)
N <- 1000
probability <- runif(N, 0, 1)
event <- rbinom(N, 1, probability)
rms::val.prob(probability, event)

By the construction of this simulation, the fit is almost the ideal. However, funky probability values result in a grossly non-ideal fit (though both have the same $c$-index, equal to the area under the ROC curve, since I applied a strictly increasing transformation to the probabilities).

library(rms)
set.seed(2023)
N <- 1000
probability <- runif(N, 0, 1)
event <- rbinom(N, 1, probability)
rms::val.prob(probability^2, event)

An advantage of doing the probability validation with rms::val.prob instead of using quantiles is that you don't have to rely on somewhat arbitrary binning; you smoothly evaluate all data all at once.

A question of mine asks how damaging to the analysis it would be to run val.prob when calibration is the truly correct action. My sentiment from the answer by EdM is that, while val.prob is not ideal, the damage is not so enormous.