Package {cifmodeling}

Title:	Visualization and Polytomous Modeling of Survival and Competing Risks
Version:	1.2.3
Description:	A publication-ready toolkit for modern survival and competing risks analysis with a minimal, formula-based interface. Both nonparametric estimation and direct polytomous regression of cumulative incidence functions (CIFs) are supported. The main functions 'cifcurve()', 'cifplot()', and 'cifpanel()' estimate survival and CIF curves and produce high-quality graphics with risk tables, censoring and competing-risk marks, and multi-panel or inset layouts built on 'ggplot2' and 'ggsurvfit'. The modeling function 'polyreg()' performs direct polytomous regression for coherent joint modeling of all cause-specific CIFs to estimate risk ratios, odds ratios, or subdistribution hazard ratios at user-specified time points. The function 'ciftest()' provides log-rank, Gray, robust score, augmented score, and multiple-direction tests for survival and competing-risk outcomes. All core functions adopt a formula-and-data syntax and return tidy and extensible outputs that integrate smoothly with 'modelsummary', 'broom', and the broader 'tidyverse' ecosystem. Key numerical routines are implemented in C++ via 'Rcpp'.
License:	MIT + file LICENSE
Encoding:	UTF-8
RoxygenNote:	7.3.3
RdMacros:	lifecycle
Config/Needs/website:	lifecycle
Depends:	R (≥ 4.1.0)
Suggests:	survival, mets, numDeriv, modelsummary, gtsummary (≥ 2.6.1), broom, broom.helpers (≥ 1.20.0), cardx (≥ 0.3.4), glue, gt, flextable, DiagrammeR, knitr, rmarkdown, testthat (≥ 3.0.0), pkgdown
Config/testthat/edition:	3
LinkingTo:	Rcpp
Imports:	Rcpp, nleqslv, boot, ggsurvfit, ggplot2, patchwork, scales, generics, lifecycle
VignetteBuilder:	knitr
BuildManual:	true
URL:	https://gestimation.github.io/cifmodeling/, https://github.com/gestimation/cifmodeling
BugReports:	https://github.com/gestimation/cifmodeling/issues
NeedsCompilation:	yes
Packaged:	2026-11-10 09:25:36 UTC; Shiro
Author:	Shiro Tanaka [aut, cre, cph], Shigetaka Kobari [ctb], Chisato Honda [ctb]
Maintainer:	Shiro Tanaka <gestimation@gmail.com>
Repository:	CRAN
Date/Publication:	2026-11-11 08:30:08 UTC

cifmodeling: Visualization and Polytomous Modeling of Survival and Competing Risks

Description

A publication-ready toolkit for modern survival and competing risks analysis with a minimal, formula-based interface. Both nonparametric estimation and direct polytomous regression of cumulative incidence functions (CIFs) are supported. The main functions 'cifcurve()', 'cifplot()', and 'cifpanel()' estimate survival and CIF curves and produce high-quality graphics with risk tables, censoring and competing-risk marks, and multi-panel or inset layouts built on 'ggplot2' and 'ggsurvfit'. The modeling function 'polyreg()' performs direct polytomous regression for coherent joint modeling of all cause-specific CIFs to estimate risk ratios, odds ratios, or subdistribution hazard ratios at user-specified time points. The function 'ciftest()' provides log-rank, Gray, robust score, augmented score, and multiple-direction tests for survival and competing-risk outcomes. All core functions adopt a formula-and-data syntax and return tidy and extensible outputs that integrate smoothly with 'modelsummary', 'broom', and the broader 'tidyverse' ecosystem. Key numerical routines are implemented in C++ via 'Rcpp'.

Author(s)

Maintainer: Shiro Tanaka gestimation@gmail.com (ORCID) [copyright holder]

Other contributors:

• Shigetaka Kobari [contributor]

• Chisato Honda [contributor]

See Also

Useful links:

• https://gestimation.github.io/cifmodeling/

• https://github.com/gestimation/cifmodeling

• Report bugs at https://github.com/gestimation/cifmodeling/issues

Create a survival or competing-risks response

Description

A lightweight response constructor used in cifcurve() and polyreg() to pass survival and competing-risks data via a model formula.

Usage

Event(time, event, allowed = getOption("cifmodeling.allowed", c(0, 1, 2)))

Arguments

Value

An object of class "Event" (a 2-column matrix) with columns time, event.

Lifecycle

See Also

polyreg() for log-odds product modeling of CIFs; cifcurve() for KM/AJ estimators; cifplot() for display of a CIF; cifpanel() for display of multiple CIFs; ggsurvfit::ggsurvfit, patchwork::patchwork and modelsummary::modelsummary for display helpers.

Examples

## event: 0=censor, 1=primary, 2=competing
data(diabetes.complications)
output <- polyreg(
  nuisance.model = Event(t, epsilon) ~ +1,
  exposure = "fruitq1",
  data = diabetes.complications,
  effect.measure1 = "RR",
  effect.measure2 = "RR",
  time.point = 8,
  outcome.type = "competing-risk"
)

Extract a ggplot for customization or export

Description

Extract a ggplot for customization or export

Usage

as_ggplot(x, ...)

Arguments

Value

A ggplot. This function does not estimate or predict anything.

Convert regression results to modelsummary input

Description

Convert regression results to modelsummary input

Usage

as_modelsummary(x, ...)

## S3 method for class 'polyreg'
as_modelsummary(
  x,
  event = c("both", "event1", "event2"),
  conf.int = x$conf.int,
  variance = NULL,
  ...
)

Arguments

Value

A modelsummary_list, or a named list of these for multiple outcomes. Each contains tidy and glance data frames. Extraction does not require modelsummary. Inference settings are recorded in glance and metadata.

Extract prediction data in the common long format

Description

Extract prediction data in the common long format

Usage

as_prediction_data(x, ...)

## S3 method for class 'polyreg'
as_prediction_data(
  x,
  exposure.level = x$estimand$code.exposure.ref,
  event = c("both", "event1", "event2"),
  ...
)

Arguments

Value

A data frame with .row, estimate, optional category, and optional conf.low/conf.high columns. No newdata prediction is added to polyreg.

Build the closed-form competing-risk augmentation

Description

This computes the working-model process

H_a(t,X)=\int_t^\tau a(u)(X-e_X(u))G_C(u)d\Lambda_1(u|X)

in each exposure-by-competing-risk-by-censoring nuisance cell, then evaluates the centered censoring-martingale projection. The returned matrix is on the score scale; its column sums are the closed-form augmentation score.

Usage

build_closed_form_augmentation(
  base,
  working,
  t,
  epsilon,
  x,
  exposure,
  strata.censor = rep.int("pooled", length(t)),
  strata.competing.risk = rep.int("pooled", length(t)),
  weights = rep.int(1, length(t)),
  code.event1 = 1L,
  code.event2 = 2L,
  code.censoring = 0L,
  prob.bound = 1e-07,
  strata.event = rep.int("pooled", length(t))
)

Build individual Fine-Gray score-scale contributions

Description

