Interactions in Cause-Specific Cox Models

Why interactions matter in cause-specific Cox models

A covariate’s association with the rate of a particular event may differ according to another predictor. For example, a treatment could have one cause-specific hazard ratio among patients with a risk factor and another among patients without it. An interaction term represents this difference in the regression model. It does not, by itself, establish that the association is causal or that the effect differs on every possible scale.

Here the outcome is a time \(T\) to the first event and its type \(J\). Fix a cause \(k\). The cause-specific hazard at time \(t\) is the instantaneous rate of a cause-\(k\) event among individuals who have not yet experienced any event just before \(t\). A cause-specific Cox model represents this rate using a baseline hazard and a covariate linear predictor. Individuals experiencing a different cause leave the event-free risk set at that time for estimation of this cause-specific hazard.

Definition: An interaction in a cause-specific Cox model is a product term in the linear predictor that allows a covariate’s cause-specific log hazard ratio to vary with another covariate. For a binary predictor \(X\) and a binary modifier \(Z\), the model is

\lambda_k(t\mid X,Z)=\lambda_{k0}(t)\exp\{\beta_X X+\beta_Z Z+\beta_{XZ}XZ\}.

Here \(\lambda_{k0}(t)\) is the cause-\(k\) baseline hazard, and \(X,Z\) are coded 0 or 1. The interaction coefficient \(\beta_{XZ}\) is the difference between the log hazard ratio for \(X\) at \(Z=1\) and that at \(Z=0\).

This interpretation depends on the chosen coding and reference values. With 0/1 coding, \(\beta_X\) is the log hazard ratio for changing \(X\) from 0 to 1 when \(Z=0\); it is not an average effect across both groups. Similarly, \(\beta_Z\) describes the contrast in \(Z\) when \(X=0\). The interaction changes the relevant contrast when the other predictor changes.

Constructing and interpreting interaction contrasts

For the binary-by-binary model, the hazard ratio comparing \(X=1\) with \(X=0\), within a fixed value of \(Z\), follows by subtracting the two linear predictors. At \(Z=0\), that contrast is \(\beta_X\); at \(Z=1\), it is \(\beta_X+\beta_{XZ}\). Thus the ratio of these two hazard ratios is \(\exp(\beta_{XZ})\). An interaction is multiplicative on the cause-specific hazard scale: it compares hazard ratios, not differences in hazards or cumulative incidences.

\frac{\lambda_k(t\mid X=1,Z=z)}{\lambda_k(t\mid X=0,Z=z)} =\exp(\beta_X+z\beta_{XZ}),\qquad \frac{HR_{X\mid Z=1}}{HR_{X\mid Z=0}}=\exp(\beta_{XZ}).

The same logic applies when predictors are continuous. If \(X\) is measured in units and \(Z\) is binary, the model term \(\beta_X X+\beta_ZZ+\beta_{XZ}XZ\) makes the log hazard ratio for a one-unit increase in \(X\), at \(Z=z\), equal to \(\beta_X+z\beta_{XZ}\). Rescaling \(X\), for example to units of ten, makes the resulting contrast more interpretable but does not change the model’s fitted values.

In general, a contrast \(c^\mathsf{T}\hat{\beta}\) has estimated variance \(c^\mathsf{T}\widehat{\mathrm{Var}}(\hat{\beta})c\). A confidence interval on the log hazard-ratio scale is the estimate plus or minus the appropriate standard-normal critical value times its standard error; exponentiating both endpoints gives an interval for the hazard ratio. This covariance calculation matters whenever a contrast combines more than one coefficient.

Worked Example: Binary Treatment and Binary Risk Group

Suppose a hypothetical cause-specific Cox fit for a target event uses \(X=1\) for treatment and \(Z=1\) for a baseline risk group. The estimates are \(\hat{\beta}_X=0.40\), \(\hat{\beta}_Z=-0.20\), and \(\hat{\beta}_{XZ}=-0.50\). Assume the model’s proportional-hazards and covariate-form assumptions are adequate for this illustration.

Among the reference risk group (\(Z=0\)), the treatment hazard ratio is \(\exp(0.40)=1.492\). Among the \(Z=1\) group, the treatment log hazard ratio is \(0.40-0.50=-0.10\), so its hazard ratio is \(\exp(-0.10)=0.9048\). The estimated ratio of these treatment hazard ratios is \(\exp(-0.50)=0.6065\).

Suppose the estimated standard error of \(\hat{\beta}_{XZ}\) is \(0.22\). The interaction test statistic is \(z=-0.50/0.22=-2.273\), giving a two-sided \(p\)-value of about \(0.0230\). A 95% confidence interval for the interaction coefficient is \(-0.50\pm1.96(0.22)=(-0.9312,-0.0688)\). Exponentiating gives a 95% interval of \((0.394,0.934)\) for the ratio of treatment hazard ratios. Under the fitted model, the estimated treatment hazard ratio is lower in the \(Z=1\) group on this multiplicative scale; the interval excludes a ratio of 1. This is an association with the cause-specific hazard, not a comparison of cumulative event probabilities.

