# External validation of a published Cox PH model

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.

• 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.
– EdM
Commented Dec 9, 2020 at 18:19
• 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.
– EdM
Commented Dec 10, 2020 at 15:29
• 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. Commented Dec 10, 2020 at 16:49
• 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. Commented Dec 10, 2020 at 17:05

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.
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.