polyreg: Fit coherent regression models of CIFs using polytomous log...

In cifmodeling: Visualization and Polytomous Modeling of Survival and Competing Risks

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

nuisance.model	A formula describing the outcome and nuisance covariates, excluding the exposure of interest. The LHS must be Event(time, status) or survival::Surv(time, status). Surv() responses must be right-censored.
exposure	A character string giving the name of the categorical exposure variable in data.
strata	Optional character string with the name of the stratification variable used to adjust for dependent censoring (default NULL).
data	A data frame containing the outcome, exposure and nuisance covariates referenced by nuisance.model.
subset.condition	Optional character string giving a logical condition to subset data (default NULL).
na.action	A function specifying the action to take on missing values (default na.omit).
code.event1	Integer code of the event of interest (default 1).
code.event2	Integer code of the competing event (default 2). Used for competing-risk models and when inferring the outcome type.
code.censoring	Integer code of censoring (default 0). The event codes used by the model must be distinct non-negative integers.
code.exposure.ref	Integer code identifying the reference exposure category (default 0).
effect.measure1	Character string specifying the effect measure for the primary event. Supported values are "RR", "OR" and "SHR".
effect.measure2	Character string specifying the effect measure for the competing event. Supported values are "RR", "OR" and "SHR".
time.point	Numeric time point at which the exposure effect is evaluated for time-point models. Required for "competing-risk" and "survival" outcomes.
outcome.type	Character string selecting the outcome type. Valid values are "competing-risk", "survival", "binomial", "proportional-survival", and "proportional-competing-risk". The default is "competing-risk". Fitting either competing-risk model requires observations of both event types. For "competing-risk", both types must be observed by time.point. If explicitly set to NULL, polyreg() attempts to infer the outcome type from the data: if the event variable has more than two distinct levels, "competing-risk" is assumed; otherwise, "survival" is assumed. Abbreviations such as "S" or "C" are accepted; mixed or ambiguous inputs trigger automatic detection from the event coding in data.
conf.int	Numeric two-sided level of CIs (default 0.95).
report.nuisance.parameter	Logical; if TRUE, the returned object includes estimates of the nuisance model parameters (default FALSE).
report.optim.convergence	Logical; if TRUE, optimization convergence summaries are returned (default FALSE).
report.sandwich.conf	Logical or NULL. When TRUE, CIs based on sandwich variance are computed. When FALSE, they are omitted (default TRUE). This CI is default for time-point models ("outcome.type=competing-risk", "survival" or "binomial") and is not available otherwise. For time-point survival and competing-risk models, the variance accounts for estimation of the censoring distribution within each supplied stratum.
report.boot.conf	Logical or NULL. When TRUE, bootstrap CIs are computed. When FALSE, they are omitted. If NULL, the function chooses based on outcome.type (default NULL). This CI is default for proportional models (outcome.type="proportional-competing-risk" or "proportional-survival").
boot.bca	Logical indicating the bootstrap CI method. Use TRUE for bias-corrected and accelerated intervals or FALSE for the normal approximation (default FALSE). For proportional models, normal intervals are centered on the mean bootstrap coefficient and use the requested conf.int. BCa limits are computed separately for each reported coefficient.
boot.multiplier	Character string specifying the wild bootstrap weight distribution. One of "rademacher", "mammen" or "gaussian" (default "rademacher").
boot.replications	Integer giving the number of bootstrap replications (default 200).
boot.seed	Numeric seed used for resampling of bootstrap.
nleqslv.method	Character string specifying the solver used in nleqslv(). Available choices are "Broyden" and "Newton".
optim.parameter1	Numeric tolerance for convergence of the outer loop (default 1e-6).
optim.parameter2	Numeric tolerance for convergence of the inner loop (default 1e-6).
optim.parameter3	Numeric constraint on the absolute value of parameters (default 100).
optim.parameter4	Integer maximum number of outer loop iterations (default 50).
optim.parameter5	Integer maximum number of nleqslv iterations per outer iteration (default 50).
optim.parameter6	Integer maximum number of iterations for the Levenberg-Marquardt routine (default 50).
optim.parameter7	Numeric convergence tolerance for the Levenberg-Marquardt routine (default 1e-10).
optim.parameter8	Numeric tolerance for updating the Hessian in the Levenberg-Marquardt routine (default 1e-6).
optim.parameter9	Numeric starting value for the Levenberg-Marquardt damping parameter lambda (default 1e-3).
optim.parameter10	Numeric upper bound for lambda in the Levenberg-Marquardt routine (default 40).
optim.parameter11	Numeric lower bound for lambda in the Levenberg-Marquardt routine (default 1e-6).
optim.parameter12	Numeric multiplicative decrement applied to lambda when the Levenberg-Marquardt step is unsuccessful (default 2).
optim.parameter13	Numeric multiplicative increment applied to lambda when the Levenberg-Marquardt step is successful (default 0.5).
data.initial.values	Optional data frame providing starting values for the optimization (default NULL).
normalize.covariate	Logical indicating whether covariates should be centered and scaled prior to optimization (default TRUE).
terminate.time.point	Logical indicating whether time points that contribute estimation are terminated by min of max follow-up times of each exposure level (default TRUE).
prob.bound	Numeric lower bound used to internally truncate probabilities away from 0 and 1 (default 1e-7).

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