Continuous predictors and covariance-aware inference

When a predictor is continuous, the reference point for the other variable still matters. Centering a continuous modifier at a meaningful value can make a main-effect coefficient easier to interpret: after centering, its zero value corresponds to that chosen reference point. Centering changes the meanings of some coefficients but leaves fitted hazard predictions unchanged, provided the model is reparameterized consistently.

Worked Example: A Continuous Predictor by Group Interaction

In a hypothetical model, let \(X\) be a biomarker measured in 10-unit increments and let \(Z=1\) identify a subgroup. The fitted terms are \(\hat{\beta}_X=0.18\), \(\hat{\beta}_Z=-0.30\), and \(\hat{\beta}_{XZ}=0.12\). The standard errors of \(\hat{\beta}_X\) and \(\hat{\beta}_{XZ}\) are \(0.20\) and \(0.15\), and their estimated covariance is \(-0.006\). Assume proportional hazards and an adequately specified linear effect of \(X\) on the log hazard in each group.

For \(Z=0\), a 10-unit increase in \(X\) has estimated cause-specific hazard ratio \(\exp(0.18)=1.197\). For \(Z=1\), the log hazard ratio is \(0.18+0.12=0.30\), giving \(\exp(0.30)=1.350\). The ratio of these hazard ratios is \(\exp(0.12)=1.128\).

For the \(Z=1\) log hazard-ratio contrast, the variance is \(0.20^2+0.15^2+2(-0.006)=0.0505\), so its standard error is \(\sqrt{0.0505}=0.2247\). A 95% interval on the log scale is \(0.30\pm1.96(0.2247)=(-0.1404,0.7404)\); after exponentiation, the hazard-ratio interval is \((0.869,2.097)\). It includes 1, so these data do not establish a nonzero log-hazard slope for \(X\) in the \(Z=1\) subgroup at the 5% level. Separately, the estimated interaction coefficient has standard error \(0.15\), so its 95% interval is \(0.12\pm1.96(0.15)=(-0.174,0.414)\), also including 0. Neither interval should be interpreted as a probability that the fixed parameter lies within it.

The test of \(\beta_{XZ}=0\) addresses whether the \(X\) slope differs between the two groups on the log cause-specific hazard scale. It is distinct from testing whether either group-specific slope is zero. In models with several interaction coefficients, the corresponding question is usually a joint test that all interaction coefficients are zero. The covariance matrix of the estimates is needed for such a test.

Interactions with categorical predictors and joint tests

For a categorical predictor with a reference level and multiple nonreference levels, create indicator variables for those levels and include the relevant product terms with the modifier. Each treatment contrast is then a combination of its indicator coefficient and its interaction coefficient. A single interaction coefficient tests one such difference; a joint test evaluates all interaction terms together.

Worked Example: Three Treatment Levels and a Joint Interaction Test

Consider a hypothetical cause-specific Cox model comparing treatments A, B, and C, with A as the reference, and a binary modifier \(Z\). Let the coefficients for B, C, \(Z\), B-by-\(Z\), and C-by-\(Z\) be \(0.35\), \(-0.10\), \(0.20\), \(-0.60\), and \(0.40\), respectively. Suppose the estimated variances of the two interaction coefficients are both \(0.04\), their covariance is zero, and the model assumptions are adequate.

At \(Z=0\), the hazard ratios relative to A are \(\exp(0.35)=1.419\) for B and \(\exp(-0.10)=0.9048\) for C. At \(Z=1\), the B-versus-A log hazard ratio is \(0.35-0.60=-0.25\), giving \(\exp(-0.25)=0.7788\). The C-versus-A log hazard ratio is \(-0.10+0.40=0.30\), giving \(\exp(0.30)=1.350\). These contrasts show why neither treatment’s reference-group coefficient alone describes its association in both modifier groups.

For the joint null hypothesis that both interaction coefficients are zero, the Wald statistic is \(W=(\hat{\beta}_{BZ}^2+\hat{\beta}_{CZ}^2)/0.04=(0.36+0.16)/0.04=13.00\). Under the null and the usual large-sample approximation, \(W\) is compared with a chi-square distribution with 2 degrees of freedom. The corresponding \(p\)-value is \(\exp(-13/2)=0.00150\). This provides evidence that at least one treatment contrast varies across \(Z\) on the multiplicative cause-specific hazard scale; it does not identify which contrast differs without examining separate estimates or tests.

The general Wald statistic for a vector of interaction estimates \(\hat{\gamma}\), with estimated covariance matrix \(\widehat V_\gamma\), is \(W=\hat{\gamma}^{\mathsf T}\widehat V_\gamma^{-1}\hat{\gamma}\). Under the joint null and regular large-sample conditions, it is approximately chi-square with degrees of freedom equal to the number of tested interaction parameters. A likelihood-ratio test can instead compare the model containing interactions with the nested model omitting them, using the difference in model deviances and the corresponding degrees of freedom.

