Diagnosing Proportionality in Cause-Specific Cox Models

Why Diagnose Proportionality?

A cause-specific Cox model summarizes a covariate’s association with the instantaneous rate of a particular event type among individuals who remain free of every event type. Its proportional-hazards assumption says that a covariate’s effect on that cause-specific hazard does not change over time. If the effect does change, a single hazard ratio can conceal an important pattern: it may average a strong early association with a weak, absent, or reversed later association.

Here the target event is cause $k$. Individuals experiencing a competing event leave the risk set at that time, as in the earlier material on cause-specific Cox models. The diagnostic question is whether the fitted covariate coefficients can reasonably be treated as constant over follow-up. A diagnostic cannot establish that the model is true; it can identify patterns that are inconsistent with the assumed form or flag aspects of the data that make the assessment inconclusive.

Definition: A cause-specific Cox model satisfies proportional hazards when its cause-specific hazard has the form $$ \lambda_k(t\mid X)=\lambda_{0k}(t)\exp(\beta^\mathsf{T}X), $$ where $\lambda_{0k}(t)$ is an unspecified baseline cause-specific hazard and the coefficient vector $\beta$ is constant over time. For a one-unit increase in covariate $X_r$, holding the other modeled covariates fixed, the cause-specific hazard ratio is $\exp(\beta_r)$ at every time $t$. A time-varying effect violates this constancy.

A useful diagnosis presumes a reasonably specified model apart from the proportionality question being examined. In particular, event causes and times must be recorded appropriately; observations must be independent; and censoring must be non-informative conditional on the modeled information in the usual Cox-model sense. The covariates must be represented adequately, since a misspecified covariate effect can produce residual patterns that look like a time trend. Competing events are handled as exits from the cause-specific risk set; the cause-specific hazard model does not require assuming independent latent failure times for the competing causes.

Schoenfeld Residuals and Their Time Pattern

For a distinct event time $t_j$ at which cause $k$ occurs, the raw Schoenfeld residual for a covariate is the observed covariate value for the individual with the event minus its model-based risk-set average. For a covariate vector $X$, that average is weighted by the fitted relative hazards. In the following expression, $R(t_j)$ is the risk set immediately before $t_j$, and $\hat\beta$ is the fitted coefficient vector.

r_j=X_j- \frac{\sum_{i\in R(t_j)}X_i\exp(\hat\beta^\mathsf{T}X_i)} {\sum_{i\in R(t_j)}\exp(\hat\beta^\mathsf{T}X_i)}

The second term is the covariate value expected at that event time under the fitted model, among those at risk. If the model’s constant coefficient adequately describes the covariate effect, residuals should not show a systematic trend with time. A scaled Schoenfeld residual is a rescaled version of this residual, using estimated model information so that its time pattern can be read on a coefficient-like scale. Software implementations may use different scaling details and time transformations; interpretation should follow the method used by the software.

The usual graphical check plots scaled Schoenfeld residuals against event time or a stated transformation of time, then adds a smooth curve and uncertainty band. A persistent curve away from a horizontal line suggests that the coefficient may vary with time. The curve is diagnostic, not a fitted effect model by itself. Wide bands, a few late events, or a sharp change in the number at risk can make its shape unreliable.

Justification: Conditional on the risk set at $t_j$, under the fitted cause-specific Cox model the event individual is selected with probability proportional to $\exp(\beta^\mathsf{T}X_i)$. Therefore the conditional mean of $X_j$ is the risk-set weighted average in the residual formula. At the true constant coefficient, the conditional mean of the raw residual is zero. A systematic association between the scaled residuals and time is therefore evidence against a constant coefficient, subject to sampling variation and correct specification of the rest of the model.

A formal test can summarize the time pattern by regressing scaled Schoenfeld residuals on a chosen function of time, or by testing their association with that function. A common choice is a rank-based time transformation; another is log time. A test for one covariate has one degree of freedom when it tests a single time trend. A global test assesses the corresponding trend across the model’s covariates jointly. The null hypothesis is no association between the residual pattern and the chosen time function; it is not a test that every possible form of time variation is absent.

