Fitting Cause-Specific Cox Models with coxph

What a Cause-Specific Cox Fit Estimates

When an individual can experience one of several event types, a cause-specific Cox model describes the instantaneous rate of one selected event among people who are still free of every event. Fitting separate models is a practical way to estimate how covariates relate to each cause-specific hazard. In R, each fit uses the ordinary Cox regression function survival::coxph, with a different event indicator for each cause.

The data should have one row per individual, a follow-up time, a status code identifying censoring or event type, and the covariates to include. Here, assume status is coded 0 for right censoring and 1, 2, and so on for event causes. The example code below uses a hypothetical dataset with those columns.

Definition: For cause $k$, define $\delta_{ik}=1$ if individual $i$ experiences cause $k$, and $\delta_{ik}=0$ otherwise. In the cause-$k$ Cox fit, an event from another cause is recorded as a non-target event: it removes the individual from the risk set at that event time, just as right censoring removes an individual from the risk set at the censoring time. The cause-specific hazard is the instantaneous rate of cause $k$ among individuals still free of all event types immediately before time $t$.

This coding defines the cause-specific hazard; it does not mean that competing events are unobserved or that they can happen after they have occurred. Once a competing event happens, that individual cannot subsequently experience the target event as a first event in this analysis. Right censoring, by contrast, means the event type and time after the recorded follow-up are unknown.

Model and Assumptions

For a selected cause $k$, the cause-specific proportional hazards model is

\lambda_k(t \mid \mathbf{x}) = \lambda_{0k}(t)\exp(\boldsymbol{\beta}_k^\mathsf{T}\mathbf{x}).

Here $\lambda_k(t \mid \mathbf{x})$ is the cause-$k$ hazard at time $t$ for covariate vector $\mathbf{x}$, conditional on being free of every event just before $t$. The unspecified baseline hazard $\lambda_{0k}(t)$ corresponds to the covariate reference $\mathbf{x}=\mathbf{0}$; $\boldsymbol{\beta}_k$ is the cause-specific coefficient vector. A one-unit increase in a covariate $x_j$, holding the other modeled covariates fixed, multiplies this hazard by $\exp(\beta_{kj})$, assuming proportional hazards.

The main assumptions to consider before interpreting a fit are:

Correct event and time coding. Each individual has the correct follow-up time and one correctly classified status value. Follow-up is measured from a common, well-defined time origin.
Noninformative right censoring, conditional on modeled covariates. Among individuals with the same modeled covariates, censoring should not selectively depend on their unobserved future target-event experience in a way that distorts the risk sets.
Proportional hazards for the selected cause. The covariate hazard ratios are constant over follow-up time under the specified model.
Appropriate model form. The included covariates, their functional forms, and any interactions should adequately represent the intended conditional association. Independent individuals are also assumed for the usual model-based standard errors; clustered or otherwise dependent observations require suitable variance methods.

The competing events themselves are not treated as ordinary loss to follow-up for defining the target cause-specific hazard: their occurrence is observed, and they end the individual's time at risk for a first event of the target cause. Interpretations concern the hazard among those still event-free, not the probability of eventually experiencing that cause.

Why the Separate Cause Fits Work

At a time when a cause-$k$ event occurs, the Cox partial likelihood compares the covariates of the individual with that event to the covariates of all individuals still at risk immediately beforehand. For a single event at time $t$, let $R(t)$ denote that risk set and let $i$ be the individual with the cause-$k$ event. The contribution, ignoring ties, is

\frac{\exp(\boldsymbol{\beta}_k^\mathsf{T}\mathbf{x}_i)} {\sum_{j \in R(t)} \exp(\boldsymbol{\beta}_k^\mathsf{T}\mathbf{x}_j)}.

This contribution follows by conditioning on one cause-$k$ event occurring in that risk set at that time. Each individual's cause-$k$ hazard has the common baseline factor $\lambda_{0k}(t)$; dividing the event individual's hazard by the sum of cause-$k$ hazards in the risk set cancels that factor. Individuals with a competing event remain in the risk set up to their event time, then leave it. The product of the contributions across target-cause event times gives the partial likelihood used to estimate $\boldsymbol{\beta}_k$. Tied event times require a tie-handling method; `coxph` defaults to Efron's method.

Worked Example: Constructing Cause-Specific Event Indicators

