Package {gcemod}

gcemod: Generalized Competing Event Models with Lunn-McNeil Testing

Description

The gcemod package fits generalized competing event (GCE) models and estimates covariate effects on \omega^{+}, the ratio of the hazard for an event of interest to the hazard for a competing event, on both the cause-specific (gcecox) and subdistribution (gcefg) hazard scales.

Details

\omega^{+} = h_{1}(t) / h_{2}(t)

is the hazard for the event of interest relative to the hazard for the competing event. Per-subject \omega = \omega^{+} / (1 + \omega^{+}) is retained as a descriptive quantity (the fraction of hazard attributable to the event of interest) for comparison of model-predicted with observed values; covariate effects are modeled on \omega^{+} only.

Confidence intervals and p-values come from the Lunn-McNeil (1995) stacked/augmented data construction: a single model on subject records duplicated by event type and stratified by event type, whose covariate-by-event-type interaction estimates the log ratio of \omega^{+}. See lunnmcneil.

Key functions

gcecox

GCE model via cause-specific Cox regression.

gcefg

GCE model via Fine-Gray regression.

gcecif

GCE model on the cumulative-incidence scale: covariate effects on the ratio \rho(t) = F_1(t)/F_2(t) by pseudo-observation regression.

lunnmcneil

Stacked-model estimator and tests.

gce_cutpoints

Risk-score cutpoints maximizing omega+ separation.

gce_alligator

Cumulative incidence by risk group.

gce_calibration

Predicted vs. observed omega+ by rank.

gce_performance

Cause-specific concordance (C-index).

gce_scaling

Covariate means/SDs behind the risk score.

gce_riskscore

Score new subjects vs. the average patient.

gce_tidy

Tidy coefficient tables, including the primary/competing/total/relative effects comparison (effects = TRUE).

Author(s)

Maintainer: Loren Mell lmell@ucsd.edu

Authors:

• Loren Mell lmell@ucsd.edu

References

Lunn M, McNeil D (1995) Applying Cox regression to competing risks. Biometrics 51:524-32. doi:10.2307/2532940

Carmona R, Gulaya S, Murphy JD, et al. (2014) Validated competing event model for the stage I-II endometrial cancer population. Int J Radiat Oncol Biol Phys 89:888-98. doi:10.1016/j.ijrobp.2014.03.047

Carmona R, Zakeri K, Green G, et al. (2016) Improved method to stratify elderly patients with cancer at risk for competing events. J Clin Oncol 34:1270-77. doi:10.1200/JCO.2015.65.0739

Mell LK, Shen H, Nguyen-Tan PF, et al. (2019) Nomogram to predict the benefit of intensive treatment for locoregionally advanced head and neck cancer. Clin Cancer Res 25:7078-7088. doi:10.1158/1078-0432.CCR-19-1832

Zakeri K, Rotolo F, Lacas B, et al. (2020) Predictive classifier for intensive treatment of head and neck cancer. Cancer 126:5263-5273. doi:10.1002/cncr.33212

Mell LK, Pugh SL, Jones CU, et al. (2024) Effects of androgen deprivation therapy on prostate cancer outcomes according to competing event risk: secondary analysis of a phase 3 randomised trial. Eur Urol 85:373-381. doi:10.1016/j.eururo.2023.01.020

Alligator plot: cumulative incidence by risk group and event type

Description

Plots the cause-specific cumulative incidence functions (CIFs) of the event of interest and the competing event over time on a single panel, with risk groups distinguished by color and event types by line type (event of interest solid, competing event dashed). In higher-risk groups the event-of-interest CIF rises above the competing-event CIF and the curves separate like an alligator's open jaws.

The x-axis is truncated at the model's readout time t, with tick marks at 1-unit (e.g. 1-year) intervals. The y-axis auto-scales its cap to the smallest of 25/50/75/100% that covers the peak cumulative incidence (so low-incidence settings are not compressed), with ticks at 10% (or 5% when capped at 25%). Numbers at risk for each group are shown in an aligned panel below the plot. Gray's test across groups is attached.

Usage

gce_alligator(
  object,
  groups = 2,
  method = c("optimal", "quantile"),
  criterion = c("omegaplus", "cif"),
  xlab = "Time (years)",
  ylab = "Cumulative incidence",
  group.lab = "Risk group",
  risk.table = TRUE,
  ymax = NULL
)

Arguments

Value