Worked Example: Interpreting a Schoenfeld Test

Suppose a hypothetical cause-specific Cox fit includes treatment and age. For treatment, a scaled Schoenfeld test using the software’s stated rank-time transformation gives $\chi^2=6.25$, with $df=1$. The model was fitted to 240 hypothetical patients, with 58 events of the target cause. Assume the model and test were applied to the same correctly coded analysis sample.

The p-value is the upper-tail probability under a chi-square distribution with one degree of freedom:

p=P(\chi^2_1\geq 6.25)\approx 0.0124

At significance level $\alpha=0.05$, $0.0124<0.05$, so the test provides evidence that the treatment residual pattern is associated with the chosen time transformation. This is evidence against a constant treatment coefficient under the fitted model, not proof that treatment effects truly vary or an estimate of the form of that variation. The residual plot should be examined to see whether the departure is gradual, localized, or driven by a few events.

Testing a Specified Time Interaction

A second diagnostic fits an extended cause-specific Cox model containing a covariate-by-time interaction. For a covariate $X_r$ and a prespecified time function $g(t)$, include both the covariate’s main effect and its interaction with time. The interaction coefficient tests whether the effect changes according to that particular function. This approach is targeted: it can have good power for a trend like the one specified, but may miss a different pattern.

\lambda_k(t\mid X)=\lambda_{0k}(t) \exp\left\{\beta^\mathsf{T}X+\gamma X_r g(t)\right\}

The instantaneous log hazard ratio for a one-unit increase in $X_r$, with other covariates held fixed, is $\beta_r+\gamma g(t)$, and its hazard ratio is $\exp\{\beta_r+\gamma g(t)\}$. The null hypothesis $H_0:\gamma=0$ corresponds to proportionality for this covariate within the specified model. A Wald test uses $\hat\gamma/SE(\hat\gamma)$; a likelihood-ratio test compares the extended model with the model omitting the interaction. For a one-parameter interaction, the Wald statistic is compared with a chi-square distribution with one degree of freedom when its usual large-sample approximation is appropriate.

The time origin and scale matter to interpretation. For example, with $g(t)=\log(t/12)$, $t$ must be measured in months and the covariate effect $\beta_r$ is the log hazard ratio at 12 months, because $g(12)=0$. The main effect must remain in the model when adding the interaction. The chosen function should be scientifically interpretable and settled before examining many alternatives; trying many transformations and reporting only the most significant result inflates false-positive risk.

Worked Example: Testing a Treatment-by-Log-Time Interaction

In a hypothetical analysis, treatment is coded $X=1$ for treated and $X=0$ for control. A cause-specific model includes treatment, other prespecified covariates, and the interaction $X\log(t/12)$, where $t$ is measured in months. The fitted treatment coefficient is $\hat\beta=0.40$, and the interaction coefficient is $\hat\gamma=-0.30$ with $SE(\hat\gamma)=0.12$. Assume the same subjects, risk sets, and tie handling are used in the model with and without the interaction.

The Wald statistic for the interaction is

z=\frac{-0.30}{0.12}=-2.50,\qquad \chi^2=z^2=6.25

The two-sided Wald p-value is approximately $0.0124$. At $\alpha=0.05$, this provides evidence of a treatment effect that changes with log time, under this specified model. At 12 months, the treatment hazard ratio is

\exp\{0.40-0.30\log(12/12)\} =\exp(0.40)\approx 1.492

At 48 months, $\log(48/12)=\log(4)\approx1.386$, giving

\exp\{0.40-0.30(1.386)\} =\exp(-0.0158)\approx0.984

Thus this model estimates a cause-specific hazard ratio of about 1.49 at 12 months and about 0.984 at 48 months, conditional on the other covariates. These are model-based associations, not causal effects. The interaction is a particular time-varying specification; it does not establish that the log-time form is the correct description.

Reading Plots and Checking the Fit

Worked Example: A Residual Plot with Sparse Late Follow-Up

