Generalized Competing Event Modeling with gcemod

gcemod

What gcemod does

gcemod estimates covariate effects on omega+ — the ratio of the hazard for an event of interest to the hazard for a competing event — on either the cause-specific (Cox) or subdistribution (Fine-Gray) scale. Confidence intervals and p-values come from the Lunn-McNeil stacked-model construction. Per-subject omega = omega+/(1+omega+) is reported as a descriptive quantity so model-predicted values can be compared with observed ones.

The workflow: (1) fit the GCE model, (2) build a risk score, (3) find cutpoints that maximize the omega+ difference between groups, (4) draw alligator plots of cumulative incidence by group, and (5) check calibration of predicted vs. observed omega+.

Two contrasting cohorts

The package bundles two competing-risks example datasets that illustrate contrasting GCE applications. They are provided for illustration only and do not represent real patients. In hn (head and neck) the event of interest (cancer recurrence) is more frequent than the competing event (omega+ > 1). In prostate the event of interest (distant metastasis or prostate-cancer death) is less frequent than the competing event (other-cause death) (omega+ < 1). Both use status: 0 = censored, 1 = event of interest, 2 = competing event, with time in years.

data(hn)
table(hn$status)
#> 
#>   0   1   2 
#> 366 391 243

1. Fit the GCE model (cause-specific)

Models are specified with a formula whose response is Surv(time, status).

fit <- gcecox(Surv(time, status) ~ age + smoker + t_cat + n_cat + p16,
              data = hn, t = 5)
fit
#> Generalized competing event model [Cox (cause-specific)]
#> Call: gcecox(formula = Surv(time, status) ~ age + smoker + t_cat + 
#>     n_cat + p16, data = hn, t = 5)
#>   n = 1000;  events of interest = 391;  competing events = 243
#> 
#> omega+ ratios, exp(coef):
#>    age smoker t_cat2 t_cat3 t_cat4 n_cat1 n_cat2 n_cat3    p16 
#> 0.9329 0.7179 1.0058 0.8073 0.9071 1.2658 1.3018 0.9140 0.6415 
#> 
#> Use summary() for standard errors, confidence intervals, and tests.

Printing gives a short display. summary() returns the full coefficient table as an object:

s <- summary(fit)
s
#> Generalized competing event model [Cox (cause-specific)]
#> Call: gcecox(formula = Surv(time, status) ~ age + smoker + t_cat + 
#>     n_cat + p16, data = hn, t = 5)
#>   n = 1000;  events of interest = 391;  competing events = 243
#>   at t = 5:  omega+(t) = 1.8031;  omega(t) = 0.6432  (average subject)
#> 
#> Coefficients (log omega+ ratio):
#>             coef exp(coef)  se(coef)      z Pr(>|z|)    
#> age    -0.069406  0.932948  0.010972 -6.326 2.52e-10 ***
#> smoker -0.331362  0.717945  0.192486 -1.721   0.0852 .  
#> t_cat2  0.005787  1.005803  0.459294  0.013   0.9899    
#> t_cat3 -0.214037  0.807318  0.448883 -0.477   0.6335    
#> t_cat4 -0.097552  0.907055  0.447594 -0.218   0.8275    
#> n_cat1  0.235729  1.265831  0.292510  0.806   0.4203    
#> n_cat2  0.263737  1.301785  0.223448  1.180   0.2379    
#> n_cat3 -0.089942  0.913984  0.337527 -0.266   0.7899    
#> p16    -0.443917  0.641519  0.254012 -1.748   0.0805 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Omega+ ratio with 95% confidence interval:
#>        exp(coef) lower .95 upper .95
#> age       0.9329    0.9131    0.9532
#> smoker    0.7179    0.4923    1.0470
#> t_cat2    1.0058    0.4088    2.4744
#> t_cat3    0.8073    0.3349    1.9460
#> t_cat4    0.9071    0.3773    2.1809
#> n_cat1    1.2658    0.7135    2.2457
#> n_cat2    1.3018    0.8401    2.0172
#> n_cat3    0.9140    0.4717    1.7711
#> p16       0.6415    0.3899    1.0554
#> 
#> Omnibus Wald test (all ratios = 1): chi-sq = 49.215 on 9 df, p = <1e-04
#> Standard errors: cluster-robust (subject)

Each p-value tests whether that covariate shifts omega+ (event of interest vs. competing event); the omnibus test asks the same jointly. The usual extractors are available:

