Diagnosing Cause-Specific Cox Fits with Schoenfeld Residuals

Why Diagnose Proportional Hazards?

A cause-specific Cox model describes the instantaneous rate of one event type among individuals who are still alive and have not yet experienced any event. Its proportional-hazards assumption says that a covariate has a constant multiplicative effect on that rate over follow-up. A treatment might, however, have a strong early effect that weakens later. If that pattern is present, a single hazard ratio can obscure an important feature of the association.

Scaled Schoenfeld residuals provide a way to look for evidence that a fitted coefficient changes with time. In R, the survival package’s cox.zph function gives tests and plots based on these residuals. If a time trend is credible and scientifically relevant, a cause-specific Cox model can include a time-varying coefficient. This tutorial focuses on those diagnostics and that extension, rather than on cumulative-incidence prediction or subdistribution-hazard models.

Model, Assumptions, and Schoenfeld Residuals

Definition: For cause $k$, a cause-specific Cox model is $$ \lambda_k(t\mid X)=\lambda_{0k}(t)\exp\{X^\mathsf{T}\beta_k\}, $$ where $\lambda_k(t\mid X)$ is the instantaneous rate of cause $k$ at time $t$, conditional on remaining event-free until $t$; $\lambda_{0k}(t)$ is an unspecified baseline hazard; and $\beta_k$ is a vector of covariate coefficients. The proportional-hazards assumption is that each coefficient is constant over time.

To fit the cause-specific model for cause $k$, code an event of cause $k$ as the event and treat other causes as censored at their event times. Those individuals leave the risk set when the competing event occurs. This defines the cause-specific hazard; it does not require competing event types to be independent. The usual Cox assumptions still matter: subjects are independent (or the dependence is handled appropriately), loss to follow-up is noninformative conditional on modeled information, the covariate effects are specified adequately, and the risk sets and event times are recorded correctly.

At an event time $t_j$ for cause $k$, a Schoenfeld residual compares the covariate value of the individual who experienced that event with the covariate value expected under the fitted Cox model among those at risk. For covariate $r$, its unscaled form is

r_{jr}=X_{i(j),r}- \frac{\sum_{i\in R(t_j)}X_{ir}\exp(X_i^\mathsf{T}\hat{\beta}_k)} {\sum_{i\in R(t_j)}\exp(X_i^\mathsf{T}\hat{\beta}_k)}.

Here $i(j)$ is the individual experiencing the cause-$k$ event, $R(t_j)$ is the risk set immediately before $t_j$, and $\hat{\beta}_k$ is the fitted coefficient vector. The residual is defined at event times, not at censored observations. Scaling adjusts residuals using information from the fitted model so their pattern can be interpreted on an approximate coefficient scale. Under proportional hazards, systematic drift in scaled residuals over time is not expected.

Worked Example: Calculating a Schoenfeld Residual

Suppose a hypothetical cause-specific Cox fit has one binary covariate, treatment, with fitted coefficient $\hat{\beta}=\log(2)$. At one event time, four individuals are at risk: two treated and two untreated. A treated individual experiences the event of interest. Calculate the unscaled treatment residual.

The model gives each treated individual relative risk $\exp\{\log(2)\}=2$, and each untreated individual relative risk $1$. The risk-set weighted mean treatment value is therefore $$ \frac{2(1)+2(0)+2(1)+2(0)}{2+1+2+1}=\frac{4}{6}=0.667. $$

The observed treatment value is $1$, so the unscaled residual is $1-0.667=0.333$. It is positive because the event occurred in a treated individual whose covariate value exceeded the risk-set weighted expectation. One residual alone does not diagnose nonproportional hazards; the relevant evidence is the pattern across event times.

Using cox.zph and Reading Its Plots

Fit the cause-specific model with an event indicator for the cause of interest. For example, if status is coded 0 for censoring and 1, 2, and 3 for event types, a model for cause 1 could be fit with the R command coxph(Surv(time, status == 1) ~ treatment + age + stage, data = d, x = TRUE). Competing events are coded as non-events in this fit and leave the risk set at their observed times. Retaining the model matrix with x = TRUE supports the subsequent diagnostic calculation.

