Package 'gcemod'

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

object	a "gcemod" object from gcecox, gcefg, or gcecif.
groups	a gce_cutpoints result, a grouping factor of length n, or a single integer giving the number of risk groups (default 2). When an integer is given, the groups are derived with gce_cutpoints using method.
method	grouping rule used when groups is an integer: "optimal" (default) or "quantile" (equal-sized groups). Ignored when groups is a cutpoints object or a factor.
criterion	optimal-search criterion when groups is an integer and method = "optimal": "omegaplus" (default, log-$\omega^{+}$ separation; not tail-sensitive) or "cif" (joint Gray statistic on the observed cumulative incidence of both events). See gce_cutpoints.
xlab, ylab	axis labels.
group.lab	legend title for the risk groups (default "Risk group").
risk.table	logical; show numbers at risk in an aligned panel below the plot (default TRUE).
ymax	optional numeric upper cap for the y-axis (0-1). If NULL (default) the cap auto-scales to 25/50/75/100% based on the peak incidence.

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

object	a "gcemod" object from gcecox, gcefg, or gcecif.
which	quantity to calibrate. For omega+ models "omegaplus" (default) or "omega"; for a gcecif fit "rho" (default) or "share". If NULL the model-appropriate default is used.
groups	number of rank groups (default 5).

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

object	a "gcemod" object from gcecox, gcefg, or gcecif. For a gcecif fit the per-group summary reports mean $\rho$ and share instead of $\omega^{+}$ / $\omega$; the "omegaplus" criterion still means separation on the risk-score (linear-predictor) scale.
groups	number of groups (default 2).
method	"optimal" (default) or "quantile".
criterion	optimal-search criterion: "omegaplus" (default, log-$\omega^{+}$ separation) or "cif" (joint Gray statistic on observed cumulative incidence of both events).
min.prop	minimum proportion of subjects per group (default 0.1).
nperm	number of permutations for the selection-adjusted p-value (default 0 = not computed). When nperm > 0, the maximally-selected joint Gray statistic is permuted (risk score vs. outcome) to give a p-value accounting for the data-driven choice of cutpoint.
labels	optional character labels for the groups (low to high risk).

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

object	a "gcemod" object.
conf.level	confidence level for the C-index interval (default 0.95).

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

object	a "gcemod" object from gcecox, gcefg, or gcecif.
newdata	a data frame containing the covariates used to fit object (same column names). Factor levels not seen in the fitting data are treated as the reference level.

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

object

a "gcemod" object from gcecox, gcefg, or gcecif.

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

object	a "gcemod" object from gcecox, gcefg, or gcecif. For a gcecif fit the relative column is the cumulative-incidence ratio (rho ratio) and the primary/competing/total columns are cause-specific hazard ratios.
effects	logical; if TRUE, return the primary/competing/total/ relative effects comparison table instead of the omega+ ratio table. Requires an object fit with gcemod >= 0.3.0 (which stores the design matrix). Default FALSE.
...	unused.

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

Time	numeric vector of follow-up times.
Ind	a data frame of event indicators: column 1 = event of interest, column 2 = competing event (0/1).
Cov	a data frame of covariates (factors are expanded automatically).
M	number of bins for the $\rho$ risk plot (default 5).
t	landmark time at which $\rho$ is evaluated.
method	estimation scale: "logit" (default; logit-link GEE on the share) or "logratio" (identity-link GEE on $\log\rho$).
ridge	non-negative ridge penalty on the covariate coefficients (not the intercept), used only with method = "logit" to stabilize weakly identified fits. Default 0.
conf.level	confidence level (default 0.95).

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

Time	numeric vector of follow-up times.
Ind	a data frame of event indicators. Column 1 = event of interest, column 2 = competing event (0/1). A third column, if present, is ignored.
Cov	a data frame of covariates.
M	number of bins for the $\omega^{+}$ risk plot (default 5).
t	time point at which per-subject $\omega$ / $\omega^{+}$ are evaluated.
conf.level	confidence level (default 0.95).
cluster.se	logical; cluster-robust SEs in the Lunn-McNeil model (default TRUE).
standardize	logical; controls only the reported omega+ ratio scale (default FALSE = per natural covariate unit; TRUE = per 1 SD). The risk score is always the normalized (mean-centered) linear predictor regardless of this setting, so cutpoints, plots, and the predicted-vs-observed omega calibration are unaffected. Covariate means/SDs are available in $scaling (see gce_scaling).

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

Time	numeric vector of follow-up times.
Ind	a data frame of event indicators. Column 1 = event of interest, column 2 = competing event (0/1).
Cov	a data frame of covariates (factors are expanded to indicators).
M	number of bins for the $\omega^{+}$ risk plot (default 5).
t	time point at which per-subject $\omega^{+}$ / $\omega$ are evaluated.
conf.level	confidence level (default 0.95).
standardize	logical; controls only the reported omega+ ratio scale (default FALSE = per natural covariate unit; TRUE = per 1 SD). The risk score is always the normalized (mean-centered) linear predictor. Covariate means/SDs are in $scaling (see gce_scaling).

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

Time	numeric vector of follow-up times.
cause	event indicator per subject: 0 = censored, and the codes in event and competing for the two event types.
Cov	a data frame of covariates (factors are expanded to indicators).
event	code marking the event of interest (default 1).
competing	code marking the competing event (default 2).
type	"coxph" (cause-specific, default) or "finegray" (subdistribution).
cluster.se	logical; cluster-robust (subject) SEs (default TRUE).
conf.level	confidence level (default 0.95).
standardize	logical; controls only the reported coefficient scale (default FALSE). FALSE reports omega+ ratios per natural covariate unit; TRUE reports them per 1 SD. This is a reparameterization, so Wald tests, p-values, and the omnibus test are identical either way. The covariate means/SDs ($center, $scale) are always returned.
...	passed to coxph.

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)

---