Estimating Adjusted Cumulative Incidence with Cause-Specific Cox Models

From Cause-Specific Hazards to Absolute Risks

A cause-specific Cox model describes the instantaneous rate of one event type among individuals who are still event-free. That rate is useful for describing covariate associations, but a hazard ratio alone is not an absolute probability of experiencing the event by a given time. To estimate that probability, we must combine the cause-specific hazards for all competing event types: every cause contributes to the probability of remaining event-free, and therefore affects the cumulative incidence of each cause.

The goal here is a covariate-adjusted cumulative incidence curve: the estimated probability of experiencing a specified cause by each time, for a person or profile with specified covariate values. This is a model-based prediction, conditional on the fitted models and their assumptions. It is not obtained by interpreting one Cox model’s hazard ratio as a risk ratio.

For continuous event-time distributions, the link between these quantities is \(F_k(t\mid z)=\int_0^t S(u\mid z)h_k(u\mid z)\,du\), where \(S(t\mid z)=\exp\{-\sum_{j=1}^K H_j(t\mid z)\}\), and \(H_j\) is the cumulative cause-specific hazard. Thus, the risk from cause \(k\) depends both on its own hazard and on the hazards of all competing causes.

Assumptions and Calculation

A cause-specific Cox fit treats an event of another cause as ending that individual’s time in the event-free risk set. This is the appropriate risk set for estimating the cause-specific hazard; it does not mean that the competing event is ordinary loss to follow-up. Loss to follow-up must still be handled under an independent-censoring assumption, typically conditional on the covariates in the model. The event cause must be recorded accurately, and the models must adequately describe the cause-specific hazards for the covariate profiles and times being predicted.

Each cause-specific model also assumes proportional hazards for its included covariates: the relative hazard associated with a fixed covariate contrast does not vary over time. If this assumption is implausible, a single constant coefficient may not support reliable predictions across the full time range. Predictions for a profile with time-varying covariates require specifying that covariate history; the fixed-profile calculation below does not silently accommodate changing covariates.

In practice, the baseline cumulative hazards are estimated from each fitted Cox model, often on a grid of event times. Multiply each cause’s baseline hazard increments by that cause’s covariate multiplier, then combine the causes. If the cumulative hazard increments over an interval are treated as arising from constant hazards within that interval, the interval calculation is exact under that piecewise-constant assumption. It is also a useful numerical approximation when the increments come from a fine event-time grid.

When \(q=0\), no cause has positive hazard in the interval, so every \(\Delta F_k=0\) and survival is unchanged. When \(q>0\), the factor \(\{1-\exp(-q)\}\) is the total probability of an event during the interval conditional on being event-free at its start; \(\Delta H_k/q\) allocates that probability to cause \(k\) in proportion to its hazard. The updated cumulative incidence is the previous cumulative incidence plus \(\Delta F_k\). This calculation preserves the probability accounting: the interval’s cause-specific risks sum to the decrease in event-free survival.

An equivalent continuous-time expression is to integrate \(S(u\mid z)h_k(u\mid z)\) over time. On a fine grid, numerical integration approximates this expression. Whichever computational convention is used, the key requirement is to update survival using the hazards of all causes, not just the cause whose risk is being reported.

Fitting Models and Predicting for a Profile

With a two-cause outcome coded as 0 for censoring, 1 for cause 1, and 2 for cause 2, a common R workflow is to fit one Cox model per cause. For example, using the survival package, the model formulas are written as `Surv(time, status == 1) ~ age + treatment` and `Surv(time, status == 2) ~ age + treatment`. In each fit, the other cause is coded as not having the event of interest, so that it leaves that cause’s event-free risk set. This is not a claim that the person remains under observation after the competing event.

After fitting, obtain each model’s estimated baseline cumulative hazard and take differences between consecutive rows to get baseline increments. In R, `basehaz(fit, centered = FALSE)` reports the cumulative hazard at the zero-valued covariate reference for the fitted linear predictor. The zero reference may not describe a realistic person; it is a mathematical reference, and the profile-specific multiplier translates it to the requested profile. Align the causes’ increments on a common time grid, using zero increments for a cause at a time where its baseline hazard does not jump. Then apply the update equations to the chosen profile and accumulate the cause-specific risks.

A practical check is that the fitted coefficients used in the multipliers correspond to the same covariate coding as the profile values. A treatment indicator coded 0/1 must be entered as 0 or 1, not as a verbal label unless the software handles the factor coding consistently. Continuous covariates must use their fitted units. Centering a continuous variable, such as age in decades from a meaningful reference age, can make the baseline reference easier to interpret without changing fitted relative hazards.

Selected Times, Diagnostics, and Limits

To report risk at a selected time \(t^\ast\), include all baseline hazard increments up to \(t^\ast\) and sum the corresponding cause-specific increments. If \(t^\ast\) falls between grid points, a stated interpolation or piecewise-constant hazard convention is needed; simply taking the nearest time without saying so can misstate the requested absolute risk. Predictions beyond the period with adequate follow-up rely increasingly on the model’s extrapolation and should be treated cautiously.

The fitted curve should be nondecreasing for each cause, and the sum of all cause-specific cumulative incidences plus event-free survival should be approximately one at every time. A decrease in a cumulative incidence, a negative survival probability, or a sum materially different from one usually signals a programming error, inconsistent time grids, or an inappropriate update formula. Minor discrepancies can result from rounding.

Other warnings concern the model rather than the arithmetic. Check proportional-hazards diagnostics for each cause-specific model, such as the Schoenfeld-residual checks used for Cox models. Inspect whether the profile’s covariate values are represented in the data; predictions far outside observed covariate ranges are extrapolations. Consider whether important predictors have been omitted and whether loss to follow-up could depend on unmodeled factors. A correct calculation cannot repair misspecified cause-specific hazards or informative censoring.

Cause-specific Cox models are often useful when the scientific question concerns instantaneous event rates among those still event-free, and their fitted hazards can be combined to produce absolute risks. Fine-Gray models instead model a subdistribution hazard directly and can relate covariates to the cumulative incidence through a different risk-set construction. Neither model’s hazard ratio is itself an absolute risk. The approach here requires fitting and combining all relevant cause-specific hazards; fitting only the model for the target cause is insufficient for its cumulative incidence.

Worked Calculation with an Estimated Hazard Grid

These estimates are point predictions. Uncertainty in regression coefficients and baseline hazards means the estimated curve itself is uncertain. Pointwise confidence intervals require an appropriate variance procedure, such as a method based on the fitted model’s covariance information or resampling. The values calculated here should not be read as certainty about an individual’s outcome, nor as evidence that changing a covariate would causally change that person’s risk. Interpretation depends on the study design, covariate adjustment, and model assumptions.