Apply cox.zph(fit, transform = "km") to the fitted model. This reports a test for each model term and a global test. The default time transformation is Kaplan–Meier based; other available choices include rank and identity transformations. The test assesses whether the scaled residuals have a systematic linear association with the chosen transformed event time. For a one-degree-of-freedom term, it is commonly reported as a chi-square statistic with one degree of freedom. A multi-level factor can have a multi-degree-of-freedom test. The global test assesses departures across the model terms jointly.

The plot from plot(cox.zph_object) displays scaled Schoenfeld residuals against transformed time, usually with a smooth curve. Under proportional hazards, the curve should be roughly horizontal; a sustained rise or fall suggests a changing coefficient. The test and plot answer related but not identical questions: the test summarizes association with the selected time transformation, while the smooth can reveal a pattern that is nonlinear or localized. Inspect the curve and its uncertainty rather than treating a p-value as a complete diagnosis.

Worked Example: Interpreting a Proportional-Hazards Diagnostic

In a hypothetical study with 240 participants and 86 events of the cause of interest, a cause-specific model includes treatment, age, and stage. The output from cox.zph using the Kaplan–Meier transformation is: treatment, $\chi^2=8.12$, $df=1$, $p=0.0044$; age, $\chi^2=1.10$, $df=1$, $p=0.294$; stage, $\chi^2=0.55$, $df=1$, $p=0.458$; global, $\chi^2=10.3$, $df=3$, $p=0.016$.

At a significance level of $0.05$, the treatment test provides evidence against a constant treatment coefficient over time, and the global test provides evidence of some departure among the modeled terms. The age and stage tests do not show statistically significant evidence of a time trend in this analysis. They do not establish that those effects are proportional: limited event information can make a test insensitive.

Suppose the treatment plot also shows a smooth curve that declines over follow-up, with uncertainty widening near the end where few events remain. The plot supports investigating a changing treatment effect, but the tail should be interpreted cautiously because it is based on sparse information. The appropriate conclusion is that the fitted proportional-hazards treatment effect is doubtful, not that the diagnostic proves a particular biological explanation.

Modeling a Time-Varying Coefficient

A direct extension replaces a constant coefficient by a function of time. For a binary treatment $Z$, one useful specification is

\lambda_k(t\mid Z,X)=\lambda_{0k}(t) \exp\left\{Z\left[\beta_0+\beta_1\log(t/t_0)\right]+X^\mathsf{T}\gamma\right\},

where $t_0$ is a chosen reference time, and the other covariates $X$ retain constant coefficients $\gamma$. The treatment hazard ratio at time $t$, comparing $Z=1$ with $Z=0$ for otherwise equal covariates, is $\exp\{\beta_0+\beta_1\log(t/t_0)\}$. At $t=t_0$, this ratio is $\exp(\beta_0)$. A positive $\beta_1$ means the log hazard ratio increases with log time; it does not necessarily mean the hazard ratio exceeds 1 at every time.

With follow-up measured in days, an R specification using $t_0=365$ days is coxph(Surv(time, status == 1) ~ treatment + tt(treatment) + age + stage, data = d, tt = function(x, t, ...) x * log(t / 365)). The term treatment estimates the effect at 365 days, and tt(treatment) estimates the change per one-unit increase in $\log(t/365)$. The time function must be defined on the observed event times; if time is zero or the chosen function is undefined there, choose a suitable positive reference and a specification appropriate to the time origin.

The model assumes that this chosen time pattern is adequate and that the remaining terms satisfy proportional hazards. A time-varying coefficient is not a general cure for misspecification: inspect whether the fitted functional form captures the diagnostic pattern, check influential observations and sparse risk sets, and consider whether continuous covariates were modeled appropriately. cox.zph is not a general diagnostic for models containing tt() terms; assess the extended model using the specified interaction, its uncertainty, and suitable checks for remaining departures.

Worked Example: Interpreting a Time-Varying Treatment Effect

Suppose a hypothetical extended cause-specific Cox model uses $t_0=365$ days and estimates $\hat{\beta}_0=-0.20$ and $\hat{\beta}_1=0.35$. Thus the treatment log hazard ratio at time $t$ is $-0.20+0.35\log(t/365)$. Calculate and interpret the treatment hazard ratio at 365 and 730 days.

At 365 days, $\log(365/365)=0$, so the log hazard ratio is $-0.20$, and the hazard ratio is $\exp(-0.20)=0.819$. At that time, the estimated cause-specific hazard is about 18% lower for treated than untreated participants with the same values of the other covariates.