A ggplot/patchwork object. The Gray's-test result, the plotted CIF data, the numbers-at-risk table, and the cutpoints are attached as attributes "grays_test", "data", "n.risk", and "cutpoints".

Examples

data(hn)
Ind <- data.frame(event = as.integer(hn$status == 1),
                  competing = as.integer(hn$status == 2))
fit <- gcecox(hn$time, Ind, hn[, c("age", "t_cat", "p16")], M = 5, t = 5)
gce_alligator(fit, groups = 2)

Calibration of predicted vs. observed omega+ by risk rank

Description

Groups subjects into equal-sized rank groups on the GCE risk score and compares the model-predicted quantity with an observed estimate in each group. For omega+ models (gcecox, gcefg) the predicted \omega^{+} (or \omega) is compared with observed cause-specific Nelson-Aalen cumulative hazards at time t: \omega^{+}_{obs} = H_1(t)/H_2(t) and \omega_{obs} = H_1(t)/(H_1(t)+H_2(t)). For a gcecif fit the predicted cumulative-incidence ratio \rho (or share \pi_1) is compared with the observed Aalen-Johansen cumulative incidence within the group: \rho_{obs} = F_1(t)/F_2(t) and \pi_{1,obs} = F_1(t)/(F_1(t)+F_2(t)). A calibration plot (predicted vs. observed with the identity line) and the underlying table are returned.

Usage

gce_calibration(object, which = NULL, groups = 5)

Arguments

Value

A ggplot object (predicted vs. observed). The calibration table is attached as attribute "data".

Examples

data(hn)
Ind <- data.frame(event = as.integer(hn$status == 1),
                  competing = as.integer(hn$status == 2))
fit <- gcecox(hn$time, Ind, hn[, c("age", "t_cat", "p16")], M = 5, t = 5)
gce_calibration(fit, which = "omegaplus", groups = 5)

Risk-score cutpoints for GCE risk stratification

Description

Splits subjects into ordered risk groups on the GCE risk score. With method = "optimal" the cutpoint(s) maximize a separation criterion by a greedy forward search; with method = "quantile" equal-sized groups are used.

The default criterion, "omegaplus", maximizes the between-group separation in log \omega^{+} (equivalently, the risk score / linear predictor). This is the additive scale on which the joint event-of-interest-vs-competing \omega^{+} contrast is defined, and – unlike separation measured on the exponential \omega^{+} scale – it is not sensitive to a few extreme high-risk subjects, so the search does not chase high-risk tails and tends to return balanced, well-separated groups. The alternative "cif" maximizes a joint Gray's test statistic (summed across both competing events) on the observed cumulative incidence.

Usage

gce_cutpoints(
  object,
  groups = 2,
  method = c("optimal", "quantile"),
  criterion = c("omegaplus", "cif"),
  min.prop = 0.1,
  nperm = 0,
  labels = NULL
)

Arguments

Value

A list with: groups (ordered factor), cutpoints (on the risk-score scale), summary (per-group n, events, mean \omega^{+} / \omega), and test. The test element reports the Lunn-McNeil \omega^{+} contrast between groups – a test that the ratio of cause-specific (or subdistribution) hazards differs across groups – with the univariable \omega^{+} ratio estimate and CI ($omegaplus), the joint test ($statistic, $df, $p.value), the per-event Gray tests for comparison ($individual_gray), and a selection-adjusted $p.perm when nperm > 0.

Examples

data(hn)
Ind <- data.frame(event = as.integer(hn$status == 1),
                  competing = as.integer(hn$status == 2))
fit <- gcecox(hn$time, Ind, hn[, c("age", "t_cat", "p16")], M = 5, t = 5)
gce_cutpoints(fit, groups = 2)

# with a selection-adjusted permutation p-value (slower)
gce_cutpoints(fit, groups = 2, nperm = 200)

Discrimination (C-index) for a GCE model

Description

Computes the cause-specific concordance (C-index) of the GCE risk score for the event of interest. For calibration (predicted vs. observed \omega^{+}), see gce_calibration.

Usage

gce_performance(object, conf.level = 0.95)

Arguments

Details

The C-index is computed with concordance using the event of interest as the outcome and the GCE risk score as the predictor. Higher risk scores denote higher risk; the statistic is oriented so that values above 0.5 indicate discrimination in the expected direction, and the "flipped" attribute of the result records whether the orientation was reversed.

Value

An object of class "gce_performance": a named numeric vector with elements estimate, se, lower, and upper, carrying attributes conf.level and flipped. Printing the object displays the C-index with its confidence interval.

