External validation of a published Cox PH model

Question

My aim is to externally validate a risk prediction model published in the medical literature that is based on a Cox regression model. I have a dataset with all the variables from the score. I read Prof. Harrell's RMS course notes and am somehow familiar with the rms package in R. I have also searched the resources at https://discourse.datamethods.org/ but I couldn't find an answer on how to proceed.

As far as I understand the documentation of val.surv, I need a fit object generated with cph which I don't have access to so I cannot use that function for external validation. I have read the article Royston P, Altman DG. External validation of a Cox prognostic model: principles and methods. BMC Med Res Methodol. 2013 Mar 6;13:33. which is great but I need some more guidance on how to perform the steps proposed in R.

The first step seems pretty clear. The authors write: A recognised approach to validation is to estimate the regression coefficient on the PI or risk score in the validation dataset [4,12], sometimes known as the ‘calibration slope’.

I calculate a linear predictor (=PI) using the published coefficients beta1-beta3 as follows (please note: I used sample data, not the larger full dataset with many more events for all the code displayed here) :

PI <- x1 * beta1 + x2 * beta2 + x3 * beta3

Then I regress on the PI using a Cox model:

S <- Surv(event = df$event, time = df$time)
cph(S ~ PI, data = df)

## Sample output
Cox Proportional Hazards Model
 
 cph(formula = S ~ PI, data = df)
 
                        Model Tests    Discrimination    
                                              Indexes    
 Obs       950    LR chi2      8.76    R2       0.045    
 Events     17    d.f.            1    Dxy      0.383    
 Center 2.5133    Pr(> chi2) 0.0031    g        0.921    
                  Score chi2   8.00    gr       2.512    
                  Pr(> chi2) 0.0047                      
 
    Coef   S.E.   Wald Z Pr(>|Z|)
 PI 0.6259 0.2264 2.76   0.0057  

I then use the coefficient of PI (in this case 0.6259) to interpret if <1 meaning discrimination is poorer, and > 1 meaning discrimination is better in the validation compared to the derivation set, correct?

And do the R2 and Dxy (and c-statistic calculated as $c = Dxy/2 + 0.5$) I get from this regression refer to the measures of discrimination I would want to interpret for the model performance of my data?

For the second step, authors say *"One reason why the slope on the PI may differ from 1 in the validation dataset is that the regression coefficients for one or more covariates may differ between the datasets. This can be tested formally (ignoring uncertainty of estimates in the derivation dataset) by running a Cox regression on the covariates x in the validation dataset, ‘offsetting’ the original PI evaluated in the validation dataset. [...] From the point of view of successful validation, the ‘best’ result is that all the coefficients $β$ are 0."

How do I do this? One approach I tried gave me an error:

cph(S ~ x1 + x2 + x3 + PI, data = df)
# Gives the error: X matrix deemed to be singular; variable PI

Any guidance on this is highly appreciated.

Update:

Following the suggestion from EdM below, I refitted the model with PI as offset term:

a <- cph(S ~ x1 + x2 + x3 + offset(IP), data = df)

print(a)
## Output:
Cox Proportional Hazards Model
 
 cph(formula = S ~ x1 + x2 + x3 + offset(IP), data = df)
 
                         Model Tests    Discrimination    
                                               Indexes    
 Obs        950    LR chi2      4.44    R2       0.023    
 Events      17    d.f.            4    Dxy      0.390    
 Center -0.6589    Pr(> chi2) 0.3494    g        0.940    
                   Score chi2   6.89    gr       2.560    
                   Pr(> chi2) 0.1417                      
 
                 Coef    S.E.   Wald Z Pr(>|Z|)
 x1              -0.0256 0.0239 -1.07  0.2839  
 x2               1.4074 0.7718  1.82  0.0682  
 x3=level1        0.1526 0.7678  0.20  0.8424  
 x3=level2       -0.1006 0.5923 -0.17  0.8651  

In this hypothetical example I see that x2 is rather far away from 0, meaning that there might be some problem there, is that correct?

Update 2:

Thanks to EdM's useful explanation, I think I have understand how to proceed. As indicated by EdM and Royston and Altman, they recommend performing a joint test for the coefficients first, which would correspond to the TOTAL row in the anova.rms function output:

anova(a)

# Output:
                Wald Statistics          Response: S 

 Factor     Chi-Square d.f. P     
 x1         1.15       1    0.2839
 x2         3.33       1    0.0682
 x3         0.15       2    0.9269
 TOTAL      5.77       4    0.2169

As EdM pointed out, the sample dataset I chose was probably not great since it contains only very few events, and therefore statistical power is very low.

Nevertheless, I hope this post and its answer will be useful for others too.

Answer 1

I have written a blog post with R code about the steps in the paper: Royston P, Altman DG. External validation of a Cox prognostic model: principles and methods. BMC Med Res Methodol. 2013 Mar 6;13:33 Here is the link

Answer 2

What you've done is correct, although there are some cautions to consider.

The first regression of your survival data against the corresponding linear predictor values, calculated from the coefficients reported in the literature and your covariate values, evaluates 2 things in particular.

One is the ability of the model to discriminate the order of events among individuals in your cohort. The C-index, as you calculate it from the Dxy value, is 0.69. That means that the order is correct for 69% of pairs of individuals for which the comparison can be made. As Royston and Altman caution, however, the c-index "is known to be be biased away from the null due to right-censoring of times to event," and you do have a very large fraction of censored cases. The very low $R^2$ also raises red flags about discrimination.

The second is the potential optimism in the reported model when applied to your data. That's evaluated in terms of the slope coefficient in that first regression. The coefficient is only 0.63, well below the value of 1 that would be expected if the model worked as well on your data as it did in the original report. I would be very worried if I got that low a slope during internal resampling validation of a model; that would suggest a highly overfit model that wouldn't generalize well to new data. In fairness, you probably can't reliably distinguish that slope in your regression example from a value of 1; the small number of events leads to a large standard error.

Your second regression, incorporating that linear predictor as an offset in the model along with your own covariate values, helps identify a potential difference between your data and those in the original literature report. A coefficient substantially different from 0 in that regression indicates a different association between that predictor and outcome from what was seen in the original report. In your case, x2 raises questions as you note. Again, with the small number of events in your example, that coefficient is not "significantly" different from 0. Furthermore, Royston and Altman suggest first performing a joint test that all covariate coefficients are 0, in this model with an offset for the linear predictor, to minimize false-positives from multiple comparisons.