Diagnosing Proportional Hazards

Why Diagnose Proportional Hazards?

A Cox model summarizes the effect of a covariate with a hazard ratio that is constant over time, conditional on the other covariates. That constancy is the proportional-hazards (PH) assumption. It is central to interpreting a fitted Cox model: if a treatment hazard ratio is reported as 0.70, the model treats that relative hazard as 0.70 throughout follow-up. A strong time pattern in the treatment effect would make that single summary inadequate.

Diagnostics do not prove that the PH assumption is true. They look for patterns that would be unexpected under a correctly specified proportional-hazards model. Schoenfeld residuals and their scaled versions focus on whether covariate effects appear to change with time. A time-interaction test provides a related, more targeted check. Neither diagnostic, by itself, establishes that the rest of the model is correctly specified.

Model and Diagnostic Setup

For a subject with covariate vector $x_i$, the Cox proportional-hazards model is

h_i(t\mid x_i)=h_0(t)\exp(x_i^\top\beta).

Here $h_i(t\mid x_i)$ is the subject's hazard at time $t$, $h_0(t)$ is an unspecified baseline hazard, and $\beta$ is a vector of constant covariate effects. For two subjects whose covariate vectors differ by one unit in a single covariate $x$, the model implies a hazard ratio $\exp(\beta_x)$ at every time, conditional on the other covariates.

The diagnostic procedures below presume that the model is fitted to the relevant event-time and right-censoring data, with the risk sets and ties handled by the chosen Cox fitting method. Usual requirements for interpreting the fitted model also matter: censoring should be non-informative conditional on the modeled covariates, subjects should be independent (or dependence handled appropriately), and covariate effects should be represented adequately apart from the particular PH question being assessed. For example, a badly chosen nonlinear functional form for a continuous covariate can produce residual patterns that are not evidence of a changing effect over time.

Definition: At an event time $t_j$, a raw Schoenfeld residual for covariate $x$ is the observed covariate value for the subject who experiences the event minus the risk-set weighted mean of that covariate under the fitted Cox model. With one event at $t_j$, $$ r_j=x_{(j)}-\frac{\sum_{i\in R(t_j)} x_i\exp(x_i^\top\hat\beta)} {\sum_{i\in R(t_j)}\exp(x_i^\top\hat\beta)}, $$ where $R(t_j)$ is the set at risk immediately before $t_j$, and $x_{(j)}$ is the event subject's covariate value. For tied events, the residual calculation must use a definition consistent with the method used to handle ties.

The residual compares the event subject with the covariate values that the fitted model expected among those still at risk. A residual above zero means the event subject's covariate value was greater than that weighted risk-set average; a residual below zero means it was smaller. One residual is noisy, so the diagnostic question is whether residuals show a systematic pattern across event times.

Why Schoenfeld Residuals Are Informative

The link to the Cox model's partial likelihood explains why these residuals are useful. For simplicity, suppose event times are distinct. At each event time, the partial-likelihood score contribution is the event subject's covariate value minus its fitted risk-set weighted mean. Summing those contributions gives the partial-likelihood score.

U(\beta)=\sum_{j:\,\text{event}}\left[ x_{(j)}- \frac{\sum_{i\in R(t_j)}x_i\exp(x_i^\top\beta)} {\sum_{i\in R(t_j)}\exp(x_i^\top\beta)} \right] =\sum_{j:\,\text{event}}r_j.

At an interior maximum of the partial likelihood, the fitted score is zero (up to numerical tolerance). Thus the residuals collectively balance at the fitted coefficients. This does not mean that each residual is zero, or that their ordering over time is irrelevant. Residuals that tend to be positive early and negative late, for example, can suggest that the model's constant effect is too small early and too large late. The total can still be close to zero because positive and negative residuals cancel.

Worked Example: Calculating a Raw Schoenfeld Residual

Suppose that at one event time four subjects are at risk, with covariate values $x=(0,1,2,1)$. For this illustrative calculation, take the fitted coefficient to be $\hat\beta=\log(2)$, and suppose the subject with $x=2$ experiences the event. There are no tied events at this time.