Setting	Codes	Meaning
competing-risk	code.event1, code.event2, code.censoring	event of interest / competing event / censoring
competing-risk (default)	code.event1 = 1, code.event2 = 2, code.censoring = 0	event of interest / competing event / censoring
survival	code.event1, code.censoring	event / censoring
survival (default)	code.event1 = 1, code.censoring = 0	event / censoring
survival (ADaM-ADTTE)	code.event1 = 0, code.censoring = 1	set to match ADaM convention
proportional-survival	code.event1, code.censoring	event / censoring
proportional-survival (default)	code.event1 = 1, code.censoring = 0	event / censoring
proportional-survival (ADaM)	code.event1 = 0, code.censoring = 1	set to match ADaM convention
proportional-competing-risk	code.event1, code.event2, code.censoring	event of interest / competing event / censoring
proportional-competing-risk (default)	code.event1 = 1, code.event2 = 2, code.censoring = 0	event of interest / competing event / censoring

Effect measures for categorical exposure

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

Argument	Applies to	Choices	Default
effect.measure1	event of interest	"RR", "OR", "SHR"	"RR"
effect.measure2	competing event	"RR", "OR", "SHR"	"RR"

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

Argument	Meaning	Default
conf.int	Wald-type CI level	0.95
report.sandwich.conf	Sandwich variance CIs	TRUE
report.boot.conf	Bootstrap CIs (used by "proportional-*" types)	NULL
boot.bca	Use BCa intervals (else normal approximation)	FALSE
boot.multiplier	Method for wild bootstrap	"rademacher"
boot.replications	Bootstrap replications	200
boot.seed	Seed for resampling	46

Optimization & solver controls (advanced)

polyreg() solves estimating equations with optional inner routines.

Argument	Role	Default
nleqslv.method	Root solver	"Newton"
optim.parameter1, optim.parameter2	Outer / inner convergence tolerances	1e-6, 1e-6
optim.parameter3	Parameter absolute bound	100
optim.parameter4	Max outer iterations	50
optim.parameter5	Max nleqslv iterations per outer	50
optim.parameter6:13	Levenberg–Marquardt controls (iterations, tolerances, lambda)	see defaults

Data handling and stability

Argument	Meaning	Default
subset.condition	Expression (as character) to subset data	NULL
na.action	NA handling function	stats::na.omit
normalize.covariate	Center/scale nuisance covariates	TRUE
terminate.time.point	Truncate support by exposure-wise follow-up maxima	TRUE
prob.bound	Truncate probabilities away from 0/1 (numerical guard)	1e-5
data.initial.values	Optional starting values data frame	NULL

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

\lifecycle

experimental

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. \Sexpr[results=rd]{tools:::Rd_expr_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")
}

cifmodeling documentation built on Nov. 11, 2026, 9:06 a.m.