4
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Your last line of code leads to the "singular" X matrix because the PI value is a linear combination of the 3 other included predictors. Try running that code without the PI value to see what your estimates of the coefficients might be, independent of what you found in the literature. With only 17 events your estimates will tend to be imprecise; the number of events, not the total number of cases, matters for precision of coefficient estimates. $\endgroup$
    – EdM
    Commented Dec 9, 2020 at 18:19
  • 1
    $\begingroup$ After looking at the Royston and Altman reference, I see that they intend the PI in your last model to be an offset, with coefficient constrained to be exactly 1. That's done with a term offset(PI) instead of just PI in your last model. Please try that and edit your question accordingly. $\endgroup$
    – EdM
    Commented Dec 10, 2020 at 15:29
  • $\begingroup$ Thanks, I was assuming that this was the reason for the error too. Regarding the number of events: I used only parts of my dataset to illustrate my problem, the full dataset fortunately has many more outcomes. I edited the question accordingly. $\endgroup$
    – b_surial
    Commented Dec 10, 2020 at 16:49
  • $\begingroup$ Sorry I didn't see your comment above before posting my response. I updated the question. Maybe you can post an answer that I can accept? You already helped me a lot. $\endgroup$
    – b_surial
    Commented Dec 10, 2020 at 17:05

2 Answers 2

1
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.