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.

gcemod-methods

Summaries, extractors, prediction, and plots for fitted models.

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

Interface

Models are specified with a formula, e.g. gcecox(Surv(time, status) ~ age + stage, data = d), and return classed objects with the usual methods (print, summary, coef, vcov, confint, predict, plot, ...); see gcemod-methods. Surv is re-exported so that the formula can be written without attaching survival.

Author(s)

Maintainer: Loren Mell lmell@ucsd.edu

Authors:

• Loren Mell lmell@ucsd.edu

Other contributors:

• Morten Overgaard moov@ph.au.dk (proposed the stacked pseudo-observation estimator used by gcecif()) [contributor]

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 xmax (by default 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,
  xmax = NULL
)

## S3 method for class 'gce_alligator'
print(x, ...)

## S3 method for class 'gce_alligator'
plot(x, ...)

Arguments

Value

An object of class "gce_alligator" with print and plot methods. It is a list holding the cumulative incidence curves ($data), the numbers at risk ($n.risk), Gray's tests across groups ($grays_test), and the cutpoints. Printing shows the group sizes and Gray's tests; plot draws the figure and returns the ggplot/patchwork object invisibly. The display arguments (xlab, ylab, group.lab, risk.table, ymax) are stored and used by plot.

Examples

data(hn)
fit <- gcecox(Surv(time, status) ~ age + t_cat + p16, data = hn, t = 5)
al <- gce_alligator(fit, groups = 2)
al
plot(al)

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, t = NULL)

## S3 method for class 'gce_calibration'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'gce_calibration'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

## S3 method for class 'gce_calibration'
plot(x, ...)

Arguments

Value

An object of class "gce_calibration" with print, plot, and as.data.frame methods. Printing shows the calibration table (group, n, predicted, observed); plot draws predicted against observed with the identity line and returns the ggplot object invisibly.

Examples

data(hn)
fit <- gcecox(Surv(time, status) ~ age + t_cat + p16, data = hn, t = 5)
cal <- gce_calibration(fit, which = "omegaplus", groups = 5)
cal
plot(cal)

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
)

## S3 method for class 'gce_cutpoints'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'gce_cutpoints'
summary(object, ...)

## S3 method for class 'summary.gce_cutpoints'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'gce_cutpoints'
plot(x, ...)

Arguments

Value

An object of class "gce_cutpoints", with print, summary, and plot methods. print shows the cutpoints, the per-group table, and the group contrast; summary returns the per-group table (a data frame) together with the test results as an object of class "summary.gce_cutpoints"; plot draws the distribution of the risk score with the cutpoints marked. The object is a list with: groups (ordered factor), cutpoints (on the risk-score scale), score (the risk score), 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)
fit <- gcecox(Surv(time, status) ~ age + t_cat + p16, data = hn, t = 5)
cp <- gce_cutpoints(fit, groups = 2)
cp
summary(cp)$groups
plot(cp)
table(cp$groups)

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

GCE contrasts over time: a visual check of relative-hazard constancy

Description

Plots a nonparametric estimate of a chosen GCE contrast over follow-up, as a visual assessment of the proportional relative-hazards (PRH) assumption and of the time range over which it holds approximately. All estimates are built from the Nelson-Aalen cause-specific cumulative hazards \Lambda_k(t) and the Aalen-Johansen cumulative incidences F_k(t):

"omega" (default)

\Lambda_1(t)/\{\Lambda_1(t)+\Lambda_2(t)\}, the bounded cumulative hazard share, plotted on [0,1].

"omegaplus"

\omega^{+}(t)=\Lambda_1(t)/\Lambda_2(t), the cause-specific relative hazard, on a log scale.

"omegatilde"

\tilde\Lambda_1(t)/\{\tilde\Lambda_1(t)+\tilde\Lambda_2(t)\}, the subdistribution analogue of "omega", with \tilde\Lambda_k(t) = -\log\{1-F_k(t)\}, on [0,1].

"omegatildeplus"

\tilde\omega^{+}(t)=\tilde\Lambda_1(t)/\tilde\Lambda_2(t), the subdistribution relative hazard, on a log scale.

"rho"

\rho(t)=F_1(t)/F_2(t), the cumulative-incidence ratio, on a log scale.

"pi"

\pi(t)=F_1(t)/\{F_1(t)+F_2(t)\} = \rho(t)/\{1+\rho(t)\}, the bounded share form of \rho (so \mathrm{logit}\,\pi = \log\rho), i.e. the probability that a realized event by t is the event of interest, on [0,1].