The fitted risk weights are $\exp(\hat\beta x_i)=2^{x_i}$, or $1,2,4,2$. The weighted mean covariate value in the risk set is

\frac{0(1)+1(2)+2(4)+1(2)}{1+2+4+2} =\frac{12}{9}=1.333.

The raw Schoenfeld residual is therefore

$$ r=2-1.333=0.667. $$

At this event time, the event subject's $x$ value is 0.667 units above the fitted risk-set weighted mean. This single positive residual does not establish a time-varying effect. The diagnosis depends on how residuals behave across many event times.

Scaled-Residual Plots and Formal Tests

Raw residuals are contrasts on the covariate scale. Scaled Schoenfeld residuals apply a covariance-based rescaling so that their values can be interpreted on a coefficient scale. In a PH model, a persistent trend in these scaled residuals against event time is evidence that a constant coefficient may not describe the data well. Software may differ in details of scaling and time transformation, so the plot's axes and test definition should be checked before interpreting its numerical values.

A common diagnostic displays scaled residuals against event time, or against a specified transformation of time, together with a smooth curve and uncertainty bands. Under PH, the smooth should not show a systematic trend away from a horizontal reference. Scatter around that reference is expected: residuals are variable, and their arrangement can be affected by the changing risk sets. A formal test commonly assesses whether the scaled residuals are associated with a chosen function of time. A test for one covariate asks whether its effect shows evidence of time dependence; a global test assesses departures across the modeled covariates jointly.

A small p-value is evidence against the tested PH specification, not a probability that the PH assumption is false. A large p-value does not prove PH: the test may have limited power, particularly with few events or a departure that the selected time transformation represents poorly. Plots are therefore important alongside tests. If many covariates are tested separately, chance small p-values become more likely; report the scope of the checks and avoid treating one isolated result as conclusive.

Worked Example: Reading a Scaled-Residual Trend

Suppose a hypothetical fitted Cox model has scaled residual summaries for a treatment indicator at five increasing transformed event times. The summaries, in coefficient-scale units, are $-0.42,-0.18,0.04,0.21,0.35$, at transformed times $0,1,2,3,4$, respectively. These are illustrative summaries, not a real study.

A simple descriptive line through these points has slope

\frac{\sum_{j=1}^5 (t_j-\bar t)(r_j-\bar r)} {\sum_{j=1}^5 (t_j-\bar t)^2} =\frac{1.93}{10}=0.193

coefficient-scale units per transformed-time unit. Here $\bar t=2$, $\bar r=0$, and the numerator is $(-2)(-0.42)+(-1)(-0.18)+(1)(0.21)+(2)(0.35)=1.93$.

The upward pattern is a reason to investigate possible time dependence; it is not a formal significance test. Five points cannot establish that the population effect changes, and the plotted summaries may have substantial uncertainty. The appropriate next check is the model's full residual plot and, if specified in a defensible way, a formal time-association test. Other misspecification, including an incorrect functional form for treatment-related covariates or influential observations, should also be considered.

Time Interactions as a Targeted Diagnostic

A second way to assess PH is to add a covariate-by-time term and test whether its coefficient is zero. For a covariate $x$ and a prespecified time function $g(t)$, the diagnostic extension is

h_i(t)=h_0(t)\exp\{x_i^\top\beta+\gamma x_i g(t)\}.

This expression is used here to formulate a diagnostic test, not as a general strategy for choosing a model after inspecting the data. When $g(t)$ is not constant over event times, $\gamma=0$ corresponds to no interaction with that chosen time function. A Wald, score, or likelihood-ratio test can assess $H_0:\gamma=0$, with the test and its standard error calculated under the fitted model. A failure to reject means there is insufficient evidence for this particular interaction; it does not rule out every possible time pattern.

The time function should have a substantive or diagnostic rationale and be chosen with care. For example, $g(t)=\log(t)$ tests for a smooth linear trend in log time, while a prespecified indicator for a time period tests a particular change across that division. Centering $g(t)$ can make the coefficient on $x$ interpretable at the centering time. Trying many time functions and reporting only the most significant one exaggerates evidence unless the exploration is clearly identified and inference accounts for it.