Suppose a hypothetical dataset records follow-up in months and uses status 0 for censoring, 1 for relapse, and 2 for death before relapse. For a relapse fit, status 1 is the event of interest; statuses 0 and 2 are non-events for that fit. For a death fit, status 2 is the event of interest; statuses 0 and 1 are non-events.

Individual	Follow-up (months)	Status	Relapse indicator	Death indicator
A	4	1	1	0
B	7	2	0	1
C	9	0	0	0

For individual B, the relapse analysis uses a time of 7 months and a relapse indicator of 0. B contributes to the relapse risk sets before month 7, then leaves at month 7 because death prevents a subsequent first relapse from being observed. Individual C is right-censored at month 9 and contributes to the relapse risk sets only up to that time. The death analysis applies the same logic with the target event changed to death.

In R, the indicators and fits can be created as follows; assume `dat` also contains numeric covariates `age` and `treatment`:

`library(survival)`
`dat$relapse <- as.integer(dat$status == 1)`
`dat$death <- as.integer(dat$status == 2)`
`fit_relapse <- coxph(Surv(time, relapse) ~ age + treatment, data = dat)`
`fit_death <- coxph(Surv(time, death) ~ age + treatment, data = dat)`

The Boolean comparisons convert each status value to 1 only for its target cause. Check that the status coding and the chosen `time` column really have the meanings assumed here. In particular, using `status != 0` as the relapse indicator would incorrectly count deaths as relapse events.

Fitting and Interpreting the Coefficients

Each cause-specific fit has its own coefficients, baseline hazard, and event indicator. Use the same covariates in each fit when a direct comparison of cause-specific associations is intended, but do not assume that coefficient estimates must be equal across causes. The example below reports hypothetical output from the relapse model; the numbers are illustrative, not the result of fitting the three-row table above.

Worked Example: Converting a Coefficient to a Hazard Ratio

Suppose the fitted relapse model reports an age coefficient of 0.0400 per year, with model-based standard error 0.0150. The estimated cause-specific hazard ratio for a one-year increase in age is

\widehat{HR} = \exp(0.0400) = 1.041.

Under the proportional hazards model, this is an estimated 4.08% higher relapse hazard per additional year of age, conditional on the other included covariates and among individuals still event-free. It is an association with the cause-specific hazard, not a 4.08 percentage-point increase in relapse probability.

A Wald 95% confidence interval is calculated on the coefficient scale and then exponentiated:

0.0400 \mathbin{\pm} 1.96(0.0150) = 0.0400 \mathbin{\pm} 0.0294 = (0.0106,\;0.0694), $$ $$ \left(\exp(0.0106),\;\exp(0.0694)\right) = (1.011,\;1.072).

Thus the hypothetical interval is consistent with a relapse hazard between about 1.1% and 7.2% higher per year of age, conditional on the modeled covariates. The calculation uses a normal approximation for the coefficient. To request this form of interval from a fitted model in R, use `exp(cbind(HR = coef(fit_relapse), confint.default(fit_relapse)))`; `confint.default` constructs the Wald interval using the estimated coefficient covariance matrix.

For the same hypothetical estimate, the Wald statistic is $z=0.0400/0.0150=2.67$, giving a two-sided normal-approximation p-value of about 0.0077. This is evidence against a zero age coefficient under the model and its assumptions; it is not the probability that the null hypothesis is true.

Extracting Cumulative Baseline Hazards

A Cox fit does not require a parametric formula for the baseline hazard, but it is possible to estimate and extract a cumulative baseline hazard. In R, `basehaz` returns a data frame with event times and a cumulative hazard column named `hazard`. Its `centered` argument matters: with `centered = FALSE`, the returned curve corresponds to covariate values zero. With the default centering, the reference is instead at the fitted model's centering values. A zero-valued reference is meaningful only if zero is a sensible value for every covariate and factor coding has a meaningful reference interpretation.

For a model with numeric age and a 0/1 treatment indicator, a basic extraction is `bh <- basehaz(fit_relapse, centered = FALSE)`. Inspect `head(bh)` and the model's contrasts before interpreting the reference curve. If age zero is not scientifically meaningful, a more interpretable reference can be created by centering age before fitting (for example, age minus 50 years), then extracting the baseline hazard at centered covariates zero. For factors, verify how the model matrix encodes the reference level; setting all dummy variables to zero typically corresponds to the reference category.

Worked Example: Using an Extracted Baseline Hazard