Consider a hypothetical scaled Schoenfeld plot for a binary biomarker in a cause-specific model. The smooth curve is near $0.50$ on the coefficient-like scale around month 3 and near $-0.05$ around month 36. Its band is fairly narrow early but widens after month 30. There are 40 target-cause events through month 24, but only 6 from months 24 through 36.

As a rough interpretation of the smooth, the implied hazard ratios at those two times are $\exp(0.50)\approx1.65$ and $\exp(-0.05)\approx0.951$. The downward curve is a potential proportionality concern: it suggests the biomarker association may weaken over time. However, the late estimate is based on fewer events and has greater uncertainty. The plot does not justify claiming a precise late effect or a true reversal unless the uncertainty and the fitted model support that conclusion.

The next checks are to confirm the event coding and the risk sets, inspect the number at risk and event counts over time, and assess whether an influential observation or an omitted covariate pattern could explain the curve. A formal interaction test may assess a prespecified trend, while the plot can reveal departures that a single chosen interaction function does not capture.

Good diagnostics combine evidence rather than treating any one p-value as a verdict. Inspect each covariate’s plot and test, then consider a global test for model-wide departures. A global test can identify that at least one coefficient may vary without identifying which one; separate tests help localize departures but introduce a multiple-testing concern. Confidence bands and event counts are essential context, especially near the tail of follow-up, where few individuals remain at risk.

Diagnostic	What it assesses	Main limitation
Scaled Schoenfeld plot	Shape of coefficient-related residual patterns over time	Smoothers and tail behavior can be unstable with few events
Schoenfeld time-trend test	Association with a specified time transformation	May have low power for departures unlike that transformation
Covariate-by-time interaction	A specified parametric pattern of changing effect	Results depend on the chosen time function and scale
Global proportionality test	Joint evidence of time trends across coefficients	Does not identify which covariate drives the signal

These checks address proportionality in a cause-specific hazard model. They are not interchangeable with diagnostics for the Fine-Gray subdistribution hazard: the models have different risk-set constructions and target different hazard quantities. A proportionality finding in one model does not establish it in the other.

Common Pitfalls and Practical Limits

Treating a non-significant result as proof of proportionality. A large p-value means the data did not provide strong evidence against the tested pattern; it does not rule out time variation. Limited events, especially late in follow-up, can make a test weak.
Reading every wiggle as a violation. Residuals vary randomly. Look for a sustained, interpretable pattern and its uncertainty, not isolated departures from zero.
Ignoring time scale and origin. With an interaction, changing the units or origin of time changes the coefficient interpretation. State the time function and units explicitly.
Confusing association with causation. A time-varying cause-specific hazard ratio remains a conditional association unless a separate causal design and assumptions justify causal interpretation.
Ignoring other model errors. Nonlinearity, omitted predictors, influential observations, or coding errors can also produce residual structure. A time trend should prompt model checking, not automatic adoption of an interaction.

Key takeaway: Use scaled Schoenfeld plots and tests to look for time structure, and use time interactions to test scientifically chosen forms of change. Interpret both in light of the cause-specific risk set, the model specification, event counts, and uncertainty.

Check Your Understanding

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

What does the null hypothesis of a test for a treatment-by-time interaction assert?
A Schoenfeld test gives $\chi^2=3.00$ with $df=1$ and $p=0.083$. What is an appropriate conclusion at $\alpha=0.05$?
For a model with $g(t)=\log(t/6)$, at what time is the treatment coefficient $\beta$ itself the log hazard ratio?
Why should a residual trend after month 30 be interpreted cautiously if only a few target-cause events occur after that time?

Solutions. 1. It asserts $\gamma=0$ for the specified time function, so the modeled treatment coefficient does not change according to that function. It does not rule out every other form of time variation. 2. Since $0.083>0.05$, fail to reject the no-trend null for this test. The result is not proof that proportionality holds. 3. At 6 months, because $\log(6/6)=0$. 4. The late smooth is informed by few events, so its uncertainty is greater and its shape may be unstable; inspect its band and event counts before interpreting the apparent trend.