Worked Example: Testing a Prespecified Time Interaction

In a hypothetical Cox analysis of 240 subjects with 96 observed events, suppose a treatment indicator has fitted coefficient $\hat\beta_x=0.50$. A prespecified diagnostic interaction uses $g(t)=\log(t)$, with time measured in years, and gives $\hat\gamma=-0.12$ with $SE(\hat\gamma)=0.05$. The test is two-sided at $\alpha=0.05$.

The Wald statistic is

z=\frac{\hat\gamma}{SE(\hat\gamma)} =\frac{-0.12}{0.05}=-2.40.

Using a standard normal reference distribution, the two-sided p-value is

p=2\{1-\Phi(2.40)\}\approx 0.0164.

Because $0.0164<0.05$, this test provides evidence against a constant treatment effect in the particular log-time interaction direction tested. The estimate says that the fitted log-hazard contrast changes by $-0.12$ for each one-unit increase in $\log(t)$, under this diagnostic specification. It does not by itself establish the right model for that changing effect, nor does it establish a causal explanation. The time scale, chosen function, censoring assumptions, and other aspects of model fit still matter.

Distinguishing Time Variation From General Misspecification

A residual trend or significant time-interaction test is a warning about the fitted model, but its source is not automatically a genuinely time-varying covariate effect. A continuous covariate modeled linearly when its relationship is nonlinear can leave systematic residual structure. Omitted covariates, influential observations, dependence between subjects, or an unsuitable time origin can also complicate interpretation. Diagnostics should therefore be read as evidence about the model as a whole, not as a mechanical label attached to one assumption.

A practical assessment asks several questions in sequence:

Inspect the residual pattern.
Look for sustained trends over time rather than isolated large residuals. Note whether the pattern is concentrated in one covariate or appears broadly across the model.

Check the diagnostic specification.
Confirm the time scale, the transformation used for the plot or test, the handling of ties, and whether the interaction function was specified before examining the result.

Consider other model problems.
Review covariate functional forms, influential observations, dependence, and the plausibility of the censoring assumption. A time trend alone does not identify which problem is responsible.

Report evidence with its limits.
State which diagnostics were used, what pattern or test result was observed, and that failure to detect a departure is not proof of proportional hazards.

Common Mistakes and Why the Diagnostics Matter

The PH assumption determines whether a single hazard ratio is a defensible summary over follow-up. Diagnosing it can therefore change how a Cox result should be communicated. Even when the data do not provide strong evidence against PH, plots and tests make the scope of that conclusion clearer: the evidence concerns the observed follow-up, available events, chosen time transformations, and fitted covariates.

Mistake: treating a non-significant test as proof of PH. A large p-value is a failure to find evidence against the tested pattern, not evidence that all time variation is absent. Consider power and inspect the plot.
Mistake: treating any residual fluctuation as a violation. Residuals are variable. Look for systematic structure and uncertainty, rather than demanding a perfectly flat sequence.
Mistake: interpreting a significant test as proof of a time-varying causal effect. The test identifies evidence against the specified constant-effect model. Functional-form errors, influential data, or other misspecification can also contribute.
Mistake: testing many time functions and presenting the smallest p-value as if it were prespecified. Such searching can inflate false-positive evidence. Disclose exploratory choices and interpret the result cautiously.

Key takeaway: Schoenfeld residuals compare each event subject with the fitted risk set. A systematic pattern in scaled residuals over time, or a supported prespecified time interaction, is evidence against the corresponding proportional-hazards specification. Neither result alone proves the source of misspecification, and no diagnostic proves the assumption holds.

Check Your Understanding

Answer each question by connecting the diagnostic to its assumptions and interpretation.

At an event time, what two quantities are contrasted to form a raw Schoenfeld residual?
Why can the sum of Schoenfeld residuals be approximately zero even when the residual plot shows a time pattern?
A scaled-residual test has $p=0.41$. What can you conclude, and what can you not conclude, about proportional hazards?
In a prespecified interaction test using $g(t)=\log(t)$, what does $H_0:\gamma=0$ assess?
Give two reasons, other than a truly time-varying covariate effect, that might produce a residual pattern over time.