Is it possible to generate a scatterplot of relative hazard for participants superimposed on Cox regression line?

I have been using the rms package in R to perform Cox regression on time-to-event data. I have used the plot(Predict()) function to create plots of regression lines in the log relative hazard and relative hazard scale against continuous covariates.

An example plot (not my own) in the log relative hazard scale is shown below: (source https://thomaselove.github.io/432-notes/cox-regression-models-for-survival-data-example-2.html).

Is it possible to superimpose scatter points of the relative hazard (compared to the reference defined in the Predict function) of every individual patient in the data frame? I am interested in doing this to visualise the degree of spread about the regression line.

Ideally, I would try and do this both in the log relative hazard and relative hazard scale. This would allow me to better assess the model for predictive accuracy. An analogous plot (not my own) for multivariable linear regression is shown below: (source https://cran.r-project.org/web/packages/ggiraphExtra/vignettes/ggPredict.html)

Many thanks!

Addendum: Thank you EdM for your kind reply. In total I have 11 predictor variables in the model with 3000 participants and 225 occurrences of the primary outcome

• Is age the only predictor in your model? If not, how many others are there and how many events are in your data set? Please provide that information by editing the question, as comments are easy to overlook and can be deleted.
– EdM
Commented Dec 26, 2022 at 17:40

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.