The base component is the known-censoring-distribution martingale score. The censoring component is the delta-method contribution from estimating the stratified censoring hazards. Their column sums equal the plug-in score and zero, respectively. The matrices are diagnostic infrastructure for the augmented branch; they do not replace Gray's classical covariance.

Usage

build_fg_score_iid(
  t,
  epsilon,
  x,
  code.event1 = 1L,
  code.event2 = 2L,
  code.censoring = 0L,
  strata = rep.int("pooled", length(t)),
  weights = rep.int(1, length(t)),
  rho = 0,
  gamma = 0,
  fh.weight = NULL,
  censoring = NULL,
  prob.bound = 1e-07,
  prob.truncation = NULL,
  strata.event = rep.int("pooled", length(t))
)

Build a finite-iteration time-directed Fine-Gray score

Description

iterations = 0 is intentionally handled by the closed-form caller. This function computes one or more empirical Appendix-E refinements and maps the resulting full-data direction to observed-data score contributions.

Usage

build_iterated_score_iid(
  base,
  working,
  t,
  epsilon,
  x,
  exposure,
  weights = rep.int(1, length(t)),
  iterations = 1L,
  code.event1 = 1L,
  code.event2 = 2L,
  code.censoring = 0L,
  prob.bound = 1e-07,
  tolerance = NULL
)

Internal Kaplan-Meier interface

Description

The numerical engine always retains fractional weighted counts. The public package interface defaults to integer reporting for plotting compatibility; internal estimation code can request the unrounded weighted totals.

Usage

calculateKM(
  t,
  d,
  w = numeric(),
  strata = integer(),
  error = "greenwood",
  count.type = c("integer", "numeric")
)

Arguments

A12 for log-rank-type score using directional finite differences with dw/dB from WeightIt

Description

Returns A12 on the "mean score" scale: (1/n) * dU_total/dB^T

Usage

calculate_A12_logrank_weightit(
  t,
  epsilon,
  strata,
  data,
  exposure,
  weightit,
  code.exposure.ref = NULL,
  prefix = "a",
  rho = 0,
  gamma = 0,
  prob.bound = 1e-07,
  fd_rel_step = 1e-06
)

Compute Gray's K-sample test

Description

Implements the predictable-weight estimating equations and covariance from Gray (1988). The pooled weight is ⁠(1 - F_1(t-))^rho * F_1(t-)^gamma⁠; gamma = 0 recovers the family used by cmprsk::cuminc().

Usage

calculate_gray(
  t,
  epsilon,
  exposure,
  code.exposure.ref = NULL,
  prefix = "a",
  weights = rep.int(1, length(t)),
  strata = rep.int(1L, length(t)),
  data,
  rho = 0,
  gamma = 0,
  prob.bound = 1e-10
)

Details

weights are frequency weights. strata defines independent analysis strata; stratum-specific scores and covariance matrices are summed.

Censoring-estimation correction to left-limit IPCW regression scores

Description

Returns q_i such that IF_beta_i = n H^-1 (s_i - mean(s) + q_i). Uses the same source cohort as KM, with unit sampling weights. The supplied scores already contain IPCW, including zero scores for unknown outcomes. This product-limit derivative includes the (1-h) factor at tied jumps. No n-by-n influence matrix is constructed: time-grouped score sums and cumulative risk contributions give O(n log n + np) work and O(np) storage.

Usage

censoring_km_score_correction(
  object,
  t,
  epsilon,
  strata,
  evaluation.time,
  score
)

These arguments are shared by cifplot(), cifpanel(), and cifcurve().

Description

These arguments are shared by cifplot(), cifpanel(), and cifcurve().

Arguments

These arguments are shared by cifplot() and cifpanel().

Description

These arguments are shared by cifplot() and cifpanel().

Arguments

Calculate the Kaplan-Meier estimator and the Aalen-Johansen estimator

Description

Core estimation routine that computes a survfit-compatible object from a formula + data interface (Event() or survival::Surv() on the LHS, and a stratification variable on the RHS if necessary). The back-end C++ routine supports both weighted and stratified data. Use this when you want numbers only (e.g. estimates, SEs, CIs and influence functions) and will plot it yourself.

Usage

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

cifcurve(
  formula,
  data,
  weights = NULL,
  n.risk.type = "weighted",
  subset.condition = NULL,
  na.action = na.omit,
  outcome.type = NULL,
  time.point = NULL,
  null.hypothesis = NULL,
  code.event1 = 1,
  code.event2 = 2,
  code.censoring = 0,
  error = NULL,
  conf.type = "arcsine-square root",
  conf.int = 0.95,
  report.influence.function = FALSE,
  report.survfit.std.err = FALSE,
  engine = "calculateAJ_Rcpp",
  prob.bound = 1e-07,
  conf.dist = c("normal", "t"),
  df.method = c("deletion", "equal")
)

Arguments

Details

Typical use cases

• When outcome.type = "survival", this is a thin wrapper around the KM estimator with the chosen variance / CI transformation.

• When outcome.type = "competing-risk", this computes the AJ estimator of CIF for code.event1. The returned ⁠$surv⁠ is 1 - CIF, i.e. in the format that ggsurvfit expects.

• Use cifplot() if you want to go straight to a figure; use cifcurve() if you only want the numbers.

Risk set display

• Set n.risk.type to control whether ⁠$n.risk⁠ reflects weighted, unweighted, or Kish effective sample size (ESS) counts. This only affects the reported counts (e.g., for plotting or debugging) and leaves estimates and SEs unchanged.

Standard error and confidence intervals

With the C++ engine, error = "fay" applies the Fay–Graubard-type sandwich correction with leverage capped at 0.75. Combining error = "if" or "fay" with conf.dist = "t" gives IF-t or FG-t, respectively. "fg" and "fay-graubard" are accepted aliases for "fay"; the returned error is always "fay" for this correction. All supplied weights are held fixed, including externally estimated IPTW; uncertainty from fitting a propensity-score model is not included. Deletion degrees of freedom use unscaled hazard-score equations and no-event-baseline hazard logits (a positive cause is the baseline at complete exit). Equal allocation is invariant to fixed invertible score transformations within this system; its small-sample coverage need not equal deletion allocation. These are working quadratic-form approximations, not exact Student laws.

With the C++ engine, IF and FG intervals use exact boundary singletons; interior zero variance, unavailable degrees of freedom, and times beyond positive-weight follow-up produce missing limits with reasons in ⁠$conf.status⁠. No normal fallback is used for a failed Student interval. Degrees of freedom vary by stratum and time; calculation can be expensive for large risk sets. The optional one-sided test retains its documented normal reference distribution. Delta and Aalen variances use all event times up to the reported time. If the Aalen formula gives a negative variance with tied events, its SE and confidence limits are missing. Returned influence functions are uncorrected case-weight derivatives (IF/n) of the event probability ⁠1 - $surv⁠, even when the SE uses FG-adjusted contributions.

One-sided normal approximation test

When both time.point and null.hypothesis are specified, cifcurve() compares the estimate at time.point with null.hypothesis using a normal approximation. The test statistic is

Z = (g(\hat{p}) - g(p_0)) / \widehat{SE}(g(\hat{p})).

