Steps to construct a nomogram
The following steps show how to build a competing risks regression model and construct a nomogram with it.
1. Fit a competing risks regression model
To start, we first want to create a competing risks regression model.
Our package heavily utilizes the model building procedures and utilities
in the rms package, so the first step is to create a Cox
Proportional-Hazards model for the event of interest (death from
prostate cancer) from the rms::cph function using the
preferred workflow:
# Register the data set
dd <- datadist(prostate.dat)
options(datadist = "dd")
# Fit the Cox-PH model for the event of interest
prostate.f <- cph(Surv(TIME_EVENT,EVENT_DOD == 1) ~ TX + rcs(PSA,3) +
BX_GLSN_CAT + CLIN_STG + rcs(AGE,3) +
RACE_AA, data = prostate.dat,
x = TRUE, y= TRUE, surv=TRUE, time.inc = 144)
prostate.f
#> Cox Proportional Hazards Model
#>
#> cph(formula = Surv(TIME_EVENT, EVENT_DOD == 1) ~ TX + rcs(PSA,
#> 3) + BX_GLSN_CAT + CLIN_STG + rcs(AGE, 3) + RACE_AA, data = prostate.dat,
#> x = TRUE, y = TRUE, surv = TRUE, time.inc = 144)
#>
#> Model Tests Discrimination
#> Indexes
#> Obs 2000 LR chi2 55.16 R2 0.029
#> Events 464 d.f. 11 R2(11,2000)0.022
#> Center -2.7576 Pr(> chi2) 0.0000 R2(11,464)0.091
#> Score chi2 56.50 Dxy 0.230
#> Pr(> chi2) 0.0000
#>
#> Coef S.E. Wald Z Pr(>|Z|)
#> TX=PI 0.5135 0.1431 3.59 0.0003
#> TX=RP 0.0594 0.1246 0.48 0.6336
#> PSA -0.0547 0.0270 -2.03 0.0425
#> PSA' 0.0884 0.0463 1.91 0.0563
#> BX_GLSN_CAT=2 0.1900 0.1137 1.67 0.0948
#> BX_GLSN_CAT=3 0.7213 0.1683 4.28 <0.0001
#> CLIN_STG=T2 -0.2888 0.1116 -2.59 0.0096
#> CLIN_STG=T3 -0.0154 0.2372 -0.06 0.9483
#> AGE -0.0429 0.0127 -3.37 0.0007
#> AGE' 0.0332 0.0158 2.10 0.0353
#> RACE_AA -0.3154 0.1429 -2.21 0.0273Then we can call crr.fit. This function takes the
rms::cph fit and uses the cmprsk::crr function
to convert it to a competing risks regression model, while maintaining
attributes from the original fit.
# Refit to a competing risks regression
prostate.crr <- crr.fit(prostate.f, cencode = 0, failcode = 1)
prostate.crr
#> convergence: TRUE
#> coefficients:
#> TX=PI TX=RP PSA PSA' BX_GLSN_CAT=2
#> 0.30480 0.06632 -0.05174 0.08788 0.12550
#> BX_GLSN_CAT=3 CLIN_STG=T2 CLIN_STG=T3 AGE AGE'
#> 0.63010 -0.26480 0.08370 -0.03062 0.01853
#> RACE_AA
#> -0.25220
#> standard errors:
#> [1] 0.13750 0.11810 0.02638 0.04513 0.11460 0.16120 0.10920 0.21500 0.01240
#> [10] 0.01593 0.14070
#> two-sided p-values:
#> TX=PI TX=RP PSA PSA' BX_GLSN_CAT=2
#> 2.7e-02 5.7e-01 5.0e-02 5.1e-02 2.7e-01
#> BX_GLSN_CAT=3 CLIN_STG=T2 CLIN_STG=T3 AGE AGE'
#> 9.3e-05 1.5e-02 7.0e-01 1.4e-02 2.4e-01
#> RACE_AA
#> 7.3e-022. Summarize and assess model output
We can see that the model object above is of class
cmprsk.
We have created generic methods to assess output:
- The
summary.cmprskmethod utilizessummary.rmsto flexibly summarize model effects and contrasts.
summary(prostate.crr)
#> Effects Response : Surv(TIME_EVENT, EVENT_DOD == 1)
#>
#> Factor Low High Diff. Effect S.E. Lower 0.95
#> PSA 5.000 9.400 4.400 -0.145130 0.074520 -0.291190
#> Hazard Ratio 5.000 9.400 4.400 0.864910 NA 0.747380
#> AGE 58.137 69.562 11.425 -0.184680 0.077432 -0.336440
#> Hazard Ratio 58.137 69.562 11.425 0.831370 NA 0.714310
#> RACE_AA 0.000 1.000 1.000 -0.252210 0.140700 -0.527970
#> Hazard Ratio 0.000 1.000 1.000 0.777080 NA 0.589800
#> TX - EBRT:RP 3.000 1.000 NA -0.066323 0.118050 -0.297700
#> Hazard Ratio 3.000 1.000 NA 0.935830 NA 0.742530
#> TX - PI:RP 3.000 2.000 NA 0.238490 0.132460 -0.021125
#> Hazard Ratio 3.000 2.000 NA 1.269300 NA 0.979100
#> BX_GLSN_CAT - 2:1 1.000 2.000 NA 0.125530 0.114610 -0.099103
#> Hazard Ratio 1.000 2.000 NA 1.133800 NA 0.905650
#> BX_GLSN_CAT - 3:1 1.000 3.000 NA 0.630090 0.161180 0.314180
#> Hazard Ratio 1.000 3.000 NA 1.877800 NA 1.369100
#> CLIN_STG - T2:T1 1.000 2.000 NA -0.264760 0.109220 -0.478820
#> Hazard Ratio 1.000 2.000 NA 0.767390 NA 0.619510
#> CLIN_STG - T3:T1 1.000 3.000 NA 0.083704 0.214990 -0.337680
#> Hazard Ratio 1.000 3.000 NA 1.087300 NA 0.713430
#> Upper 0.95
#> 0.00092617
#> 1.00090000
#> -0.03291100
#> 0.96762000
#> 0.02355300
#> 1.02380000
#> 0.16505000
#> 1.17950000
#> 0.49809000
#> 1.64560000
#> 0.35016000
#> 1.41930000
#> 0.94600000
#> 2.57540000
#> -0.05069400
#> 0.95057000
#> 0.50509000
#> 1.65710000- The
anova.cmprskmethod utilizesanova.rmsto test for statistical significance of model terms.
anova(prostate.crr)
#> Wald Statistics Response: Surv(TIME_EVENT, EVENT_DOD == 1)
#>
#> Factor Chi-Square d.f. P
#> TX 5.21 2 0.0739
#> PSA 3.85 2 0.1458
#> Nonlinear 3.79 1 0.0515
#> BX_GLSN_CAT 15.29 2 0.0005
#> CLIN_STG 6.88 2 0.0320
#> AGE 9.27 2 0.0097
#> Nonlinear 1.35 1 0.2445
#> RACE_AA 3.21 1 0.0730
#> TOTAL NONLINEAR 5.16 2 0.0758
#> TOTAL 44.64 11 <.0001Obviously there is a lot of work that would go into determining that
we have the right model specifications, covariates, term specifications,
etc. We’ll assume that the current model is the one we are moving
forward with. See the cmprsk package to learn more about
the competing risks modeling framework used here, how to interpret
parameters, etc.
3. Validate model predictions
Suppose we are particularly interested in accurately predicting the 10-year risk of death from prostate cancer.
Generate cross-validated predictions
First, we can use the tenf.crr function to generate
K-fold cross-validated predictions for each observation in the
development data set. The default is to set fold=10 for
10-fold cross-validation, so we’ll leave it at that:
set.seed(123)
prostate.dat$preds.tenf <- tenf.crr(prostate.crr, time = time_of_interest)
#> 1 2 3 4 5 6 7 8 9 10
str(prostate.dat$preds.tenf)
#> num [1:2000] 0.374 0.376 0.277 0.372 0.394 ...This vector shows us, for each patient, the estimated (out of sample) probability of death from prostate cancer given the covariates in the model, adjusting for the possibility of death from other causes as a competing event.
Compute the concordance index (C-Index)
The cindex function computes the concordance
index for binary, time-to-event, or competing risks outcomes. In
this case, we are interested in the competing risks version, which
utilizes the event times, and only considers pairwise comparisons
involving those with the event of interest when testing for
concordance.
For example, if Patient A experienced the event of interest at 5 years and Patient B experienced the event of interest at 10 years, the algorithm checks to see if the predicted risk for Patient A is larger than that of Patient B–in which case they would be considered concordant. There are similar comparisons if one of the patients had a different event status. However, if Patient A nor Patient B experienced the event of interest (i.e., they either experienced a competing event or were censored), they would not be compared.
cindex(
prob = prostate.dat$preds.tenf,
fstatus = prostate.dat$EVENT_DOD,
ftime = prostate.dat$TIME_EVENT,
type = "crr",
failcode = 1
)
#> N n usable concordant cindex
#> 2.000000e+03 2.000000e+03 5.652590e+05 3.228440e+05 5.711435e-01The output displays the number of observations (in the original input
N, and those used n), the number of pairwise
comparisons made (usable), the number of concordant pairs
(concordant), and the concordance index
(cindex), which is the proportion of the usable pairs that
were concordant. The C-index ranges from 0.50, a coin flip, to 1,
perfect discrimination. In this case, our model has relatively poor
discrimination ability.
Assess model calibration
We can also assess how calibrated our model predictions are
by comparing the predicted probabilities from the model to the actual
event rate observed. The groupci function does this by
breaking up the input risk distribution into groups (e.g., by quantiles,
user-defined groups, etc.) and computing the observed cumulative
incidence within each group at the time point in which the predictions
were made for. The goal is for these quantities to follow 45-degree line
such that the risks produced by the models align with the true rates at
which events occur.
groupci(
x = prostate.dat$preds.tenf,
ftime = prostate.dat$TIME_EVENT,
fstatus = prostate.dat$EVENT_DOD,
g = 10, # Deciles
u = time_of_interest,
failcode = 1,
xlab = "Predicted 10-year prostate cancer-specific mortality",
ylab = "Actual 10-year prostate cancer-specific mortality"
)
#> x n events ci std.err
#> [1,] 0.2195250 200 41 0.2137948 0.03444284
#> [2,] 0.2614243 200 42 0.2594906 0.04058236
#> [3,] 0.2844747 200 45 0.2740857 0.04014348
#> [4,] 0.3045506 200 48 0.3427236 0.04718743
#> [5,] 0.3239753 200 50 0.3807385 0.04866854
#> [6,] 0.3418128 200 43 0.2773428 0.04117005
#> [7,] 0.3595614 200 48 0.3746006 0.04898192
#> [8,] 0.3808118 200 44 0.4879819 0.06502863
#> [9,] 0.4101135 200 49 0.3324923 0.04835294
#> [10,] 0.4827395 200 54 0.4596141 0.05677256There might be some evidence that our model slightly overestimates the actual risk of 10-year prostate cancer-specific mortality for those most at risk.
4. Create the nomogram
Once the model is built, validated, and ready to use in practice, we
can construct the nomogram. The nomogram.crr function takes
our original model fit and maps it to a diagram that allows it to be
printed and used in a clinical (or any other) setting to quickly
“plug-in” model covariates and generate a risk estimate.
# Set some nice display labels (also see ?Newlevels)
prostate.g <-
Newlabels(
fit = prostate.crr,
labels =
c(
TX = "Treatment options",
PSA = "PSA (ng/mL)",
BX_GLSN_CAT = "Biopsy Gleason Score Sum",
CLIN_STG = "Clinical Stage",
AGE = "Age (Years)",
RACE_AA = "Race"
)
)
# Construct the nomogram
nomogram.crr(
fit = prostate.g,
failtime = time_of_interest,
lp = FALSE,
xfrac = 0.65,
fun.at = seq(0.2, 0.45, 0.05),
funlabel = "Predicted 10-year risk"
)
It is a point system-based diagram in which we can extract the individual contribution of each risk factor (top axis), then add them up to obtain the final risk estimate (using the “Total Points” axis, then drawing a line straight down to the risk axis). Also notice that even with non-linear terms in the original model, they are able to be displayed with axes adjusted accordingly. Interactions can be handled as well.