See Also

gce_calibration

Examples

data(hn)
Ind <- data.frame(event = as.integer(hn$status == 1),
                  competing = as.integer(hn$status == 2))
fit <- gcecox(hn$time, Ind, hn[, c("age", "t_cat", "p16")], M = 5, t = 5)
gce_performance(fit)

Normalized risk score and relative omega+ for new subjects

Description

Computes the GCE risk score for new subjects using the covariate means/SDs and coefficients stored in a fitted model. The score is the normalized (mean-centered) linear predictor, so exp(riskscore) is the subject's omega+ relative to the average subject in the fitting data (a value of 1 = same as average, 2 = twice the average omega+). For cause-specific (gcecox) models the absolute predicted omega+ and omega at the model's time t are also returned.

Usage

gce_riskscore(object, newdata)

Arguments

Value

A data frame with riskscore and exp(riskscore) relative to the average subject: omegaplus_rel for omega+ models (with absolute omegaplus and omega at time t for cause-specific fits), or rho_rel for a gcecif fit.

See Also

gce_scaling

Examples

data(prostate)
Ind <- data.frame(as.integer(prostate$status == 1), as.integer(prostate$status == 2))
fit <- gcecox(prostate$time, Ind, prostate[, c("age", "psa", "gleason")], M = 5, t = 10)
gce_riskscore(fit, newdata = prostate[1:5, ])

Covariate scaling (means and SDs) used by the GCE risk score

Description

Returns each covariate's mean and standard deviation (from the data used to fit the model), together with the per-SD log-omega+ coefficient and omega+ ratio. These are the quantities needed to reproduce the normalized risk score for new subjects (see gce_riskscore) and to interpret how far a subject sits from the average patient.

Usage

gce_scaling(object)

Arguments

Value

A data frame with columns term, mean, sd, coef_perSD (log ratio per 1 SD) and the per-SD ratio (omegaplus_ratio_perSD for omega+ models, rho_ratio_perSD for a gcecif fit).

See Also

gce_riskscore

Examples

data(prostate)
Ind <- data.frame(as.integer(prostate$status == 1), as.integer(prostate$status == 2))
fit <- gcecox(prostate$time, Ind, prostate[, c("age", "psa", "gleason")], M = 5, t = 10)
gce_scaling(fit)

Tidy a GCE model into a data frame of estimates

Description

With effects = FALSE (the default) the omega+ ratio estimates are returned. With effects = TRUE a covariate-effects comparison table is returned that juxtaposes, for each covariate, the hazard ratio for the primary (event-of-interest) hazard, the competing-event hazard, the total (composite, any-event) hazard, and the relative hazard ratio (the omega+ ratio) from the GCE model, each with a confidence interval. This lets the investigator see whether a covariate acts mainly on the event of interest, the competing event, both, or on their balance.

Usage

gce_tidy(object, effects = FALSE, ...)

Arguments

Details

The individual- and total-event columns are estimated on the same hazard scale as the fitted model: cause-specific Cox models for a gcecox fit, and Fine-Gray subdistribution models for a gcefg fit. The total (composite) event has no competing event, so its subdistribution and cause-specific hazards coincide and a Cox model on the composite endpoint is used in both cases.

Value

If effects = FALSE, a data frame with columns term, log_ratio, ratio, conf.low, conf.high, p.value. If effects = TRUE, a data frame with columns term, primary, competing, total, and relative, each a formatted "HR (low, high)" string.

Examples

data(hn)
Ind <- data.frame(event = as.integer(hn$status == 1),
                  competing = as.integer(hn$status == 2))
fit <- gcecox(hn$time, Ind, hn[, c("age", "t_cat", "p16")], M = 5, t = 5)
gce_tidy(fit)
gce_tidy(fit, effects = TRUE)

Generalized competing event model on the cumulative-incidence-ratio scale

Description

Fits a generalized competing event (GCE) model on the cumulative incidence scale, estimating covariate effects on the ratio of cumulative incidences \rho(t) = F_1(t) / F_2(t) at a landmark time t, where F_1 and F_2 are the cumulative incidence functions of the event of interest and the competing event. \rho(t) is the odds that a subject's realized event by time t is the event of interest rather than the competing event, so gcecif() is a cumulative-incidence-scale companion to the cause-specific (gcecox) and subdistribution (gcefg) omega+ models and returns a compatible "gcemod" object.