where g() is the transformation specified by conf.type. For survival outcomes, \hat p is the Kaplan-Meier survival estimate and the default one-sided alternative is S(t) > S_0(t). For competing-risk outcomes, the reported estimate is the Aalen-Johansen cumulative incidence estimate F_1(t) for code.event1, but the transformation for the one-sided test is applied to 1 - F_1(t) and 1 - F_{10}(t), consistently with the survfit-compatible ⁠$surv⁠ component. Thus, conf.type = "log-log" corresponds to \log[-\log(1 - F_1(t))] for competing-risk outcomes. The default one-sided alternative is F_1(t) < F_{10}(t), equivalently 1 - F_1(t) > 1 - F_{10}(t). When strata are present, the test is calculated separately within each stratum.

Value

An object inheriting from "cifcurve" and "survfit". For outcome.type="survival", ⁠$surv⁠ is the survival function. For outcome.type="competing-risk", ⁠$surv⁠ equals 1 - CIF for code.event1. SE and CIs are provided per error, conf.type and conf.int. This enables an independent use of standard methods for survfit such as:

• summary(): time-by-time estimates with SEs and CIs

• plot(): base R stepwise survival/CIF curves

• mean(): restricted mean survival estimates with CIs

• quantile(): quantile estimates with CIs

If null.hypothesis is specified, the returned object additionally contains a one.sided.p list with the one-sided normal approximation test results. Its table marks strata beyond positive-weight follow-up with test.status = "outside_observed_followup".

Note that ⁠$n.risk⁠, ⁠$n.event⁠, and ⁠$n.censor⁠ are rounded up to the nearest integer regardless of whether the data is weighted or not. Some methods (e.g. residuals.survfit) may not be supported. IF/FG results also record conf.dist, df, df.reason, critical.value, and conf.status; Student results record df.method. These vectors follow the fitted curve's stratum/time ordering. FG results record the fixed leverage cap in fay.bound. Influence matrices have one row per original observation within each stratum, including zero-weight rows. summary(fit, times = ..., extend = TRUE) extracts the saved intervals, probability-scale SEs and time-specific inference diagnostics. Unsupported times retain missing limits. Use this summary rather than indexing added inference vectors on copies whose time grid was expanded by survfit0(). logse records the scale of the stored std.err. The internal reporting component preserves the native inference grid for this extraction.

Lifecycle

See Also

polyreg() for log-odds product modeling of CIFs; cifplot() for display of a CIF; cifpanel() for display of multiple CIFs; ggsurvfit::ggsurvfit, patchwork::patchwork and modelsummary::modelsummary for display helpers.

Examples

data(diabetes.complications)
output1 <- cifcurve(Event(t,epsilon) ~ fruitq,
                    data = diabetes.complications,
                    outcome.type="competing-risk")
cifplot(output1,
        outcome.type = "competing-risk",
        type.y = "risk",
        add.risktable = FALSE,
        label.y = "CIF of diabetic retinopathy",
        label.x = "Years from registration")

output2 <- cifcurve(Event(t,epsilon) ~ fruitq,
                    data = diabetes.complications,
                    outcome.type = "competing-risk",
                    time.point = 8,
                    null.hypothesis = 0.30)
output2$one.sided.p

Create a flowchart of exclusions, treatment groups, and outcome status

Description

cifflowchart() creates a lightweight flowchart of exclusions, treatment groups, and outcome status with DiagrammeR. It can be used as a simple CONSORT-like flowchart helper, but it is not intended to fully automate a CONSORT-compliant diagram.

Usage

cifflowchart(
  formula,
  data,
  time.point = NULL,
  withdraw.consent = NULL,
  ineligible = NULL,
  pre.exclude = NULL,
  post.exclude = NULL,
  subset.condition = NULL,
  na.action = na.pass,
  outcome.type = NULL,
  code.event1 = 1,
  code.event2 = 2,
  code.censoring = 0,
  label.strata = NULL,
  order.strata = NULL,
  label.events = NULL,
  title = NULL,
  percent = TRUE,
  digits = 1,
  ...
)

Arguments

Value

A DiagrammeR::grViz() htmlwidget.

Lifecycle

Examples

if (requireNamespace("DiagrammeR", quietly = TRUE)) {
  data(diabetes.complications)

  cifflowchart(
    Event(t, epsilon) ~ fruitq1,
    data = diabetes.complications,
    time.point = 8,
    outcome.type = "competing-risk"
  )
}

Arrange multiple survival and CIF plots in a panel display

Description

cifpanel() is the panel-building counterpart of cifplot(). It takes one or more model formulas (or, alternatively, one formula and several event-coding specifications) and returns a multi-panel figure, typically as a patchwork-compatible object. Most display options (axis labels, marks, style, ggsave options) are shared with cifplot(), but per-panel legends and number-at-risk tables are suppressed to avoid duplicated display. Typical use cases are:

• Compare CIF (event 1) vs CIF (event 2) in a 1×2 layout.

• Compare survival/CIF curves across strata with a shared legend and matched axes.

• Display a plot with an enlarged y-axis inside a full-scale plot.

Usage

cifpanel(
  plots = NULL,
  formula = NULL,
  formulas = NULL,
  data = NULL,
  weights = NULL,
  subset.condition = NULL,
  na.action = na.omit,
  outcome.type = NULL,
  code.event1 = 1,
  code.event2 = 2,
  code.censoring = 0,
  code.events = NULL,
  error = NULL,
  conf.type = NULL,
  conf.int = NULL,
  n.risk.type = c("weighted", "unweighted", "ess"),
  type.y = NULL,
  label.x = NULL,
  label.y = NULL,
  label.strata = NULL,
  order.strata = NULL,
  level.strata = NULL,
  limits.x = NULL,
  limits.y = NULL,
  breaks.x = NULL,
  breaks.y = NULL,
  add.conf = NULL,
  add.risktable = NULL,
  add.estimate.table = NULL,
  symbol.risk.table = NULL,
  font.size.risk.table = NULL,
  add.censor.mark = NULL,
  shape.censor.mark = NULL,
  size.censor.mark = NULL,
  add.competing.risk.mark = NULL,
  competing.risk.time = NULL,
  shape.competing.risk.mark = NULL,
  size.competing.risk.mark = NULL,
  add.intercurrent.event.mark = NULL,
  intercurrent.event.time = NULL,
  shape.intercurrent.event.mark = NULL,
  size.intercurrent.event.mark = NULL,
  add.quantile = NULL,
  level.quantile = NULL,
  rows.columns.panel = c(1, 1),
  inset.panel = FALSE,
  title.panel = NULL,
  subtitle.panel = NULL,
  caption.panel = NULL,
  tag.panel = NULL,
  title.plot = NULL,
  style = "classic",
  palette = NULL,
  linewidth = 0.8,
  linetype = FALSE,
  font.family = "sans",
  font.size = 8,
  legend.position = "top",
  legend.collect = TRUE,
  inset.left = 0.6,
  inset.bottom = 0.05,
  inset.right = 0.98,
  inset.top = 0.45,
  inset.align.to = c("panel", "plot", "full"),
  inset.legend.position = NULL,
  print.panel = FALSE,
  filename.ggsave = NULL,
  width.ggsave = NULL,
  height.ggsave = NULL,
  dpi.ggsave = 300,
  survfit.info = NULL,
  axis.info = NULL,
  visual.info = NULL,
  panel.info = NULL,
  style.info = NULL,
  inset.info = NULL,
  print.info = NULL,
  ggsave.info = NULL,
  engine = "cifplot",
  ...
)