Model conditions, diagnostics, and alternatives

Interaction interpretation relies on the model being appropriate for the target cause-specific hazard. In particular, the Cox model assumes that the covariate effects represented in the linear predictor do not vary with time, unless time-varying terms are included. It also assumes that the chosen functional forms adequately represent continuous predictors and that the risk-set and censoring assumptions needed for the fitted model are reasonable. An interaction does not repair a nonlinear main effect or a time-varying association.

Useful checks include examining proportional-hazards diagnostics for patterns that suggest time-varying coefficients, and checking whether continuous covariates have a plausible functional form within the relevant groups. Sparse failures in one combination of predictor levels can produce unstable estimates and wide intervals; cross-tabulate event counts and inspect coefficient standard errors and convergence warnings. If the groups have little covariate overlap, some contrasts may be estimated largely by model extrapolation. A small interaction \(p\)-value does not resolve these problems.

A model without interaction imposes a common log hazard ratio across modifier levels. The interaction model relaxes that restriction but estimates more parameters, so it can be less precise, especially with few events. Cause-specific Cox and Fine-Gray regression also target different hazard quantities: the former concerns the instantaneous cause-specific rate among those still event-free, while the latter models a subdistribution hazard. An interaction in one model is not automatically the same scientific contrast as an interaction in the other. Selecting between them depends on the estimand, not merely on which interaction test is significant.

Interactions are scale-dependent. A zero product coefficient means no modification of the log cause-specific hazard ratio (equivalently, a common multiplicative hazard ratio); it does not imply no modification of absolute risk or cumulative incidence. To communicate results, report the coding and units, the group-specific contrasts with uncertainty, and the interaction test when appropriate. Avoid presenting a main-effect coefficient as if it applied uniformly when an interaction is in the model.

Common mistake: interpreting a main effect without its reference

Incorrect reasoning: “The treatment coefficient is 0.40, so treatment multiplies the cause-specific hazard by \(\exp(0.40)\) in everyone.”
Correction: With a treatment-by-group product term and 0/1 coding, \(\exp(0.40)\) is the treatment hazard ratio only in the group coded \(Z=0\). For \(Z=1\), add the interaction coefficient before exponentiating.
Incorrect reasoning: “A significant interaction means treatment benefits one group.”
Correction: The interaction test detects a difference between log hazard ratios. Direction and magnitude require the group-specific estimates, and a cause-specific hazard association is not itself a causal effect or a cumulative-incidence difference.

A useful justification for the central interpretation is direct subtraction. Under the model, the log hazard at \(X=1\) minus the log hazard at \(X=0\), holding \(Z=z\) fixed, equals \(\beta_X+\beta_{XZ}z\), because the baseline hazard and the terms involving only \(Z\) cancel. Exponentiating gives the stratum-specific hazard ratio. Subtracting the \(Z=0\) contrast from the \(Z=1\) contrast leaves \(\beta_{XZ}\), and exponentiation therefore gives the ratio of hazard ratios. This derivation also clarifies why the interaction coefficient itself is not a hazard ratio.

Key takeaway: In a cause-specific Cox model, product terms allow covariate hazard-ratio contrasts to vary across values or levels of another predictor. Interpret each contrast from the full linear predictor, use its coefficient covariance for uncertainty, and state that the result concerns the cause-specific hazard scale.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

In a model with binary \(X\) and \(Z\), coefficients \(\beta_X=0.25\) and \(\beta_{XZ}=0.30\), what are the hazard ratios for \(X=1\) versus \(X=0\) at \(Z=0\) and \(Z=1\)?
What contrast does \(\exp(\beta_{XZ})\) represent in that model?
Why is \(\beta_X\) not a treatment effect that applies to both groups when the model includes \(XZ\)?
A two-parameter interaction test gives \(W=5.99\). What is its approximate \(p\)-value under a chi-square distribution with 2 degrees of freedom, and what does the test address?
Does an interaction in a cause-specific Cox model directly establish modification of cumulative incidence? Explain.

Solutions. (1) At \(Z=0\), the hazard ratio is \(\exp(0.25)=1.284\). At \(Z=1\), it is \(\exp(0.25+0.30)=\exp(0.55)=1.733\). (2) It is the ratio of the \(X\) hazard ratio at \(Z=1\) to that at \(Z=0\), equal to \(\exp(0.30)=1.350\). (3) \(\beta_X\) is the log hazard ratio at the reference value \(Z=0\); the product term changes the contrast at \(Z=1\). (4) For 2 degrees of freedom, the upper-tail probability is \(\exp(-5.99/2)\approx0.0500\). This tests whether both interaction coefficients are zero jointly, not whether every covariate association is zero. (5) No. It concerns modification on the cause-specific hazard scale; cumulative incidence also depends on the competing-event process.