How to make calibration plot for survival data without binning data?

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 calibrate for 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 provides more detail on the HARE method, with a reference to the original paper. (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 hare functionality.

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, and (more cryptically) in his associated course notes, chapter 17, pages 18-19.

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 calibrate for 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.

• Thanks. If possible, could you elaborate on the logic of using the 'hare' method. I see the benefits listed, but I don't understand what is going on 'under the hood'. – julieth Apr 8 '16 at 22:58
• @julieth, I've tried to open the hood for you in an addition to my answer. – EdM Apr 10 '16 at 15:44