Arguments

Details

Overview

cifpanel() composes multiple survival/CIF plots into a single figure. For each panel, it estimates curves via cifcurve() and renders them with cifplot(). You can supply a single formula reused across panels or a list in formulas (one per panel). When both are provided, formulas wins.

Outcome type & event coding

• Use outcome.type to set per-panel estimator ("survival"=KM, "competing-risk"=AJ).

• Alternatively, pass code.events per panel to infer the type:

  • length 2 = survival: c(event1, censor)

  • length 3 = competing-risk: c(event1, event2, censor)

• If code.events is NULL, code.event1, code.event2, code.censoring are combined into code.events = list(c(code.event1, code.event2, code.censoring)) with NA values dropped.

• If outcome.type is NULL, the function infers each panel from its code.events[[i]] length. When both are given, outcome.type takes precedence.

• Control risk-set displays via n.risk.type, which is recycled per panel and forwarded to cifcurve() to decide which risk set size populates ⁠$n.risk⁠ (e.g., weighted vs. unweighted counts).

Panel-wise vs shared arguments

Panel layout is specified by length-2 vector rows.columns.panel. This function can also automatically determine the panel count in the following order: (1) if plots is supplied, its length defines the number of plots, (2) else if formulas is supplied, its length defines the number of plots, (3) else if code.events is supplied, its length defines the number of plots together with formula, and (4) otherwise rows.columns.panel=c(1,1).

Many arguments accept a scalar (recycled to all panels) or a list/vector (one entry per panel). Precedence: panel-wise explicit values > shared scalar > internal defaults. Length-1 inputs are recycled.

Grid vs inset composition

• Grid mode (inset.panel = FALSE, default): plots are arranged with patchwork::wrap_plots() and plot_layout(). If legend.collect = TRUE, legends are collected across panels where possible.

• Inset mode (inset.panel = TRUE): the second plot is overlaid into the first using patchwork::inset_element(). Only the first two plots are used; extra plots are ignored. Control the inset box with inset.left, inset.bottom, inset.right, inset.top, and its reference frame via inset.align.to ("panel", "plot", or "full").

Advanced panel controls (forwarded to cifplot())

The following arguments allow per-panel control by supplying vectors/lists, or shared control by supplying scalars. They are forwarded to cifplot().

• formula or formulas: one formula or a list of formulas; each entry creates a panel.

• data, outcome.type, code.events, type.y: recycled across panels unless a list is supplied for per-panel control.

• rows.columns.panel: specification of grid layout by c(rows, cols).

• inset.panel: inset layout.

• title.panel, subtitle.panel, caption.panel, title.plot: overall titles and captions.

• tag.panel: panel tag style (e.g., "A", "a", "1").

• label.x, label.y, limits.x, limits.y, breaks.x, breaks.y: shared axis control unless a list is supplied for per-panel control.

Scale & labels

Plot layers (toggles)

(Time marks inputs such as competing.risk.time / intercurrent.event.time can be given via ... if needed; names must match strata labels.)

Legend & annotations

• legend.position: "top", "right", "bottom", "left", or "none" (applies to all panels).

• Grid mode: legend.collect = TRUE attempts a shared legend.

• Panel annotations: title.panel, subtitle.panel, caption.panel.

• Tagging: tag.panel is passed to patchwork::plot_annotation().

• In inset mode, title.plot = c(title_base, title_inset) labels the two plots.

Export (optional)

If filename.ggsave is non-NULL, the composed panel is saved with ggsave() using width.ggsave, height.ggsave, and dpi.ggsave. Otherwise, the function returns objects without saving.

Notes

• Mixed panel types are supported (e.g., AJ in panel 1; KM in panel 2).

• If formulas is shorter than the grid capacity, empty slots are ignored.

• When supplying vectors/lists per panel, their lengths must match the number of panels; length-1 inputs are recycled; otherwise an error is thrown.

• For CIF displays, set type.y = "risk". For survival scale, use type.y = NULL or ⁠= "surv"⁠. For ADaM-style data, use code.events=c(0,1) or code.event1 = 0, code.censoring = 1.

• Additional graphical options (e.g., theme) can be added post-hoc to each element of list.plot or to the composed patchwork.

Value

A "cifpanel" object (returned invisibly), which is a list with at least the following elements:

• list.plot: a list of ggplot objects, one per panel

• patchwork: a patchwork object representing the composed panel

• plot: reserved for backwards compatibility (always NULL)

• metadata fields mirroring those in cifplot() (such as information on the fitted curves and display settings)

When print.panel = TRUE, the patchwork object is printed in interactive sessions in addition to being returned.

Lifecycle

See Also

polyreg() for log-odds product modeling of CIFs; cifcurve() for KM/AJ estimators; cifplot() for display of a CIF; ggsurvfit::ggsurvfit, patchwork::patchwork and modelsummary::modelsummary for display helpers.

Examples

data(diabetes.complications)
output1 <- cifpanel(
  title.panel = "A comparison of cumulative incidence of competing events",
  rows.columns.panel = c(1,2),
  formula = Event(t, epsilon) ~ fruitq,
  data = diabetes.complications,
  outcome.type = "competing-risk",
  code.events = list(c(1,2,0), c(2,1,0)),
  label.y = c("Diabetic retinopathy", "Macrovascular complications"),
  label.x = "Years from registration",
  subtitle.panel = "Stratified by fruit intake",
  caption.panel  = "Data: diabetes.complications",
  title.plot = c("Diabetic retinopathy", "Macrovascular complications"),
  legend.position = "bottom",
  legend.collect=TRUE
)
print(output1)

output2 <- cifplot(Event(t,epsilon) ~ fruitq,
                   data = diabetes.complications,
                   outcome.type="competing-risk",
                   code.event1=2,
                   code.event2=1,
                   add.conf = FALSE,
                   add.risktable = FALSE,
                   label.y="CIF of macrovascular complications",
                   label.x="Years from registration")
output3 <- cifplot(Event(t,epsilon) ~ fruitq,
                   data = diabetes.complications,
                   outcome.type="competing-risk",
                   code.event1=2,
                   code.event2=1,
                   add.conf = FALSE,
                   add.risktable = FALSE,
                   label.y="",
                   label.x="",
                   limits.y=c(0,0.15))
output4 <- list(a = output2$plot, b = output3$plot)
output5 <- cifpanel(plots = output4,
         inset.panel = TRUE,
         inset.left = 0.40, inset.bottom = 0.45,
         inset.right = 1.00, inset.top = 0.95,
         inset.align.to = "plot",
         inset.legend.position = "none",
         legend.position = "bottom")
print(output5)

Generate a survival/CIF curve with marks that represent censoring, competing risks and intercurrent events

Description

