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 colour 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 or gcefg.
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 $\omega^{+}$ (or $\omega$) with an observed estimate in each group. Observed values use cause-specific Nelson-Aalen cumulative hazards within the group at time t: $\omega^{+}_{obs} = H_1(t)/H_2(t)$ and $\omega_{obs} = H_1(t)/(H_1(t)+H_2(t))$. A calibration plot (predicted vs. observed with the identity line) and the underlying table are returned.

Usage

gce_calibration(object, which = c("omegaplus", "omega"), groups = 5)

Arguments

object	a "gcemod" object from gcecox or gcefg.
which	"omegaplus" (default) or "omega".
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 or gcefg.
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, nperm = 200)

---

Discrimination (C-index) for a GCE model

Description

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

Value

Invisibly, a named numeric vector cindex (estimate, se, lower, upper). The value is also printed.

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.
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, omegaplus_rel (= exp(riskscore)), and, for cause-specific models, omegaplus and omega (absolute, at time t).

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

Value

A data frame with columns term, mean, sd, coef_perSD (log-omega+ ratio per 1 SD) and omegaplus_ratio_perSD.

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 or gcefg.
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 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)

---