I am hoping to get advice on how to construct a confidence score for a Cox PH survival model. I used an electronic medical record (EMR) dataset to train a survival model. For each patient, a time-weighted average of all past clinical data is used to predict survival time. Now, for a patient coming into clinic, I would like to display to the doctor two things:

- a predicted survival curve
- a confidence score for this patient, indicating how likely is the true survival distribution to be close to the predicted one.

I am aware that bootstrapping can be used to generate confidence intervals for predictive models, as in this link:

But for my use case, I suspect that using data from outside the Cox model could improve the prediction of confidence. With EMR data, some patients will have large amounts of missing data. So the model performance will vary depending on how much data is available for each patient. For example, take *n _{lab}*, the number of laboratory results in the EMR for this patient. Patients with high

*n*will have improved predictive performance, but if

_{lab}*n*were included in the Cox model, its coefficient may be close to 0 and its value may not affect the width of the bootstrapped confidence intervals for the predicted survival curve.

_{lab}One solution could be to fit a second model, which predicts the Cox modelâ€™s squared deviance residual, or maybe log likelihood, for each observation. When the second model predicts a large residual or low log likelihood, we would say that the Cox modelâ€™s confidence score is low. But, this seems hacky so I would love to hear if there is a more theoretically justified solution.

Here is some sample code which illustrates the concept of the second model that predicts the performance of the Cox model.

```
library(survival)
library(rms)
data(kidney)
cph.model <- cph(Surv(time, status) ~ frail,
data=kidney,surv=TRUE,x=TRUE,y=TRUE)
kidney$resid_squared <- residuals(cph.model, type='deviance')**2
conf.model <- lm(resid_squared ~ frail+age, data=kidney)
conf.scores <- predict(conf.model, kidney)
```