If the cause-specific relative hazard is constant in time then \Lambda_1(t)=\omega^{+}\Lambda_2(t) and, because F_k(t)=\int_0^t S(u^-)\lambda_k(u)\,du, also F_1(t)=\omega^{+}F_2(t): the "omega", "omegaplus", "rho", and "pi" displays are then exactly flat, and drift in any of them delineates where a single time-constant summary is and is not adequate. The subdistribution displays are flat only approximately (exactly so when both incidences are small, since \tilde\Lambda_k \approx F_k).

With groups, one curve is drawn per risk group. Approximately parallel group curves on a log-scaled contrast indicate a time-constant group contrast—the modeling assumption of the stacked estimators—even where the individual curves drift.

Usage

gce_omegatime(
  object,
  which = c("omega", "omegaplus", "omegatilde", "omegatildeplus", "rho", "pi"),
  groups = NULL,
  method = "optimal",
  xmax = NULL
)

## S3 method for class 'gce_omegatime'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'gce_omegatime'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

## S3 method for class 'gce_omegatime'
plot(x, ...)

Arguments

Value

An object of class "gce_omegatime" with print, plot, and as.data.frame methods. The step-function estimates (group, time, value) are in $data. Printing gives, per group, the range of the contrast and its value at the readout time t; plot draws the curves (the readout time is marked with a dashed vertical line) and returns the ggplot object invisibly.

See Also

gce_cutpoints, gce_alligator

Examples

