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().
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 = 1e-06,
optim.parameter10 = 40,
optim.parameter11 = 0.025,
optim.parameter12 = 2,
optim.parameter13 = 0.5,
data.initial.values = NULL,
normalize.covariate = TRUE,
terminate.time.point = TRUE,
prob.bound = 1e-07
)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().
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).
A character string giving the name of the categorical exposure
variable in data.
Optional character string with the name of the stratification
variable used to adjust for dependent censoring (default NULL).
A data frame containing the outcome, exposure and nuisance
covariates referenced by nuisance.model.
Optional character string giving a logical condition to subset
data (default NULL).
A function specifying the action to take on missing values (default na.omit).
Integer code of the event of interest (default 1).
Integer code of the competing event (default 2).
Integer code of censoring (default 0).
Integer code identifying the reference exposure
category (default 0).
Character string specifying the effect measure for the
primary event. Supported values are "RR", "OR" and
"SHR".
Character string specifying the effect measure for the
competing event. Supported values are "RR", "OR" and
"SHR".
Numeric time point at which the exposure effect is evaluated for
time-point models. Required for "competing-risk" and "survival"
outcomes.
Character string selecting the outcome type. Valid values are
"competing-risk", "survival", "binomial", "proportional-survival",
and "proportional-competing-risk". The default is "competing-risk".
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.
Numeric two-sided level of CIs (default 0.95).
Logical; if TRUE, the returned object
includes estimates of the nuisance model parameters (default FALSE).
Logical; if TRUE, optimization
convergence summaries are returned (default FALSE).
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.
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").
Logical indicating the bootstrap CI method.
Use TRUE for bias-corrected and accelerated intervals or FALSE
for the normal approximation (default FALSE).
Character string specifying the wild bootstrap weight distribution.
One of "rademacher", "mammen" or "gaussian" (default "rademacher").
Integer giving the number of bootstrap replications
(default 200).
Numeric seed used for resampling of bootstrap.
Character string specifying the solver used in
nleqslv(). Available choices are "Broyden" and "Newton".
Numeric tolerance for convergence of the outer loop
(default 1e-6).
Numeric tolerance for convergence of the inner loop
(default 1e-6).
Numeric constraint on the absolute value of
parameters (default 100).
Integer maximum number of outer loop iterations
(default 50).
Integer maximum number of nleqslv
iterations per outer iteration (default 50).
Integer maximum number of iterations for the
Levenberg-Marquardt routine (default 50).
Numeric convergence tolerance for the
Levenberg-Marquardt routine (default 1e-10).
Numeric tolerance for updating the Hessian in the
Levenberg-Marquardt routine (default 1e-6).
Numeric starting value for the Levenberg-Marquardt
damping parameter lambda (default 1e-6).
Numeric upper bound for lambda in the
Levenberg-Marquardt routine (default 40).
Numeric lower bound for lambda in the
Levenberg-Marquardt routine (default 0.025).
Numeric multiplicative increment applied to lambda
when the Levenberg-Marquardt step is successful (default 2).
Numeric multiplicative decrement applied to lambda
when the Levenberg-Marquardt step is unsuccessful (default 0.5).
Optional data frame providing starting values for
the optimization (default NULL).
Logical indicating whether covariates should
be centered and scaled prior to optimization (default TRUE).
Logical indicating whether time points
that contribute estimation are terminated by min of max follow-up times
of each exposure level (default TRUE).
Numeric lower bound used to internally truncate probabilities away
from 0 and 1 (default 1e-5).
lifecycle::badge("experimental")
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. Additionally, polyreg() supports direct binomial regression
for survival outcomes and the Richardson model for binomial outcomes,
both of which use log odds products.
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.
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 |
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.
| 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 |
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 |
| 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 |
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.
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).
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.
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)
Run the code above in your browser using DataLab