Effects are estimated by jackknife pseudo-observation regression. Two equivalent parameterizations are provided (both estimate the effect on \log \rho(t)): method = "logit" (default) regresses pseudo-values of the bounded share \pi_1 = F_1/(F_1+F_2) by a logit-link GEE, and method = "logratio" regresses pseudo-values of \log \rho directly by an identity-link GEE (ordinary least squares). Both use a robust (sandwich) variance. A GCE risk score is formed from the model linear predictor, so exp(riskscore) is a subject's \rho relative to the average subject (see gce_riskscore, gce_scaling).

Usage

gcecif(
  Time,
  Ind,
  Cov,
  M = 5,
  t,
  method = c("logit", "logratio"),
  ridge = 0,
  conf.level = 0.95
)

Arguments

Details

Because F_k(t) = \int_0^t S(u^-)\,h_k(u)\,du with S the all-cause survivor, \rho(t) is a survival-weighted time-average of the instantaneous cause-specific hazard ratio \omega^{+}(u); under a time-constant \omega^{+} it equals \omega^{+} exactly. The ratio of cumulative incidences is a non-collapsible, marginal-over-time measure, so exp(coef) is interpreted as a conditional association with the realized event mix, not a mechanism-level per-cause effect.

The estimand is defined for a genuine competing-risks setting, in which both F_1(t) and F_2(t) are appreciable at the landmark t. If the competing event is effectively absent, or t precedes appreciable competing-event incidence, \log \rho is weakly identified; the optional ridge penalty stabilizes such fits.

Value

An object of class "gcemod" (subclass "gcecif") with: coef (log \rho ratios), result (ratio table with CIs and p-values), omnibus (Wald test that all \rho ratios equal 1), riskscore (linear predictor), scaling (covariate means/SDs and per-SD effects), F1, F2, rho0 (baseline \rho at the mean covariate profile), rho (per-subject fitted \rho), pseudo (per-subject pseudo-values), and the design matrix, time, and cause coding used by gce_tidy.

References

Andersen PK, Klein JP, Rosthoj S (2003). Generalized linear models for correlated pseudo-observations, with applications to multi-state models. Biometrika 90:15-27.

See Also

gcecox, gcefg, gce_tidy, gce_riskscore

Examples

data(prostate)
Ind <- data.frame(event = as.integer(prostate$status == 1),
                  competing = as.integer(prostate$status == 2))
Cov <- prostate[, c("age", "psa", "gleason", "t2b", "comorbidity")]
fit <- gcecif(prostate$time, Ind, Cov, t = 10)
summary(fit)
gce_tidy(fit)

Generalized competing event model via cause-specific Cox regression

Description

Fits a generalized competing event (GCE) model on the cause-specific hazard scale and estimates covariate effects on \omega^{+} (event of interest vs. competing event) using the Lunn-McNeil stacked model (see lunnmcneil), with 95% confidence intervals, Wald p-values, and an omnibus test. A GCE risk score is formed from the model linear predictor, and per-subject \omega = \omega^{+}/(1+\omega^{+}) is returned for comparison with observed values.

Usage

gcecox(
  Time,
  Ind,
  Cov,
  M = 5,
  t,
  conf.level = 0.95,
  cluster.se = TRUE,
  standardize = FALSE
)

Arguments

Value

An object of class "gcemod" (subclass "gcecox") with: coef (log \omega^{+} ratios), result (ratio table with CIs and p-values), omnibus (Lunn-McNeil global test), riskscore (linear predictor), omega / omegaplus (per-subject values at t), omegaplot and omegatimeplot (ggplot objects), and the fitted Lunn-McNeil model.

See Also

gcefg, lunnmcneil, gce_cutpoints

Examples

data(hn)
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 <- gcecox(hn$time, Ind, Cov, M = 5, t = 5)
summary(fit)

Generalized competing event model via Fine-Gray regression

Description

Fits a generalized competing event (GCE) model on the subdistribution hazard scale and estimates covariate effects on \omega^{+} using a Lunn-McNeil-style stacked Fine-Gray model (see lunnmcneil with type = "finegray"), giving 95% confidence intervals, Wald p-values, and an omnibus test from one joint model. A GCE risk score is formed from the model linear predictor, and per-subject \omega^{+} / \omega at time t are returned for comparison with observed values.

Usage

gcefg(Time, Ind, Cov, M = 5, t, conf.level = 0.95, standardize = FALSE)