Suppose a hypothetical fitted model uses age centered at 50 years and a treatment indicator coded 0 for control and 1 for active treatment. Suppose `basehaz(fit, centered = FALSE)` returns a cumulative baseline hazard of 0.240 at 12 months. This is the estimated cumulative cause-specific hazard through 12 months for the reference covariates: age 50 and control treatment.

Suppose the fitted treatment coefficient is $-0.4308$, so the treatment hazard ratio is $\exp(-0.4308)=0.6500$. At age 50, the estimated 12-month cumulative cause-specific hazard for active treatment is scaled by this relative hazard:

\widehat{H}_k(12 \mid \text{age}=50,\text{ active}) = 0.240\exp(-0.4308) = 0.240(0.6500) = 0.156.

This calculation relies on proportional hazards and on the stated coding and reference values. The value 0.156 is a cumulative cause-specific hazard, not a 15.6% probability of relapse. Transforming it to a cause-specific survival quantity as $\exp(-0.156)=0.8556$ gives the fitted survival function associated with that cause-specific hazard; it does not by itself give the real-world probability of remaining free of all causes or the probability of relapse by 12 months.

The R extraction and hazard-ratio calculation can be kept explicit: `bh <- basehaz(fit_relapse, centered = FALSE)`; select the row with `time == 12` if that event time is present; then multiply its `hazard` value by `exp(coef(fit_relapse)["treatment"])` for the treatment comparison when age is at its reference value. If month 12 is not an event time, the extracted step-function estimate must be interpreted using the software's event-time convention rather than assuming a new jump exactly at month 12.

Practical Checks and Common Failure Modes

Before interpreting results, verify that the code implements the intended risk sets and model:

Tabulate `status` and each constructed indicator. For cause $k$, the indicator should equal 1 exactly when status equals $k$, and should be 0 for both other causes and censoring.
Check the number of target events and inspect warnings, convergence, and coefficient standard errors. Sparse causes or covariate groups with few events can produce unstable estimates or very wide intervals.
Confirm that all fitted causes use the same time origin, units, and covariate definitions when comparing their results.
Check factor reference levels, missing-data exclusions, and model-matrix coding before interpreting coefficients or a baseline curve.
Investigate whether the proportional hazards assumption and chosen covariate functional forms are plausible. Formal residual-based diagnostics are a separate step; a successful software fit alone does not establish these assumptions.

A frequent error is to code every observed event as the event indicator in every model. For a relapse model, setting the indicator to 1 for both relapse and death no longer estimates the relapse cause-specific hazard; it fits a composite endpoint. The correction is to create a cause-specific indicator for each fit and verify that competing events have indicator 0 while remaining in the risk set until their observed event times.

Separate cause-specific Cox models are useful when the goal is to describe or compare covariate associations with each cause-specific hazard. They are not interchangeable with a Fine–Gray subdistribution-hazard model, which targets a different quantity, and neither set of hazard ratios alone gives adjusted cumulative incidence. Estimating event probabilities requires combining information across causes; a cause-specific hazard ratio is not a risk ratio.

Key takeaway: Fit one Cox model per cause by making that cause the only event coded as 1. Competing events end follow-up for that cause-specific risk set, while the fitted hazard ratio describes a conditional instantaneous rate. Interpret `basehaz` only after identifying the covariate reference represented by its cumulative hazard curve.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A status variable is coded 0 for censoring, 1 for relapse, and 2 for death. What event indicator should be used for a death cause-specific model?
A fitted coefficient for a binary covariate is $-0.20$. Calculate its hazard ratio and describe its meaning under the model.
`basehaz(fit, centered = FALSE)` returns a cumulative hazard of 0.10. What must you know before calling this a baseline hazard at a particular patient profile?
Why is a competing event coded as 0 in a cause-specific model, and why is that not the same as treating it as an unobserved future event?

Solutions. (1) Use an indicator equal to 1 when status is 2 and 0 otherwise; relapse and censoring are both non-events for this death fit. (2) $\exp(-0.20)=0.819$, so the covariate is associated with an estimated 18.1% lower cause-specific hazard, holding other modeled covariates fixed, under proportional hazards. (3) Identify the covariate values represented by zero, including factor reference coding and any centering; the curve is not automatically for an arbitrary patient. (4) The person is at risk for the target cause until the competing event occurs, then leaves the risk set because a first target event can no longer occur. The competing event is observed and has its own type, unlike right censoring, which leaves the later event status unknown.