This function generates a survival or CIF curve from a unified formula–data interface or from an existing survfit object. When a formula is supplied, the LHS is typically Event() or survival::Surv(), and the RHS specifies an optional stratification variable. In addition to the curves themselves, cifplot() can add numbers-at-risk tables, tables of point estimates and CIs, censoring marks, competing-risk marks, and intercurrent-event marks.

For usual single-panel mode, the function returns an object whose plot component is a regular ggplot object that can be further modified. Use as_ggplot(result) + ggplot2::labs(...) or result$plot + ...; the wrapper itself is not a ggplot. For more complex multi-panel displays, cifplot() can internally call cifpanel() via several “panel modes” (per event, per variable, or censoring-focused).

Usage

cifplot(
  formula_or_fit,
  data = NULL,
  weights = NULL,
  subset.condition = NULL,
  na.action = na.omit,
  outcome.type = "competing-risk",
  code.event1 = 1,
  code.event2 = 2,
  code.censoring = 0,
  code.events = NULL,
  error = NULL,
  conf.type = "arcsine-square root",
  conf.int = 0.95,
  n.risk.type = c("weighted", "unweighted", "ess"),
  type.y = NULL,
  label.x = "Time",
  label.y = NULL,
  label.strata = NULL,
  level.strata = NULL,
  order.strata = NULL,
  limits.x = NULL,
  limits.y = NULL,
  breaks.x = NULL,
  breaks.y = NULL,
  use.coord.cartesian = FALSE,
  add.conf = TRUE,
  add.risktable = TRUE,
  add.estimate.table = FALSE,
  symbol.risk.table = "square",
  font.size.risk.table = 3,
  add.censor.mark = TRUE,
  shape.censor.mark = 3,
  size.censor.mark = 2,
  add.competing.risk.mark = FALSE,
  competing.risk.time = list(),
  shape.competing.risk.mark = 16,
  size.competing.risk.mark = 2,
  add.intercurrent.event.mark = FALSE,
  intercurrent.event.time = list(),
  shape.intercurrent.event.mark = 1,
  size.intercurrent.event.mark = 2,
  add.quantile = FALSE,
  level.quantile = 0.5,
  panel.per.event = FALSE,
  panel.censoring = FALSE,
  panel.per.variable = FALSE,
  panel.mode = "auto",
  rows.columns.panel = NULL,
  style = "classic",
  palette = NULL,
  linewidth = 0.8,
  linetype = FALSE,
  font.family = "sans",
  font.size = 12,
  legend.position = "top",
  print.panel = FALSE,
  filename.ggsave = NULL,
  width.ggsave = 6,
  height.ggsave = 6,
  dpi.ggsave = 300,
  survfit.info = NULL,
  axis.info = NULL,
  visual.info = NULL,
  panel.info = NULL,
  style.info = NULL,
  inset.info = NULL,
  print.info = NULL,
  ggsave.info = NULL,
  ...
)

Arguments

Details

Typical use cases

• Draw one survival/CIF curve set by exposure groups (e.g., treatment vs control).

• Call cifpanel() with a simplified code to create a panel displaying plots of multiple stratified survival/CIF curves or CIF curves for each event type.

• Add CIs and censor/competing-risk/intercurrent-event marks.

• Add number-at-risk table to display the number at risk or the estimated survival probabilities or CIFs and CIs at each point in time.

Key arguments shared with cifcurve()

• Outcome type and estimator

  • outcome.type = "survival": Kaplan-Meier estimator

  • outcome.type = "competing-risk": Aalen-Johansen estimator

• Confidence intervals

  • conf.int sets the two-sided level (default 0.95)

  • conf.type chooses the transformation ("arcsine-square root", "plain", "log", "log-log", "logit", or "none")

  • error chooses the estimator for SE ("greenwood", "tsiatis" or "if" for survival curves and "delta", "aalen" or "if" for CIFs)

• Risk sets used in tables

  • n.risk.type controls whether ⁠$n.risk⁠ reflects weighted, unweighted, or effective sample size counts when building risk tables (forwarded to cifcurve()); when a fitted survfit object is supplied, existing risk sets are used as-is.

Key arguments for cifplot()

• Data visualization

  • add.conf adds CIs on the ggplot2-based plot

  • add.competing.risk.mark and add.intercurrent.event.mark adds symbols to describe competing risks or intercurrent events in addition to conventional censoring marks with add.censor.mark

  • add.risktable adds numbers at risk

  • add.estimate.table adds time-by-time estimates and CIs

  • add.quantile adds a reference line at a chosen quantile level

• Plot customization

  • type.y chooses y-axis ("surv" for survival and "risk" for 1-survival/CIF)

  • limits.x, limits.y, breaks.x, breaks.y: numeric vectors for axis control

  • style specifies the appearance of plot ("classic", "bold", "framed", "grid", "gray" or "ggsurvfit")

  • palette specifies color of each curve (e.g. ⁠palette=c("blue1", "cyan3", "navy", "deepskyblue3"))⁠)

• Panel display

  • panel.per.variable produces multiple survival/CIF curves per stratification variable specified in the formula

  • panel.per.event produces CIF curves for each event type

  • panel.censoring produces the Kaplan–Meier curves for (event, censor) and (censor, event) so that censoring patterns can be inspected

  • panel.mode uses automatic panel mode

When panel.per.event = TRUE, two panels are created with code.events = list(c(e1, e2, c), c(e2, e1, c)), where code.events = c(e1, e2, c) is the input coding for event1, event2, and censoring. Common legend is collected by default (legend.collect = TRUE).

Numeric stratification variables are normalized automatically. Columns with fewer than nine distinct numeric values are coerced to factors; columns with nine or more distinct numeric values are split at the median into “Below median” and “Above median” strata.

Advanced control not required for typical use

The arguments below fine-tune internal estimation and figure appearance. Most users do not need to change these defaults.

Graphical layers

Time for marks

Appearance of marks

Panel display

Axes and legend

Export

Notes

• For CIF displays, set type.y = "risk". For survival scale, use type.y = NULL or ⁠= "surv"⁠. For a cumulative hazard plot, use type.y = "cumhaz". To generate a log-log plot, use type.y = "cloglog".

• Event coding can be controlled via code.event1, code.event2, code.censoring. For ADaM-style data, use code.event1 = 0, code.censoring = 1.

• Per-stratum time lists should have names identical to plotted strata labels.

Value

A "cifplot" object (a list) with at least the following elements:

• plot: a ggplot object containing the main plot

• patchwork: reserved for compatibility with panel displays (typically NULL for single-panel plots)

• survfit.info, axis.info, visual.info, panel.info, style.info, inset.info, print.info, ggsave.info: internal lists storing the fitted curves and display settings

• version: a character string giving the cifmodeling version used

• call: the original function call

When a panel mode is active and print.panel = TRUE, the panel is also printed in interactive sessions.

Lifecycle

See Also

polyreg() for log-odds product modeling of CIFs; cifcurve() for KM/AJ estimators; cifpanel() for display of multiple CIFs; ggsurvfit::ggsurvfit, patchwork::patchwork and modelsummary::modelsummary for display helpers.

Examples