data(hn)
fit <- gcecox(Surv(time, status) ~ age + t_cat + p16, data = hn, t = 5)
ot <- gce_omegatime(fit)
ot
plot(ot)
plot(gce_omegatime(fit, which = "omegaplus", groups = 2, xmax = 10))

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)
fit <- gcecox(Surv(time, status) ~ age + t_cat + p16, data = hn, 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. This is a convenience wrapper around predict(), which offers the same quantities individually (type = "lp", "risk", "ratio", "share") and at arbitrary times. 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)
fit <- gcecox(Surv(time, status) ~ age + psa + gleason, data = prostate, 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)
fit <- gcecox(Surv(time, status) ~ age + psa + gleason, data = prostate, 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)
fit <- gcecox(Surv(time, status) ~ age + t_cat + p16, data = hn, 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 a stacked pseudo-observation regression: each subject contributes one jackknife pseudo-observation for F_1(t) and one for F_2(t) (both mean-type functionals of the Aalen-Johansen estimator, for which conditional unbiasedness of pseudo-observations is established), and a single log-link GEE is fit to the stacked pairs with an event-type indicator. The main effects model \log F_2(t \mid z); the event-type interaction coefficients estimate \log \rho(t \mid z) = \log F_1(t \mid z) - \log F_2(t \mid z) directly. Because each subject appears twice, a cluster-robust (subject) sandwich variance is used, which captures the dependence between a subject's two pseudo-observations induced by the shared all-cause survivor. This is the cumulative-incidence-scale analogue of the stacked Lunn-McNeil construction used on the two hazard scales. A GCE risk score is formed from the interaction linear predictor, so exp(riskscore) is a subject's \rho relative to the average subject (see gce_riskscore, gce_scaling).

Usage

gcecif(formula, ...)

## S3 method for class 'formula'
gcecif(
  formula,
  data,
  subset,
  na.action,
  t,
  event = NULL,
  competing = NULL,
  conf.level = 0.95,
  contrasts = NULL,
  model = TRUE,
  ...
)

## Default S3 method:
gcecif(formula, Ind, Cov, M = NULL, t, conf.level = 0.95, ...)

Arguments

Details

The model is specified with a formula whose response is Surv(time, status), as for gcecox; the landmark time t is part of the estimand and must be supplied. The default method is the matrix interface of earlier versions, retained for backward compatibility.

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 estimating function is linear in the two per-cause (mean-type) pseudo-observations, so its conditional unbiasedness follows from the standard pseudo-observation theory (Graw et al. 2009; Overgaard et al. 2017); pseudo-observations of the nonlinear ratio itself are not used. The working model comprises a log-linear model for F_2(t \mid z) (the nuisance main effects, reported as $coef_f2) and a log-linear model for \rho(t \mid z) (the target interactions, reported as $coef).

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 and the fit may fail or return very large standard errors; this should be read as a signal that the competing-risks premise does not hold.

Value

An object of class c("gcecif", "gcemod"); see gcemod-methods for the available methods. The list contains coef and var (log \rho ratios and their covariance matrix), omnibus (Wald test that all \rho ratios equal 1), riskscore (interaction linear predictor, mean-centered), scaling (covariate means/SDs and per-SD effects), F1, F2, rho0 (marginal F_1(t)/F_2(t)), rho (per-subject fitted \rho), pseudo (a two-column data frame of per-cause pseudo-observations), coef_rho and coef_f2 (the full log-linear models for \rho(t \mid z) and the nuisance F_2(t \mid z), intercepts included), converged, and the model terms, xlevels, contrasts, and model frame.

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.

Graw F, Gerds TA, Schumacher M (2009). On pseudo-values for regression analysis in competing risks models. Lifetime Data Analysis 15:241-255.

Overgaard M, Parner ET, Pedersen J (2017). Asymptotic theory of generalized estimating equations based on jack-knife pseudo-observations. The Annals of Statistics 45:1988-2015.

Overgaard M, Parner ET, Pedersen J (2018). Estimating the variance in a pseudo-observation scheme with competing risks. Scandinavian Journal of Statistics 45:923-940.

See Also

gcecox, gcefg, gcemod-methods, gce_tidy

Examples

data(prostate)
fit <- gcecif(Surv(time, status) ~ age + psa + gleason + t2b + comorbidity,
              data = prostate, t = 10)
summary(fit)
predict(fit, newdata = prostate[1:3, ], type = "ratio")

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 confidence intervals, Wald tests, and an omnibus test. A GCE risk score is formed from the model linear predictor.

Usage

gcecox(formula, ...)

## S3 method for class 'formula'
gcecox(
  formula,
  data,
  subset,
  na.action,
  t = NULL,
  event = NULL,
  competing = NULL,
  conf.level = 0.95,
  cluster.se = TRUE,
  standardize = FALSE,
  contrasts = NULL,
  model = TRUE,
  ...
)

## Default S3 method:
gcecox(
  formula,
  Ind,
  Cov,
  M = NULL,
  t = NULL,
  conf.level = 0.95,
  cluster.se = TRUE,
  standardize = FALSE,
  ...
)

Arguments

Details

The model is specified with a formula whose response is Surv(time, status); see the Interfaces section of lunnmcneil for the accepted status codings. The fitted object stores the model terms, factor levels, and contrasts, so that predict() can evaluate the risk score and the predicted \omega^{+} for new data.

Fitting does not depend on a time point. The argument t only sets the default readout time for quantities that are functions of time: the predicted absolute \omega^{+}(t) and \omega(t) (predict(type = "ratio"), type = "share"), calibration (gce_calibration), and the time axis of the displays. It can be overridden in each of those functions.

The default method is the matrix interface of earlier versions, gcecox(Time, Ind, Cov, t = ), retained for backward compatibility.

Value

An object of class c("gcecox", "gcemod"). The components most often needed are available through methods rather than by name; see gcemod-methods (print, summary, coef, vcov, confint, logLik, nobs, formula, terms, model.frame, model.matrix, predict, fitted, update, plot). The list also contains coef and var (log omega+ ratios and their covariance matrix on the reported scale), riskscore (mean-centered linear predictor), omegaplus and omega (per-subject predictions at t), omnibus, lunnmcneil (the stacked-model fit), scaling, and the model terms, xlevels, contrasts, and model frame.

See Also

gcefg, gcecif, lunnmcneil, gcemod-methods, gce_cutpoints

Examples

data(hn)
fit <- gcecox(Surv(time, status) ~ age + smoker + t_cat + n_cat + p16,
              data = hn)
fit
summary(fit)
confint(fit)

# risk score and predicted omega+ at 5 years for new subjects
predict(fit, newdata = hn[1:3, ], type = "lp")
predict(fit, newdata = hn[1:3, ], type = "ratio", times = 5)

# matrix interface
Ind <- data.frame(event = as.integer(hn$status == 1),
                  competing = as.integer(hn$status == 2))
gcecox(hn$time, Ind, hn[, c("age", "smoker", "p16")], t = 5)

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 \tilde\omega^{+} using a Lunn-McNeil-style stacked Fine-Gray model (see lunnmcneil with type = "finegray"), giving confidence intervals, Wald tests, and an omnibus test from one joint model. A GCE risk score is formed from the model linear predictor.

Usage

gcefg(formula, ...)

## S3 method for class 'formula'
gcefg(
  formula,
  data,
  subset,
  na.action,
  t = NULL,
  event = NULL,
  competing = NULL,
  conf.level = 0.95,
  standardize = FALSE,
  se = c("robust", "boot"),
  B = 200,
  contrasts = NULL,
  model = TRUE,
  ...
)

## Default S3 method:
gcefg(
  formula,
  Ind,
  Cov,
  M = NULL,
  t = NULL,
  conf.level = 0.95,
  standardize = FALSE,
  se = c("robust", "boot"),
  B = 200,
  ...
)

Arguments

Details

The interface, the role of t, and the returned object are as for gcecox. Predicted absolute \tilde\omega^{+}(t) values are obtained from two separate Fine-Gray fits (crr), stored as $fit1 and $fit2.

Value

An object of class c("gcefg", "gcemod") with the same components and methods as a gcecox fit, on the subdistribution scale; see gcemod-methods.

See Also

gcecox, lunnmcneil, gcemod-methods

Examples

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

Methods for fitted generalized competing event models

Description

Methods for objects of class "gcemod" returned by gcecox, gcefg, and gcecif, and for the "gce_lunnmcneil" objects returned by lunnmcneil.

Usage

## S3 method for class 'gcemod'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'gce_lunnmcneil'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'gcemod'
summary(object, conf.level = object$conf.level, ...)

## S3 method for class 'gce_lunnmcneil'
summary(object, conf.level = object$conf.level, ...)

## S3 method for class 'summary.gcemod'
print(
  x,
  digits = max(3L, getOption("digits") - 3L),
  signif.stars = getOption("show.signif.stars"),
  ...
)

## S3 method for class 'gcemod'
coef(object, ...)

## S3 method for class 'gce_lunnmcneil'
coef(object, ...)

## S3 method for class 'gcemod'
vcov(object, ...)

## S3 method for class 'gce_lunnmcneil'
vcov(object, ...)

## S3 method for class 'gcemod'
confint(object, parm, level = 0.95, exponentiate = FALSE, ...)

## S3 method for class 'gce_lunnmcneil'
confint(object, parm, level = 0.95, exponentiate = FALSE, ...)

## S3 method for class 'gcemod'
nobs(object, ...)

## S3 method for class 'gce_lunnmcneil'
nobs(object, ...)

## S3 method for class 'gcemod'
formula(x, ...)

## S3 method for class 'gce_lunnmcneil'
formula(x, ...)

## S3 method for class 'gcemod'
terms(x, ...)

## S3 method for class 'gce_lunnmcneil'
terms(x, ...)

## S3 method for class 'gcemod'
model.frame(formula, ...)

## S3 method for class 'gcemod'
model.matrix(object, ...)

## S3 method for class 'gcemod'
logLik(object, ...)

## S3 method for class 'gce_lunnmcneil'
logLik(object, ...)

## S3 method for class 'gcemod'
predict(
  object,
  newdata = NULL,
  type = c("lp", "risk", "ratio", "share"),
  times = NULL,
  ...
)

## S3 method for class 'gcemod'
fitted(object, type = c("lp", "risk", "ratio", "share"), ...)

## S3 method for class 'gcemod'
plot(
  x,
  which = c("riskscore", "omegatime", "alligator", "calibration"),
  bins = 5,
  ...
)

Arguments

Details

print

a short display: the call, the sample and event counts, and the estimated ratios exp(coef).

summary

returns an object of class "summary.gcemod" (or "summary.gce_lunnmcneil") holding the full coefficient table (estimate, ratio, standard error, Wald z, p-value), the confidence intervals for the ratios, the omnibus Wald test, and the baseline quantity of the model (\omega^{+}(t) of the average subject, or F_1(t), F_2(t), and \rho(t)). Its print method formats the table with printCoefmat.

coef, vcov, confint

the log ratios (log \omega^{+} ratios; log \rho ratios for gcecif), their covariance matrix, and Wald confidence intervals, all on the reported coefficient scale (per unit, or per SD if the model was fitted with standardize = TRUE). confint(exponentiate = TRUE) returns intervals for the ratios.

logLik

the maximized log partial likelihood of the stacked model for gcecox fits, and the corresponding weighted log partial pseudo-likelihood for gcefg fits (useful for comparing nested specifications descriptively; it is not a genuine likelihood). gcecif is estimated by generalized estimating equations and has no likelihood, so logLik signals an error.

nobs, formula, terms, model.frame, model.matrix

the number of subjects and the stored model specification. update() works through the stored call.

predict, fitted

see below.

plot

dispatches to the display functions; see which.

Value

summary returns an object of class "summary.gcemod" or "summary.gce_lunnmcneil"; coef a named numeric vector; vcov a matrix; confint a two-column matrix; logLik an object of class "logLik"; predict and fitted a numeric vector (or matrix, for several times); plot a ggplot object, invisibly. print methods return their argument invisibly.

Prediction

predict() rebuilds the design matrix for newdata from the stored terms, factor levels, and contrasts. type selects:

"lp"

the GCE risk score: the linear predictor centered at the covariate means of the fitting data.

"risk"

exp(lp), the subject's ratio relative to the average subject.

"ratio"

the absolute predicted ratio: \omega^{+}(t) for gcecox (the marginal Nelson-Aalen ratio H_1(t)/H_2(t) times exp(lp)), \tilde\omega^{+}(t) for gcefg, and \rho(t) for gcecif.

"share"

ratio / (1 + ratio): \omega(t), or the share \pi_1(t) for gcecif.

times defaults to the readout time t of the fit. With several times a matrix with one column per time is returned. For gcecif the landmark time is part of the fitted model and cannot be changed. fitted() is predict() on the fitting data.

Examples

data(hn)
fit <- gcecox(Surv(time, status) ~ age + smoker + p16, data = hn, t = 5)

fit                       # short display
s <- summary(fit)         # full table, returned as an object
s$coefficients
coef(fit); vcov(fit)
confint(fit, exponentiate = TRUE)
logLik(fit); nobs(fit); formula(fit)

head(fitted(fit))                                  # risk score
predict(fit, hn[1:3, ], type = "ratio", times = c(2, 5))

fit2 <- update(fit, . ~ . - smoker)
plot(fit, which = "riskscore", bins = 5)

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.

lunnmcneil() is the estimation engine behind gcecox and gcefg; use it directly when only the coefficient table and tests are needed, and the GCE fitters when the risk score, predictions, and displays are wanted as well.

Usage

lunnmcneil(formula, ...)

## S3 method for class 'formula'
lunnmcneil(
  formula,
  data,
  subset,
  na.action,
  event = NULL,
  competing = NULL,
  type = c("coxph", "finegray"),
  cluster.se = TRUE,
  conf.level = 0.95,
  standardize = FALSE,
  se = c("robust", "boot"),
  B = 200,
  contrasts = NULL,
  ...
)

## Default S3 method:
lunnmcneil(
  formula,
  cause,
  Cov,
  event = 1,
  competing = 2,
  type = c("coxph", "finegray"),
  cluster.se = TRUE,
  conf.level = 0.95,
  standardize = FALSE,
  se = c("robust", "boot"),
  B = 200,
  ...
)

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 stacked model ($fit), the omega+ coefficient table ($omegaplus), the coefficient vector and its covariance matrix on the reported scale ($coef, $var), the reference (competing-event) coefficients ($beta_competing), the omnibus test ($omnibus), the model terms, xlevels, contrasts and formula, and the hazard scale ($type). Methods are available for print, summary, coef, vcov, confint, nobs, logLik, formula, and terms.

Interfaces

The formula method is the recommended interface. The response is written Surv(time, status), where status is either numeric with 0 = censored and the codes given by event and competing (defaults 1 and 2), or a factor whose first level denotes censoring (the convention of Surv for multi-state data), in which case event and competing default to the second and third levels. The right-hand side accepts the usual model-formula syntax (factors, interactions, transformations, splines). The model terms, factor levels, and contrasts are stored in the fitted object.

The default method is the matrix interface of earlier versions (lunnmcneil(Time, cause, Cov)); it is retained for backward compatibility and returns the same object.

References

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

See Also

gcecox, gcefg, gcemod-methods

Examples

data(hn)
fit <- lunnmcneil(Surv(time, status) ~ age + smoker + t_cat + n_cat + p16,
                  data = hn)
fit
summary(fit)
confint(fit)

# subdistribution scale
lunnmcneil(Surv(time, status) ~ age + smoker + p16, data = hn,
           type = "finegray")

# matrix interface
lunnmcneil(hn$time, hn$status, hn[, c("age", "smoker", "p16")])

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)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

survival

Surv()