coef(fit)
#>          age       smoker       t_cat2       t_cat3       t_cat4       n_cat1 
#> -0.069405699 -0.331362307  0.005786625 -0.214037222 -0.097552275  0.235728618 
#>       n_cat2       n_cat3          p16 
#>  0.263736512 -0.089942229 -0.443917094
confint(fit, exponentiate = TRUE)
#>            2.5 %    97.5 %
#> age    0.9130994 0.9532283
#> smoker 0.4923194 1.0469729
#> t_cat2 0.4088479 2.4743690
#> t_cat3 0.3349304 1.9459649
#> t_cat4 0.3772599 2.1808539
#> n_cat1 0.7134949 2.2457449
#> n_cat2 0.8401181 2.0171505
#> n_cat3 0.4716667 1.7710954
#> p16    0.3899365 1.0554184

The Fine-Gray version uses the same interface:

fit_fg <- gcefg(Surv(time, status) ~ age + smoker + t_cat + n_cat + p16,
                data = hn, t = 5)
summary(fit_fg)

2. Risk score and predictions

The risk score is the linear predictor centered at the covariate means, so the average subject has score 0. predict() returns it, and the predicted omega+ and omega at any time, for the fitting data or for new subjects.

new <- hn[1:3, ]
predict(fit, new, type = "lp")
#>         1         2         3 
#> 0.3111360 1.6653168 0.5548256
predict(fit, new, type = "ratio", times = c(2, 5))   # omega+(t)
#>         t=2      t=5
#> 1  4.033511 2.461148
#> 2 15.624155 9.533470
#> 3  5.146552 3.140298
predict(fit, new, type = "share")                    # omega at t = 5
#>         1         2         3 
#> 0.7110785 0.9050645 0.7584715

3. Cutpoints maximizing the omega+ difference

cuts <- gce_cutpoints(fit, groups = 2, method = "optimal")
cuts
#> GCE risk groups: 2 groups by optimal cutpoint search (criterion: omegaplus)
#> Risk-score cutpoint: -0.02023
#> 
#>  group   n events_of_interest competing_events mean_omegaplus mean_omega
#>    Low 484                181              167          1.106     0.5071
#>   High 516                210               76          3.396     0.7458
#> 
#> Lunn-McNeil omega+ contrast (relative cause-specific hazards):
#>   chi-sq = 31.833 on 1 df, p = <1e-04

4. Alligator plot (cumulative incidence by risk group)

Risk groups are distinguished by color and event types by line type (event of interest solid, competing event dashed) on a single panel. The display functions return objects; plot() draws them.

plot(gce_alligator(fit, groups = cuts))

5. Calibration of predicted vs. observed omega+

cal <- gce_calibration(fit, which = "omegaplus", groups = 5)
cal
#> Calibration of predicted vs. observed omega+ at t = 5 (5 risk-rank groups)
#> 
#>  group   n predicted observed
#>      1 200    0.7158   0.7135
#>      2 207    1.2761   1.7564
#>      3 195    1.8475   3.5769
#>      4 198    2.6248   2.5259
#>      5 200    5.0017   3.0124
plot(cal)

6. Is the relative hazard constant over time?

plot(gce_omegatime(fit, which = "omegaplus", groups = cuts, xmax = 10))

Discrimination

gce_performance(fit)
#> C-index (event of interest): 0.502 (95% CI 0.471-0.532)  [risk-score orientation flipped]

The competing-dominant contrast: prostate

Running the same workflow on prostate shows the opposite regime. Here the competing event (other-cause death) dominates, so omega+ ratios and the risk score point the other way, and the alligator jaws open with the competing-event CIF on top.

data(prostate)
fit_p <- gcecox(Surv(time, status) ~ age + psa + gleason + t2b + comorbidity,
                data = prostate, t = 5)
summary(fit_p)
plot(fit_p, which = "alligator", groups = 2)

The matrix interface

Earlier versions took the follow-up times, a two-column data frame of event indicators, and a covariate data frame. That interface is still available and gives the same fit:

Ind <- data.frame(event     = as.integer(hn$status == 1),
                  competing = as.integer(hn$status == 2))
Cov <- hn[, c("age", "smoker", "t_cat", "n_cat", "p16")]
fit_m <- gcecox(hn$time, Ind, Cov, t = 5)
all.equal(coef(fit_m), coef(fit))
#> [1] TRUE

References

Lunn M, McNeil D (1995) Applying Cox regression to competing risks. Biometrics 51:524–32.

Carmona R, et al. (2014) Validated competing event model for the stage I-II endometrial cancer population. Int J Radiat Oncol Biol Phys 89:888–98.

Carmona R, et al. (2016) Improved method to stratify elderly patients with cancer at risk for competing events. J Clin Oncol 34:1270–77.

Mell LK, et al. (2019) Nomogram to predict the benefit of intensive treatment for locoregionally advanced head and neck cancer. Clin Cancer Res 25:7078–7088.

Zakeri K, et al. (2020) Predictive classifier for intensive treatment of head and neck cancer. Cancer 126:5263–5273.

Mell LK, et al. (2024) Effects of androgen deprivation therapy on prostate cancer outcomes according to competing event risk. Eur Urol 85:373–381.