data(diabetes.complications)
cifplot(Event(t,epsilon) ~ fruitq,
        data = diabetes.complications,
        outcome.type="competing-risk",
        add.risktable = FALSE,
        label.y='CIF of diabetic retinopathy',
        label.x='Years from registration')

Tests for survival and cumulative incidence curves

Description

Provides a formula-based interface for log-rank, Fleming-Harrington, Gray, and augmented score tests. test = "auto" selects log-rank for survival outcomes and the augmented score test for competing-risk outcomes. The "early" and "late" presets select that same outcome-specific test with Fleming–Harrington parameters ⁠(rho, gamma) = (1, 0)⁠ and ⁠(0, 1)⁠, respectively. Analysis strata can be specified independently with strata. For augmented tests, working Aalen-Johansen distributions are estimated within exposure-by-strata.competing.risk cells, censoring distributions are estimated within strata.censor, and the closed-form augmentation is used. Setting iteration to a positive integer applies that many Scheike fixed-point refinements. Finite iterates are anchored as the closed-form score plus the change in the AIPWCC image from the seed to the refined direction, preserving the zeroth-step identity in the observed sample.

Usage

ciftest(
  formula,
  data,
  weights = NULL,
  subset.condition = NULL,
  na.action = na.omit,
  outcome.type = "auto",
  test = "auto",
  code.event1 = 1,
  code.event2 = 2,
  code.censoring = 0,
  rho = NULL,
  gamma = NULL,
  iteration = 0L,
  tolerance = NULL,
  strata = NULL,
  strata.censor = NULL,
  strata.competing.risk = NULL,
  tau = NULL,
  prob.bound = 1e-07,
  prob.truncation = NULL,
  ...
)

Arguments

Value

An object inheriting from "ciftest" and "htest", with a test-specific subclass. An ordinary unweighted log-rank result also inherits from "survdiff" and contains a fully compatible object in survdiff. Every score-based branch returns an analysis-row by contrast matrix of individual null-score contributions in score.iid. Augmented results include score.iid.base, score.iid.censor, and score.iid.augment matrices and use their summed empirical cross-product as vcov.score. Standard Gray results retain the classical Gray covariance. If the optional standard-Gray score-IID diagnostic cannot be constructed, the Gray test is still returned; its score-IID matrices contain NA and the reason is recorded in diagnostics$score.iid.error. Positive iteration values return the requested finite-iteration score, its subject-level score matrix, anchoring diagnostics, and fixed-point diagnostics. The event-time FH process actually used by the fit is returned in diagnostics$fh.weight.process, including whether it came from the pooled or event-stratified AJ left limit or the native Gray subdistribution construction.

Build a finite-iteration time-directed Fine-Gray score

Description

iterations = 0 is intentionally handled by the closed-form caller. This function computes one or more empirical Appendix-E refinements and maps the resulting full-data direction to observed-data score contributions.

Usage

ciftest_iteration_setup(
  base,
  working,
  t,
  epsilon,
  x,
  exposure,
  weights = rep.int(1, length(t)),
  code.event1 = 1L,
  code.event2 = 2L,
  code.censoring = 0L,
  prob.bound = 1e-07,
  prob.truncation = NULL
)

Multiple-direction quadratic score test

Description

Combines several Fleming–Harrington score directions in the quadratic form proposed by Ditzhaus and Friedrich. Each component is fitted with the same outcome-specific ciftest() branch. Classical log-rank components use the joint, finite-population-corrected hypergeometric covariance, with cross-direction blocks weighted by w_r(t) * w_s(t). Classical Gray components analogously use the joint Gray (1988) covariance, including the target-event, competing-event, and censoring contributions. Score and augmented components use the cross-product of the stacked individual score contributions. The default chi-square reference uses the numerical rank of the joint covariance. This function does not perform the permutation calibration from the original two-sample survival proposal.

Usage

ciftest_mdir(
  formula,
  data,
  directions = c("early", "late", "unweighted"),
  weights = NULL,
  subset.condition = NULL,
  outcome.type = "auto",
  test = "auto",
  code.event1 = 1,
  code.event2 = 2,
  code.censoring = 0,
  iteration = 0L,
  tolerance = NULL,
  strata = NULL,
  strata.censor = NULL,
  strata.competing.risk = NULL,
  tau = NULL,
  prob.bound = 1e-07,
  prob.truncation = NULL,
  na.action = stats::na.omit
)

Arguments

Value

An object inheriting from "ciftest_mdir", "ciftest", and "htest". Component fits are available in components; the unmodified individual contributions are in score.iid, and the covariance actually used by the quadratic statistic is in vcov.score.

References

Ditzhaus M, Friedrich S. More powerful logrank permutation tests for two-sample survival data. arXiv preprint (2018).

Data from a cohort study of patients with type 2 diabetes

Description

Anonymized data from a cohort study of patients with type 2 diabetes followed for ocular and macro-vascular complications.

Usage

data(diabetes.complications)

Format

A data frame with 978 observations and 19 variables:

t

Follow-up time in years.

epsilon

Event type indicator (0 = censored, 1 = diabetic retinopathy, 2 = macro-vascular complication).

fruit

Fruit intake (g/day).

fruitq

Quartile of fruit intake.

fruitq1

Binary indicator for low fruit intake.

strata

Stratum used for inverse probability of censoring weights.

age

Age at baseline (years).

sex

Sex coded as 0 = woman, 1 = man.

bmi

Body mass index at baseline.

hba1c

Hemoglobin A1c (%).

diabetes_duration

Duration of diabetes (years).

drug_oha

Indicator for oral hypoglycemic agent use.

drug_insulin

Indicator for insulin use.

sbp

Systolic blood pressure (mmHg).

ldl

Low-density lipoprotein cholesterol (mg/dL).

hdl

High-density lipoprotein cholesterol (mg/dL).

tg

Triglycerides (mg/dL).

current_smoker

Indicator for current smoking status.

alcohol_drinker

Indicator for current alcohol drinking.

ltpa

Leisure-time physical activity (METs).

Details

The variables include follow-up time, cause-specific event indicators, exposure indicators for fruit intake, censoring strata, and a set of covariates used in the package vignettes.

Source

Anonymized example data from the Japan Diabetes Complications Study. The original cohort analysis of fruit intake and incident diabetic retinopathy was reported by Tanaka S, Yoshimura Y, Kawasaki R, Kamada C, Tanaka S, Horikawa C, Ohashi Y, Araki A, Ito H, Akanuma Y, Yamada N, Yamashita H, Sone H, for the Japan Diabetes Complications Study Group. Fruit intake and incident diabetic retinopathy with type 2 diabetes. Epidemiology. 2013, 24(2), 204–211. doi:10.1097/EDE.0b013e318281725e

Examples

data(diabetes.complications)
str(diabetes.complications)

Estimate a stratified censoring Kaplan-Meier distribution

Description

Uses the package Kaplan-Meier counts, retaining both left and right limits. The default keeps the existing censoring-first convention for ciftest. With tie.order="event_first", actual failures leave before the censoring jump, consistently with observation of failures when T <= C.

Usage

