I am interested in doing this to visualise the degree of spread about the regression line... This would allow me to better assess the model for predictive accuracy.
To answer the question you posed, a model that inherits from the basic coxph
object of the R survival
package, like a cph
object from the rms
package, contains linear.predictors
for the cases (log hazard relative to some particular choice of reference case). So you can simply extract those values for plotting.
That's not, however, the best way to "assess the model for predictive accuracy." If you do that for a single predictor like age, particularly in a model with 11 predictors some of which are likely to be correlated with age, you won't know how much the other predictors are posing problems. It's also not clear how that would deal with censored survival time properly.
Furthermore, the interest in a survival model is typically in the ability to discriminate between individuals in terms of survival order or the precision of predictions at a specific time of interest. Remember that the log-hazard is with respect to a reference baseline hazard, which typically is associated with a wide distribution of survival times even at the reference condition.
With only 225 events, the best way to "assess the model for predictive accuracy" overall is to validate and calibrate the modeling process by resampling from your data, repeating the modeling, and comparing predictions on the resampled data and the full data set. As you are already using the rms
package you can use the validate()
function to evaluate overfitting/optimism and measures of the model's ability to discriminate among cases. The calibrate()
function documents how well predicted and observed probabilities of survival agree at any time point of interest.
If you are particularly interested in the quality of the fit for a continuous predictor like age, the best choice is to evaluate the martingale residuals, as outlined in this answer and explained in detail in the section on assessment of Cox-model fits in Frank Harrell's online Regression Modeling Strategies book. As Therneau and Grambsch explain in Section 4.2.1, they bear some similarity to ordinary residuals in least-square models:
The martingale residual is really a simple difference $O - E$ between the observed number of events for an individual and the conditionally expected number given the fitted model, followup time, and the observed course of any time-varying covariates.
You can't use martingale residuals exactly like ordinary residuals of linear regression, but they can be useful for evaluating how well your model has captured the shape of the association between a continuous predictor and outcome. They allow you to take all the other predictors into account. If a smoothed plot of martingale residuals from your model against the values of the predictor is reasonably flat, then you've done OK.
If not, you can repeat the model with the particular predictor omitted, then produce a smoothed plot of martingale residuals against the values of the omitted predictor. The shape of that plot can indicate the shape of the association. Better, you can let the data tell you the shape of the association by fitting the predictor flexibly, as with a regression spline.