Arguments

Value

An object of class "gcemod" (subclass "gcefg") with the same components as gcecox on the subdistribution scale.

See Also

gcecox, lunnmcneil

Examples

data(hn)
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 <- gcefg(hn$time, Ind, Cov, M = 5, t = 5)
summary(fit)

Head-and-neck cancer competing-risks example dataset

Description

An example head-and-neck cancer competing-risks dataset used to illustrate the gcemod functions. The data are provided for illustration only and do not represent real patients. The event of interest is cancer recurrence and the competing event is death without recurrence; recurrence is the more frequent event (recurrence-dominant, \omega^{+} > 1). Contrast with prostate.

Usage

hn

Format

A data frame with 1000 rows and 11 variables:

id

sequential subject identifier

time

follow-up time in years

status

competing-risks status: 0 = censored, 1 = cancer recurrence (event of interest), 2 = death without recurrence (competing event)

age

age at diagnosis (years)

ps

ECOG performance status (factor: 0, 1, 2)

female

sex indicator (1 = female)

smoker

smoking indicator (1 = more than 10 pack-years)

t_cat

T category (factor: 1-4)

n_cat

N category (factor: 0-3)

p16

p16 status (1 = positive)

site

primary site (factor: oropharynx, larynx, hypopharynx, oralcavity)

Lunn-McNeil stacked-model estimation of omega+ ratios

Description

Estimates covariate effects on \omega^{+} (the ratio of the hazard for the event of interest to the hazard for the competing event) from a single stacked ("augmented") model, following Lunn and McNeil (1995). Each subject contributes one record per event type; the model is stratified by event type and the covariate-by-event-type interaction estimates \log(\omega^{+}\text{ ratio}) with model-based standard errors, per-covariate Wald tests, and an omnibus multivariate Wald test.

Usage

lunnmcneil(
  Time,
  cause,
  Cov,
  event = 1,
  competing = 2,
  type = c("coxph", "finegray"),
  cluster.se = TRUE,
  conf.level = 0.95,
  standardize = FALSE,
  ...
)

Arguments

Details

Two hazard scales are supported:

type = "coxph"

Cause-specific hazards. The augmented data set duplicates each subject into an "event" stratum and a "competing" stratum; the interaction coefficient is exactly the difference in cause-specific log-hazards, i.e. \log(\omega^{+}). This is the original Lunn-McNeil construction and the interaction MLE equals the difference of two separate cause-specific Cox fits.

type = "finegray"

Subdistribution hazards. Risk-set-weighted data are generated for each event type with finegray, stacked, and fit with a single weighted, event-type-stratified coxph. The interaction estimates the difference in subdistribution log-hazards. Because subjects recur across the weighted risk sets, cluster-robust (subject) standard errors are used. This is an extension of Lunn-McNeil to the subdistribution scale.

Subject records are duplicated in both cases, so cluster-robust standard errors (clustered on subject) are the default.

Value

An object of class "gce_lunnmcneil": a list with the fitted model ($fit), the omega+ coefficient table ($omegaplus), the reference (competing-event) coefficients ($beta_competing), the omnibus test ($omnibus), the design column names ($terms), and the hazard scale ($type).

References

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

Examples

data(hn)
Cov <- hn[, c("age", "smoker", "t_cat", "n_cat", "p16")]
lunnmcneil(hn$time, hn$status, Cov, type = "coxph")
lunnmcneil(hn$time, hn$status, Cov, type = "finegray")

Prostate cancer competing-risks example dataset

Description

An example prostate cancer competing-risks dataset used to illustrate the gcemod functions. The data are provided for illustration only and do not represent real patients. The event of interest is distant metastasis or prostate-cancer death and the competing event is death from other causes; the competing event is the more frequent (competing-death-dominant, \omega^{+} < 1). Contrast with hn.

Usage

prostate

Format

A data frame with 1000 rows and 9 variables:

id

sequential subject identifier

time

follow-up time in years

status

competing-risks status: 0 = censored, 1 = distant metastasis or prostate-cancer death (event of interest), 2 = death from other causes (competing event)

age

age at diagnosis (years)

ps

performance status (factor: 0, 1)

psa

pre-treatment PSA (ng/mL, capped at 20)

gleason

Gleason grade group (factor: "<=6", "3+4", "4+3", ">=8")

t2b

stage indicator (1 = T2b or higher)

comorbidity

significant comorbidity indicator (1 = yes)