estimate_censoring_km(
  t,
  epsilon,
  code.censoring = 0L,
  strata = rep.int("pooled", length(t)),
  weights = rep.int(1, length(t)),
  prob.bound = 1e-07,
  prob.truncation = NULL,
  censoring.event = epsilon == code.censoring,
  tie.order = c("censoring_first", "event_first")
)

Estimate working Aalen-Johansen distributions by exposure and stratum

Description

The package Aalen-Johansen engine is evaluated separately within every observed exposure-by-stratum cell. A second call with the two event codes exchanged supplies the competing-event cumulative incidence.

Usage

estimate_working_aj(
  t,
  epsilon,
  exposure,
  strata = rep.int("pooled", length(t)),
  weights = rep.int(1, length(t)),
  code.event1 = 1L,
  code.event2 = 2L,
  code.censoring = 0L,
  prob.bound = 1e-07,
  prob.truncation = NULL
)

Extract per-stratum event times from a formula and data

Description

Creates a list of event times that can be passed to downstream visualization or analysis functions such as competing.risk.time or intercurrent.event.time in cifplot() and cifpanel(). Event types are specified by event 1, event 2, censoring, or user-specified codes.

Usage

extract_time_to_event(
  formula,
  data,
  weights = NULL,
  subset.condition = NULL,
  na.action = na.omit,
  which.event = c("event2", "event1", "censor", "censoring", "user_specified"),
  code.event1 = 1,
  code.event2 = 2,
  code.censoring = 0,
  code.user.specified = NULL,
  read.unique.time = TRUE,
  drop.empty = TRUE
)

Arguments

Details

This function is typically used internally by plotting and model functions, but can also be called directly to inspect the per-stratum event-time structure of a data frame.

Value

A named list of numeric vectors, where each element corresponds to a stratum and contains the event times of the selected type.

Lifecycle

See Also

polyreg() for log-odds product modeling of CIFs; cifcurve() for KM/AJ estimators; cifplot() for display of a CIF; cifpanel() for display of multiple CIFs; ggsurvfit::ggsurvfit, patchwork::patchwork and modelsummary::modelsummary for display helpers.

Examples

data(diabetes.complications)
output <- extract_time_to_event(Event(t,epsilon) ~ fruitq,
                                data = diabetes.complications,
                                which.event = "event2")
cifplot(Event(t,epsilon) ~ fruitq,
        data = diabetes.complications,
        outcome.type="competing-risk",
        add.conf=FALSE,
        add.risktable=FALSE,
        add.censor.mark=FALSE,
        add.competing.risk.mark=TRUE,
        competing.risk.time=output,
        label.y="CIF of diabetic retinopathy",
        label.x="Years from registration")

Choose a safe finite-difference step so w +/- h*delta stays >= 0

Description

Choose a safe finite-difference step so w +/- h*delta stays >= 0

Usage

fd_step_safe(w, delta, rel = 1e-06, max_tries = 20L)

Extract Mparts from a weightit object (robustly)

Description

Extract Mparts from a weightit object (robustly)

Usage

get_weightit_mparts(weightit)

Fit coherent regression models of CIFs using polytomous log odds products

Description

polyreg() fits regression models of CIFs, targeting familiar effect measures (risk ratios, odds ratios and subdistribution hazard ratios). Modeling the nuisance structure using polytomous log odds products ensures that the sum of cause-specific CIFs does not exceed one, and enables coherent modelling of the multiplicative effects.

This function follows a familiar formula–data workflow: the outcome and covariates other than the exposure are specified through a formula in nuisance.model (with Event() or survival::Surv() on the LHS), and the exposure of interest is given by a separate variable name in exposure. The fitted object contains tidy summaries of exposure effects (point estimates, SEs, CIs, and p-values) and can be summarised with summary.polyreg() or formatted with external tools such as modelsummary::modelsummary().

Usage

polyreg(
  nuisance.model,
  exposure,
  strata = NULL,
  data,
  subset.condition = NULL,
  na.action = na.omit,
  code.event1 = 1,
  code.event2 = 2,
  code.censoring = 0,
  code.exposure.ref = 0,
  effect.measure1 = "RR",
  effect.measure2 = "RR",
  time.point = NULL,
  outcome.type = "competing-risk",
  conf.int = 0.95,
  report.nuisance.parameter = FALSE,
  report.optim.convergence = FALSE,
  report.sandwich.conf = TRUE,
  report.boot.conf = NULL,
  boot.bca = FALSE,
  boot.multiplier = "rademacher",
  boot.replications = 200,
  boot.seed = 46,
  nleqslv.method = "Newton",
  optim.parameter1 = 1e-06,
  optim.parameter2 = 1e-06,
  optim.parameter3 = 100,
  optim.parameter4 = 50,
  optim.parameter5 = 50,
  optim.parameter6 = 50,
  optim.parameter7 = 1e-10,
  optim.parameter8 = 1e-06,
  optim.parameter9 = 0.001,
  optim.parameter10 = 40,
  optim.parameter11 = 1e-06,
  optim.parameter12 = 2,
  optim.parameter13 = 0.5,
  data.initial.values = NULL,
  normalize.covariate = TRUE,
  terminate.time.point = TRUE,
  prob.bound = 1e-07
)

Arguments

Details

Overview

polyreg() implements log odds product modeling for CIFs at user-specified time points, focusing on multiplicative effects of a categorical exposure, or constant effects over time like Cox regression and Fine-Gray models. It estimates multiplicative effects such as risk ratios, odds ratios, or subdistribution hazard ratios, while ensuring that the probabilities across competing events sum to one. This is achieved through reparameterization using polytomous log odds products, which fits so-called effect-measure models and nuisance models on multiple competing events simultaneously (Tanaka and Scheike 2026). Additionally, polyreg() supports direct binomial regression for survival outcomes and the Richardson model for binomial outcomes, both of which use log odds products. With more than two exposure levels, the single-outcome model uses the generalized log odds product, the sum of log odds across exposure levels.

Key arguments

• nuisance.model: a formula with Event() or survivai::Surv() describing the outcome and nuisance covariates, excluding the exposure of interest.

• exposure: name of the categorical exposure variable

• effect.measure1 and effect.measure2: the effect measures for event1 and event2 ("RR", "OR" or "SHR").

• outcome.type: type of the outcome variable ("competing-risk", "survival", "binomial", "proportional-survival" or "proportional-competing-risk").

• time.point: time point(s) at which the exposure effect is evaluated. Required for "competing-risk" and "survival" outcomes.

• strata: name of the stratification variable used for IPCW adjustment for dependent censoring.

Outcome type and event status coding

The outcome.type argument must be set to:

• Effects on cumulative incidence probabilities at a specific time: "competing-risk".

• Effects on a risk at a specific time: "survival".

• Common effects on cumulative incidence probabilities over time: "proportional-competing-risk".

• Common effects on a risk over time: "proportional-survival".

• Effects on a risk of a binomial outcome: "binomial".

Effect measures for categorical exposure

Choose the effect scale for event 1 and (optionally) event 2:

• RR: risk ratio at time.point or common over time.

• OR: odds ratio at time.point or common over time.

• SHR: subdistribution hazard ratio or common over time.

Inference and intervals (advanced)

