This release redesigns the user interface. Estimates are unchanged: every estimate, standard error, test, risk score, cutpoint, and calibration value agrees with 0.3.1 on the bundled data.
gcecox(), gcefg(),
gcecif(), and lunnmcneil() are now generics
with a formula method, which is the documented
interface:
fit <- gcecox(Surv(time, status) ~ age + smoker + t_cat + p16, data = hn)The response is Surv(time, status), with
status either numeric (0 = censored, 1 = event of interest,
2 = competing event; other codes via event /
competing) or a factor whose first level is censoring. The
right-hand side takes the usual formula syntax (factors, interactions,
transformations, splines), and data, subset,
na.action, and contrasts are supported.
Surv() is re-exported.
Fitted objects store the model terms, factor levels,
and contrasts, so the design matrix is rebuilt correctly for new
data.
The matrix interface of earlier versions
(gcecox(Time, Ind, Cov, t = ),
lunnmcneil(Time, cause, Cov)) is retained as the default
method and gives the same results as before.
print() is now a short display (call, counts,
estimated ratios). summary() returns an object of class
"summary.gcemod" (or "summary.gce_lunnmcneil")
with the full coefficient table, confidence intervals, omnibus test, and
baseline quantity; its print method uses printCoefmat().
Previously the two printed the same table and summary()
returned an unclassed list.
New methods for "gcemod" objects:
coef(), vcov(), confint(),
nobs(), formula(), terms(),
model.frame(), model.matrix(),
logLik(), predict(), fitted(),
and plot(); update() works through the stored
call. "gce_lunnmcneil" objects gain summary(),
coef(), vcov(), confint(),
nobs(), logLik(), formula(), and
terms().
predict() returns the risk score
(type = "lp"), the ratio relative to the average subject
("risk"), or the absolute predicted ratio and share
("ratio", "share") at one or several
times, for the fitting data or newdata.
gce_riskscore() is now a wrapper around it.
gce_cutpoints() returns an object of class
"gce_cutpoints" with print(),
summary(), and plot() methods (it was an
unclassed list; its elements are unchanged).
gce_calibration(), gce_omegatime(), and
gce_alligator() now return classed objects holding the
computed results, with print() and plot()
methods (and as.data.frame() for the first two).
Behaviour change: calling these functions at the prompt
prints a numerical summary; wrap the call in plot() to draw
the figure, or use
plot(fit, which = "calibration" | "omegatime" | "alligator").
The plotted data, previously attached as attributes of the ggplot
object, are now list elements ($table, $data,
$n.risk, $grays_test).
gce_calibration() gains t, and
gce_alligator() gains xmax, to set the time
independently of the fit.
M is deprecated. It only set the
number of bins of a plot stored inside the fitted object and had no
effect on estimation. It is ignored with a warning in the matrix
interface; the binned ratio-versus-risk-score plot is now
plot(fit, which = "riskscore", bins = ). The stored ggplot
objects $omegaplot, $omegatimeplot, and
$ratioplot are removed.
t is optional for gcecox() and
gcefg(). Estimation on the two hazard scales does
not involve a time point; t only sets the default readout
time for time-dependent predictions and displays, and defaults to the
largest event time. It remains required for gcecif(), where
it defines the estimand.
Coefficient tables are stored at full precision
($result was rounded to five decimals); rounding now
happens only when printing.
gce_omegatime() no longer emits “Ignoring unknown
labels: colour” when drawn without groups.
predict() and gce_riskscore() signal an
error for factor levels not present in the fitting data (they were
silently treated as the reference level).
In "gce_lunnmcneil" objects, the design column names
are now $covnames; $terms holds the model
terms object.
The coefficient table stored in a gcefg() fit
(fit$result) now heads its first column
exp(coef) (omega-tilde+ ratio); it was previously labelled
as the cause-specific omega+ ratio. The values are
unchanged.
gcecif() estimator replaced. The
cumulative-incidence-ratio model is estimated by a stacked
pseudo-observation regression, wherein each subject contributes one
jackknife pseudo-observation for each of F1(t) and
F2(t) — mean-type functionals of the Aalen-Johansen
estimator, for which conditional unbiasedness of the 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 F2(t | z); the event-type interactions estimate
log rho(t | z) directly; the variance is cluster-robust on
subject, capturing 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 by gcecox() and gcefg().
The previous estimator regressed jackknife pseudo-values of the
nonlinear ratio itself, so estimates could be biased for large covariate
effects. The method and ridge arguments are
removed; $pseudo is now a two-column data frame of
per-cause pseudo-observations, and $coef_f2 and
$converged are added. The interface, risk score, and
downstream tooling (gce_cutpoints(),
gce_alligator(), gce_calibration(),
gce_riskscore(), gce_tidy()) are
unchanged.
New gce_omegatime() PRH diagnostic
plot. Plots a nonparametric estimate of a chosen GCE contrast
over follow-up, for the whole sample or within risk groups
(groups, resolved as in gce_cutpoints()). The
which argument selects the contrast: "omega"
(default; the Nelson-Aalen cumulative hazard share
Lambda1/(Lambda1+Lambda2)), "omegaplus" (Lambda1/Lambda2,
log scale), "omegatilde" / "omegatildeplus"
(the subdistribution analogues built from -log(1 - Fk) with
Aalen-Johansen CIFs), "rho" (F1/F2, log scale), or
"pi" (F1/(F1+F2), the bounded share form of rho, so
logit(pi) = log(rho)). Under a time-constant cause-specific relative
hazard the omega, omega+, rho, and pi displays are exactly flat, so a
flat curve is a direct visual check of the proportional relative-hazards
assumption — and drift shows the time range over which it is
approximately valid; with groups, parallel curves on a log scale
indicate a time-constant group contrast. xmax truncates the
display to avoid over-reading the sparse late-follow-up tail. Successor
to the omega-versus-time display in gcerisk.
Bootstrap variance option for the subdistribution
scale. gcefg() and
lunnmcneil(type = "finegray") gain
se = c("robust", "boot") and B. The default
cluster-robust variance treats the IPCW censoring weights as fixed,
which errs on the conservative side; se = "boot" performs a
subject bootstrap that re-estimates the weights within each resample, so
the variance fully accounts for their estimation.
Morten Overgaard, who proposed the corrected gcecif() estimator, is credited as a package contributor.
First CRAN release.
gce_tidy(effects = TRUE) covariate-effects
table. gce_tidy() has an effects
argument. With effects = TRUE it returns a comparison table
that, for each covariate, juxtaposes the hazard ratio for the primary
event, the competing event, the total (composite) event, and the
relative hazard ratio (the omega+ ratio) from the GCE model, each with a
confidence interval. The individual/total columns are estimated on the
model’s own scale — cause- specific Cox for gcecox(),
Fine-Gray subdistribution for gcefg() — so the reader can
see whether a covariate acts on the event of interest, the competing
event, both, or their balance. gcecox() and
gcefg() retain the model design matrix on the fitted object
($design) to support this.gcecif() cumulative-incidence-ratio GCE
model. A third GCE estimator alongside gcecox()
and gcefg(), on the cumulative-incidence scale: covariate
effects on the ratio of cumulative incidences rho(t) = F1(t)/F2(t) at a
landmark time t — the odds that a subject’s realized event
by time t is the event of interest, a cumulative-incidence-scale
analogue of omega+. Uses jackknife pseudo-observation regression of the
Aalen-Johansen estimator with a robust variance;
method = "logit" (default) fits the bounded share via a
logit-link GEE, method = "logratio" fits log(F1/F2) via an
identity-link GEE. Returns a gcemod object, so it works
with gce_tidy(), gce_scaling(), and
gce_riskscore() like the other fits. The full
methodological development is described in a companion manuscript (in
preparation).gcecox(),
gcefg(), and lunnmcneil() report omega+ ratios
per natural covariate unit by default; standardize = TRUE
reports them per 1 SD instead (a reparameterization, so p-values and the
omnibus test are identical). The GCE risk score is always the
normalized, mean-centered linear predictor, so the average patient has
risk score 0 and cutpoints / plots / calibration are unaffected by the
reporting choice.gce_scaling() and
gce_riskscore(). gce_scaling()
returns each covariate’s mean, SD, and per-SD omega+ ratio.
gce_riskscore(fit, newdata) scores new subjects with the
stored means/SDs and coefficients, returning the risk score and
omegaplus_rel = the subject’s omega+ relative to the
average subject in the fitting data (with absolute omega+/omega for
cause-specific models). $scaling, $center,
$scale, $coef_raw, and $coef_std
are stored on the fitted object.criterion = "cif" (joint Gray
statistic across both events) is available as an alternative.gce_cutpoints()$test reports the Lunn-McNeil test that the
ratio of cause-specific (or subdistribution) hazards differs between
groups, with the univariable omega+ ratio estimate and CI – the GCE test
that separate per-event Gray/log-rank tests miss. The per-event Gray
tests are also returned ($individual_gray) so the
difference is visible, along with a selection-adjusted permutation
p-value ($p.perm, via nperm > 0).gce_alligator()).
Draws the cumulative incidence of the event of interest and the
competing event on a single panel, with risk groups distinguished by
color and event types by line type (event of interest solid, competing
event dashed). By default it forms 2 groups at the optimal omega+
cutpoint (method = "optimal", log-scale;
criterion selectable) and annotates the cutpoint value and
group sizes. The x-axis is truncated at the readout time t
with 1-unit ticks; the y-axis auto-caps at the smallest of 25/50/75/100%
covering the peak incidence, with 5%/10% ticks (override with
ymax), so low-incidence plots are not compressed; and
numbers at risk are shown in an aligned panel below the plot (via
patchwork). Groups of 2-3 are labeled Low/High or
Low/Intermediate/High by default.hn,
prostate). Two competing-risks example datasets
show opposite GCE regimes: hn (head and neck, n = 1000) is
recurrence-dominant (omega+ > 1), while prostate (n =
1000) is competing-death-dominant (omega+ < 1). Both code
status as 0 = censored, 1 = event of interest, 2 =
competing, with time in years and sequential subject IDs.
They are provided for illustration only and do not represent real
patients. Categorical covariates (stage, site, Gleason, performance
status) are factors. All examples and the vignette use these
datasets.Refocused the package on modeling omega+ (event of interest vs. competing event).
$omega) for comparing model-predicted
with observed values.gce_tidy() no longer takes which (omega+
only). gce_cutpoints() now targets omega+ separation rather
than a log-rank split.lunnmcneil()
gains type = "finegray": a stacked, risk-set-weighted
(survival::finegray) model giving omega+ CIs and p-values
on the subdistribution scale from one joint model, replacing the
bootstrap in gcefg().gce_cutpoints(method = "optimal") finds cutpoints
maximizing the between-group omega+ separation (greedy maximum
between-group sum of squares).gce_alligator() — cumulative incidence of the event of
interest and the competing event within risk groups, with Gray’s
test.gce_calibration() — predicted vs. observed omega+ (or
omega) by rank group, user-specified number of groups.First release. gcemod is a successor to
gcerisk.
gceriskVariance for the omega ratio. Where
gcerisk estimates the variance of
log(omega ratio) by sampling independent normal deviates
for the cause-specific and all-cause coefficients and taking the
variance of their difference, gcemod uses a case-resampling
(nonparametric) bootstrap that refits the models on resampled subjects.
(This omega-ratio regression was removed in 0.2.0.)
Hypothesis testing added (Lunn-McNeil). The
omega-plus ratio (event of interest vs. competing event) is estimated
from a single stacked/augmented Cox model in the style of Lunn &
McNeil (1995). The event-type interaction coefficient is exactly
log(omega+ ratio) and carries a model-based (robust)
standard error, a per-covariate Wald test, and an omnibus multivariate
Wald test of no differential covariate effect across event types.
Because omega is a monotone transform of omega-plus
(omega = omega+/(1 + omega+)), this test applies to
both.
Cumulative hazard at time t.
gcerisk estimates cumulative hazards by regressing
-log(S(t)) on time with lm() and
extrapolating, which forces a linear cumulative hazard.
gcemod reads the Nelson-Aalen cumulative hazard
(cause-specific model) or the subdistribution cumulative hazard derived
from the Fine-Gray CIF (-log(1 - F(t))) directly at the
requested time via step functions.
Edge cases. Formula construction now uses
reformulate() instead of manual string pasting, which
failed for a single covariate (the 1:(p-1) loop). NA
handling (complete-case with a warning), single-covariate models, and
non-convergent bootstrap replicates (skipped with a warning) are
handled.
lunnmcneil() exposes the stacked-model estimator and
tests directly.gce_cutpoints() derives risk-stratification groups from
a GCE risk score (quantile groups, or an optimal two-group split by
log-rank scan).gce_performance() reports the cause-specific
concordance (C-index) and a grouped calibration table for the event of
interest.print(), summary(), and
gce_tidy() methods for compact, report-ready output.