At 730 days, $\log(730/365)=\log(2)=0.693$. The log hazard ratio is $-0.20+0.35(0.693)=0.0426$, giving a hazard ratio of $\exp(0.0426)=1.044$. The fitted association therefore changes from an estimated lower cause-specific hazard at 365 days to a slightly higher one at 730 days. These are conditional cause-specific hazard comparisons, not probabilities of experiencing the event.

For illustration, suppose the estimated variances are $\operatorname{Var}(\hat{\beta}_0)=0.0400$, $\operatorname{Var}(\hat{\beta}_1)=0.0225$, and $\operatorname{Cov}(\hat{\beta}_0,\hat{\beta}_1)=-0.0120$. At 730 days, the estimated variance of the log hazard ratio is $$ 0.0400+0.693^2(0.0225)+2(0.693)(-0.0120)=0.0342, $$ so its standard error is $\sqrt{0.0342}=0.185$. A 95% Wald confidence interval for the log hazard ratio is $0.0426\pm1.96(0.185)=(-0.320,\,0.405)$. Exponentiating gives a 95% confidence interval of approximately $0.726$ to $1.50$ for the 730-day hazard ratio. This interval includes 1, so the estimate at that time is not precise evidence of a difference in cause-specific hazards.

To test whether the time interaction adds evidence beyond the constant-effect model, compare the nested models using a likelihood-ratio test, provided both use the same observations and compatible tie handling. If their partial log likelihoods are $\ell_0$ and $\ell_1$, the statistic is $2(\ell_1-\ell_0)$. For a hypothetical statistic of 6.80 with one additional coefficient, the reference distribution is chi-square with one degree of freedom and $p\approx0.0091$. This supports the added time interaction, subject to the model and testing assumptions; it does not validate every aspect of the time-varying model.

Choosing a Response and Recognizing Failure Modes

A time interaction is useful when the time pattern itself is important and can be represented sensibly. A stratified Cox model is another option when a categorical variable violates proportional hazards but is a nuisance: it allows separate baseline hazards across strata, but does not estimate a single coefficient for that stratifying variable. These approaches answer different questions. A plotted trend may also suggest a different time function, but trying many functions and reporting only the most favorable test can exaggerate evidence; specify and justify the function, and report the exploration.

A cause-specific Cox model addresses the cause-specific hazard, not the cumulative probability of that cause by a given time. In particular, a hazard ratio below 1 does not directly imply a particular reduction in cumulative incidence, which also depends on the other causes. The earlier material on cause-specific hazards gives the surrounding interpretation; here the focus is whether the modeled coefficient is stable over time.

Do not read a nonsignificant test as proof. The test may have low power, especially with few cause-specific events or little late follow-up.
Do not diagnose from one residual. Schoenfeld residuals vary around their expectation; assess the pattern across event times and its uncertainty.
Do not equate a trend with causation. A changing fitted association does not establish why the association changes.
Check the risk sets and covariate form. Coding errors, omitted nonlinear effects, influential events, and sparse time regions can distort diagnostics.
Do not overinterpret the tail. Few late events yield unstable smooths and time-specific estimates.

Key takeaway: Scaled Schoenfeld residuals help assess whether cause-specific Cox coefficients are stable over time. When a credible departure matters, specify and interpret a time-varying coefficient explicitly, with uncertainty, rather than reporting a single hazard ratio as if it applied throughout follow-up.

Check Your Understanding

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

In a cause-specific Cox fit for cause 1, what happens to a participant’s risk-set membership when cause 2 occurs?
A one-degree-of-freedom cox.zph test gives $p=0.18$. Does that establish proportional hazards? Why or why not?
In the model with treatment coefficient $\beta_0+\beta_1\log(t/t_0)$, what does $\exp(\beta_0)$ represent?
Name one limitation of using a time-varying coefficient and one alternative when a categorical nuisance variable is nonproportional.

Solutions. (1) The participant leaves the risk set at the cause-2 event time; that event is treated as a non-event for cause 1. (2) No. It indicates that this test did not find statistically significant evidence against proportionality at the chosen level, but limited power or an untested pattern may remain. (3) It is the treatment cause-specific hazard ratio at the reference time $t_0$, conditional on the other covariates. (4) The chosen time function may not capture the true pattern, and the added coefficient can be imprecise; stratification allows separate baseline hazards for categories but does not estimate a coefficient for the stratifying variable.