Optimization & solver controls (advanced)

polyreg() solves estimating equations with optional inner routines.

Data handling and stability

Downstream use

polyreg() returns an object of class "polyreg" that contains regression coefficients (coef), variance-covariance matrix (vcov) and a list of event-wise tidy and glance tables (summary). Users should typically access results via the S3 methods:

• coef() — extract regression coefficients.

• vcov() — extract the variance–covariance matrix (sandwich or bootstrap, depending on outcome.type and the ⁠report.*⁠ arguments).

• nobs() — number of observations used in the fit.

• summary() — print an event-wise, modelsummary-like table of estimates, CIs and p-values, and return the underlying list of tidy/glance tables invisibly.

For backward compatibility, components named coefficient and cov may also be present and mirror coef and vcov, respectively. The summary component can be passed to external functions such as modelsummary() for further formatting, if desired.

Reproducibility and conventions

• If convergence warnings appear, relax/tighten tolerances or cap the parameter bound (optim.parameter1–3) and inspect the output with report.optim.convergence = TRUE.

• If necessary, modify other optim.parameter, provide user-specified initial values, or reduce the number of nuisance parameters (e.g., provide a small set of time points contributing to estimation when using "proportional-survival" or "proportional-competing-risk").

• Set boot.seed for reproducible bootstrap results.

• Match CDISC ADaM conventions via code.event1 = 0, code.censoring = 1 (and, if applicable, code.event2 for competing events).

Value

A list of class "polyreg" containing the fitted exposure effects and supporting results. Key components and methods include:

• coef: regression coefficients on the chosen effect-measure scale

• vcov: variance–covariance matrix of the regression coefficients

• diagnostic.statistics: a data frame with inverse probability weights, influence function contributions, and predicted potential outcomes

• summary: event-wise tidy/glance summaries used by summary.polyreg() or modelsummary::modelsummary()

• additional elements storing convergence information and internal tuning parameters.

Standard S3 methods are available: coef.polyreg(), vcov.polyreg(), nobs.polyreg(), and summary.polyreg().

Lifecycle

References

Tanaka S, Scheike TH. Estimating effects on cumulative incidence probabilities by direct polytomous regression and polytomous log-odds product. Biometrics. 2026, 82(3), ujag160. doi:10.1093/biomtc/ujag160

See Also

cifcurve() for KM/AJ estimators; cifplot() for display of a CIF; cifpanel() for display of multiple CIFs; ggsurvfit::ggsurvfit, patchwork::patchwork and modelsummary::modelsummary for display helpers.

Examples

data(diabetes.complications)
output <- polyreg(
  nuisance.model = Event(t, epsilon) ~ +1,
  exposure = "fruitq1",
  data = diabetes.complications,
  effect.measure1 = "RR",
  effect.measure2 = "RR",
  time.point = 8,
  outcome.type = "competing-risk"
)

coef(output)
vcov(output)
nobs(output)
summary(output)

if (requireNamespace("gtsummary", quietly = TRUE) &&
    utils::packageVersion("gtsummary") >= "2.6.1") {
  gtsummary::tbl_regression(output, exponentiate = TRUE)
  gtsummary::tbl_regression(output, exponentiate = TRUE, event = "event2")
}

Methods for polyreg objects

Description

S3 methods to extract coefficients, variance-covariance matrix, sample size, formatted summaries, and tidy/glance/augment from objects returned by polyreg().

Usage

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

## S3 method for class 'polyreg'
vcov(
  object,
  type = c("default", "sandwich", "bootstrap"),
  variance = NULL,
  ...
)

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

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

## S3 method for class 'summary.polyreg'
print(x, digits = 3, ...)

effect_label.polyreg(
  x,
  event = c("event1", "event2"),
  add.time.point = TRUE,
  add.outcome = TRUE,
  add.exposure.levels = TRUE,
  add.conf = TRUE,
  add.p = TRUE,
  value.time = NULL,
  unit.time = NULL,
  digits = 2,
  p_digits = 2,
  p_cut = 0.05,
  ...
)

## S3 method for class 'polyreg'
tidy(
  x,
  event = c("event1", "event2", "both"),
  conf.int = x$conf.int,
  variance = NULL,
  exponentiate = FALSE,
  conf.level = NULL,
  ...,
  vcov = NULL,
  coef_rename = FALSE
)

## S3 method for class 'polyreg'
glance(x, event = c("event1", "event2"), ..., gof_map = NULL)

## S3 method for class 'polyreg'
augment(x, ...)

## S3 method for class 'polyreg'
tbl_regression(x, ..., event = c("event1", "event2"))

Arguments

Value

• coef.polyreg() returns a numeric vector of regression coefficients.

• vcov.polyreg() returns a variance-covariance matrix.

• nobs.polyreg() returns the number of observations.

• summary.polyreg() returns a list of tidy and glance summaries by event.

• print.summary.polyreg() is called for its side effect of printing a formatted, modelsummary-like table to the console and returns x invisibly.

• tidy.polyreg() returns a data frame of tidy coefficients by event.

• glance.polyreg() returns a data frame of model-level summaries by event.

• augment.polyreg() returns an augmented data frame containing diagnostics, weights, and predicted CIFs.

See Also

polyreg() for log odds product modeling of CIFs

Evaluate a censoring Kaplan-Meier distribution

Description

Evaluate a censoring Kaplan-Meier distribution

Usage

predict_censoring_km(
  object,
  time,
  strata = rep.int("pooled", length(time)),
  side = c("left", "right"),
  use = c("used", "raw")
)

Evaluate a working Aalen-Johansen distribution

Description

Evaluate a working Aalen-Johansen distribution

Usage

predict_working_aj(
  object,
  time,
  exposure,
  strata = rep.int("pooled", length(time)),
  side = c("left", "right"),
  use = c("used", "raw")
)

Data from a prostate cancer trial in Byer & Green (1980)

Description

Anonymized data from a randomized clinical trial of prostate cancer published in Byer & Green (1980).

Usage

data(prostate)

Format

A data frame with 502 observations and 18 variables, including:

dtime

Follow-up time in months.

status

Event status ("alive", "dead - prostatic ca", "dead - other ca", "dead - heart or vascular", "dead - cerebrovascular").

rx

Treatment assignment to diethylstilbestrol (DES) or a placebo.

age

Age at baseline (years).

wt

Weight in pounds.

pf

Performance status.

hx

History of cardiovascular disease.

sbp

Systolic blood pressure.

dbp

Diastolic blood pressure.

ekg

Electrocardiogram category.

hg

Hemoglobin level.

sz

Size of the primary tumor.

sg

Stage/grade of disease.

ap

Serum acid phosphatase.

bm

Bone metastases indicator.

stage

Clinical stage.

sdate

Start date.

patno

Patient number.

Details

The dataset records follow-up for cause of death together with treatment assignment and baseline characteristics. It is used in the package documentation to illustrate stratified cumulative incidence analyses.

Source

Byer, D. P. & Green, S. B. (1980), 'Prognostic variables for survival in a randomized comparison of treatments for prostatic cancer', Bulletin du Cancer 67, 477-488

Examples

data(prostate)
head(prostate)