To make a calibration plot for survival probabilities estimated from a Cox model, one can divide the estimated risk into groups, calculate the average risk within a group, and then compare this to the Kaplan-Meier estimate. What are alternative approaches that do not require binning? What are the specific steps required to implement such a technique, and what is the logic behind it? Some discussion is present in this question.
A straightforward way to calibrate Cox survival models is to use the
calibrate function provided by the
rms package in R, as in the page that you linked. This package provides a
cph method for Cox models that is designed to work with the calibration and validation methods that it provides for several types of regression models.
Quoting from the manual page: "[
calibrate u]ses bootstrapping or cross-validation to get bias-corrected (overfitting- corrected) estimates of predicted vs. observed values." Instead of using binned Kaplan-Meier estimates to provide "observed" values, however, the default
cmethod="hare" setting in
cph models uses regression-spline interpolation estimates that allow for non-proportional hazards and nonlinearity while take censoring into account. This adaptive modeling of the observed data allows for a continuous calibration plot for a specific survival time.
This outline by Kooperberg et al in the Encyclopedia of Statistical Sciences provides more detail on the HARE method; if you can't get access to that, there's a proof version currently available here. The original paper describing the method is Kooperberg et al., J Am Stat Assoc 90: 78-94, 1995. (If you set
cmethod="KM" in the call to
calibrate you will get comparison against binned Kaplan-Meier estimates.)
Before you do this, make sure to read the manual pages to ensure that you generate the
cph model with the correct settings that allow use of these functions, and that the
polspline package is installed to provide the
Frank Harrell, the author of the
rms package, compares the binned Kaplan-Meier and continuous HARE approach on pages 506-9 of Regression Modeling Stratgies, second edition.
What's going on "under the hood"
The calibration problem for survival analysis is that observations are events while we are trying to calibrate probabilities of the events. So some type of interpolation is needed. For Cox models, a particular survival time is taken for analysis.
Consider the binned KM calibration. You start by grouping cases together by predicted probabilities, then for each group of similar predicted probabilities you plot the KM survival curve and interpolate among those few cases at the chosen survival time to get the estimated "observed" survival probability for that group.
In the approach used by
cph models, you instead first interpolate the hazard among all the cases as a general function of the predictor variables and time. This general function, as provided by
hare, allows each predictor variable and time to be modeled as a linear spline, and allows for pairwise interactions among the variables (including time as a variable). The collection of splines and interactions provides a set of basis functions each of which is a function of predictor variables and time. You then look for the combination of these basis functions that best fit the survival data, with the coefficients weighting the basis functions estimated by maximizing partial likelihood, similarly to how a Cox model finds coefficients for the predictor variables themselves. Model complexity (e.g., which basis functions to include, number of spline knots) is selected by a stepwise addition and deletion process that is nicely explained in the outline noted above.
The result of the
hare process is a single function providing the hazard as a function of all the predictor variables and time, in a form that will typically be much more complicated than the proportional-hazards Cox model. I like to think of this as starting with an interpolation among all the cases at all times, along with their associated values of predictor variables, rather than the interpolation to a particular time for a subset of cases used in the binned KM calibration scheme.
With the hazard function provided by
hare, you can now take any combination of predictor variable values, and compare the "predicted" Cox probability of survival at a particular time against the interpolated estimate of "observed" survival probability provided by
hare. (Note that, in principle, you are not limited to the cases at hand for this comparison; any reasonable combination of predictor variable values can be examined.) The calibration curve is then a smoothed plot of "observed" versus "predicted" survival probability, at the desired time, among all the cases. The
calibrate method in
rms repeats this process for multiple bootstrapped samples from the cases to gauge how well the results will generalize to the population from which the cases were drawn.
If you want not only to look under the hood but also to deconstruct the engine, that's easy in R. When the
rms package is loaded, type
rms:::calibrate.cph at the command prompt to get the code for this
calibrate method. The wrapper for
hare is available by typing
hare at the prompt when the
polspline package is loaded; much of the work is done by compiled functions whose source code is available from CRAN.