Tutorials › Multi-Omics Data Integration › Diagnosing Proportional Hazards In Omics Models

Clinical Outcomes: Survival and Longitudinal Data · Tutorial 436 of 597

Diagnosing Proportional Hazards In Omics Models

You will learn what each proportional-hazards diagnostic can show, how to calculate and interpret simple examples, and how to distinguish evidence of nonproportionality from other model problems.

Advanced 10 min read

What You'll Learn

  • State the proportional-hazards assumption for an omics-augmented Cox model.
  • Interpret scaled Schoenfeld residual patterns and their time-trend tests.
  • Use log-minus-log survival plots to assess whether group curves are approximately parallel.
  • Test a prespecified covariate-by-time interaction as a diagnostic.
  • Recognize censoring, model-form, and multiple-testing issues that can mislead diagnostics.

Why Check Proportional Hazards?

A Cox model can summarize how omics measurements and clinical covariates relate to an event rate while leaving the baseline hazard unspecified. Its standard interpretation, however, depends on the proportional-hazards (PH) assumption: the hazard ratio for a given covariate contrast stays constant over follow-up. If that assumption is poor, a single coefficient may hide an effect that changes over time. Diagnostics help identify this problem; they do not by themselves establish its cause.

For an omics-augmented model, a changing pattern could reflect a genuinely time-dependent association, but it could also arise from an incorrectly specified covariate effect, influential observations, sparse risk sets, or chance. The aim is therefore to look for consistent evidence using complementary checks, then interpret it in light of the study design and model specification.

Definition: A Cox model with covariate vector \(z\) has hazard \(h(t\mid z)=h_0(t)\exp(\beta^\top z)\), where \(h_0(t)\) is an unspecified baseline hazard and \(\beta\) is a vector of coefficients. The PH assumption requires that each coefficient is constant over time. For two covariate vectors \(z_a\) and \(z_b\), the hazard ratio is \(\exp\{\beta^\top(z_a-z_b)\}\), which does not depend on \(t\).

As in the earlier material on the Cox partial likelihood, estimation uses comparisons among subjects still at risk at observed event times. The diagnostics below require correctly constructed risk sets and a meaningful time origin. They also rely on censoring being independent of the event time conditional on the variables in the analysis, or on an appropriate account of the censoring and sampling design. These conditions do not guarantee PH; they make the fitted model and its diagnostics interpretable.

Scaled Schoenfeld Residuals

At an observed event time \(t_j\), the Schoenfeld residual for the subject who experiences the event is the observed covariate vector minus its risk-set average under the fitted model. With no tied events, its unscaled form is

$$ r_j=z_{(j)}-\frac{\sum_{i\in R(t_j)}z_i\exp(\hat\beta^\top z_i)} {\sum_{i\in R(t_j)}\exp(\hat\beta^\top z_i)}. $$

Here \(R(t_j)\) is the set at risk just before \(t_j\), and \(z_{(j)}\) is the event subject's covariate vector. The risk-set average is weighted toward subjects with higher fitted relative hazard. Tied-event handling depends on the partial-likelihood method used; residuals should be obtained consistently with that fit.

Scaled Schoenfeld residuals adjust these event-time residuals using the estimated coefficient covariance, placing them on a scale related to the coefficient. Software commonly plots each covariate's scaled residual against event time or a transformation of time and may add a smooth trend. Under PH, there should be no systematic trend in the residuals over time. Random scatter around zero is expected; individual points need not equal zero.

Worked Example: Reading a Schoenfeld Residual

Suppose a hypothetical fitted Cox model has one binary omics indicator \(x\), with \(\hat\beta=\log(2)\). At an event time, the risk set contains three subjects with \(x=(0,1,1)\), and the subject with \(x=1\) has the event. Calculate the unscaled Schoenfeld residual.

The fitted relative hazard weights are \(\exp(\hat\beta x)\), giving weights \(1,2,2\). The weighted risk-set mean is

$$ \frac{0(1)+1(2)+1(2)}{1+2+2}=\frac{4}{5}=0.800. $$

The event subject's observed value is 1, so \(r_j=1-0.800=0.200\). This positive residual says that the event subject's \(x\) value is above the model-weighted risk-set average at this event time. One positive residual is not evidence of a PH violation. The diagnostic question is whether scaled residuals show a systematic pattern across event times.

In a hypothetical residual plot, suppose early scaled residuals are mostly above zero and late residuals mostly below zero. That pattern would suggest that the fitted coefficient may not represent a constant effect. The plot is a prompt for further checking, not proof: sparse late risk sets or a few influential events could produce a similar appearance.

A common formal companion is the Grambsch–Therneau test, which assesses association between scaled Schoenfeld residuals and a chosen transformation of event time. The null hypothesis is no systematic time trend for the tested coefficient; a global version tests the coefficients jointly. The chosen time transformation affects sensitivity to different trend shapes. A test can flag departures without explaining them, and a nonsignificant result does not establish that PH holds.

Log-Minus-Log Plots

A log-minus-log plot compares estimated survival curves for groups. For a categorical covariate, estimate each group's survival function \(\hat S(t)\), then plot \(\log[-\log\{\hat S(t)\}]\) against \(\log(t)\), for positive event times. Under PH, the transformed curves should be approximately parallel: they can be separated, but their separation should remain roughly constant over time.

The reason follows from the model. For a group with fixed covariate vector \(z\), the Cox model implies \(S(t\mid z)=S_0(t)^{\exp(\beta^\top z)}\), where \(S_0(t)\) is baseline survival. Taking the specified transformation gives

$$ \log[-\log\{S(t\mid z)\}] =\log[-\log\{S_0(t)\}]+\beta^\top z. $$

The first term is common to all groups. The covariate contrast shifts the transformed curve by a constant, which is why nonparallel shapes can indicate nonproportional hazards. The argument also explains a limitation: curves based on Kaplan–Meier estimates are usually unadjusted. For continuous omics measurements, categorizing values for a plot discards information and can make results depend on arbitrary cut points. Use this plot as a visual diagnostic, not as the primary analysis of a continuous feature.

Worked Example: Comparing Log-Minus-Log Curves

Consider hypothetical Kaplan–Meier estimates for low and high expression groups, with independent censoring within each group. At 2 months, suppose estimated survival is 0.80 for the low group and 0.90 for the high group. At 4 months, the estimates are 0.60 and 0.75, respectively. Calculate the transformed values and compare their separation.

At 2 months, the low-group value is \(\log[-\log(0.80)]=\log(0.2231)\approx-1.500\); for the high group it is \(\log[-\log(0.90)]=\log(0.1054)\approx-2.250\). The gap, low minus high, is approximately 0.750.

At 4 months, the corresponding values are \(\log[-\log(0.60)]=\log(0.5108)\approx-0.672\) and \(\log[-\log(0.75)]=\log(0.2877)\approx-1.246\). The gap is approximately 0.574. The gap changes by about \(0.750-0.574=0.176\), so these two points are not perfectly parallel. This small example alone cannot establish a violation: sampling variability, especially with few people at risk, can alter the apparent gap. The complete plot and risk-set sizes matter.

Time-Interaction Checks

A prespecified time interaction can turn a visual suspicion into a targeted test. For a covariate \(x\), one diagnostic model adds \(xg(t)\) to the linear predictor:

$$ h(t\mid x,w)=h_0(t)\exp\{\beta x+\gamma xg(t)+\theta^\top w\}. $$

Here \(w\) denotes other covariates and \(g(t)\) is a chosen time function, such as \(\log(t/t_0)\) for a positive reference time \(t_0\). The PH null for this interaction is \(H_0:\gamma=0\); evidence against it indicates that the association represented by \(x\) changes with the chosen time function, conditional on the other model terms. This test can miss departures that do not resemble the selected \(g(t)\). Choosing the function after inspecting many possibilities also makes the nominal p-value overly optimistic unless that search is accounted for.

Worked Example: Testing a Time Interaction

In a hypothetical model, time is measured in months and \(t_0=12\). The fitted expression for a one-unit increase in an omics score is \(\hat\beta x+\hat\gamma x\log(t/12)\), with \(\hat\beta=0.80\), \(\hat\gamma=-0.35\), and \(SE(\hat\gamma)=0.16\). Assume the model and its standard-error approximation are suitable for a Wald test. Test the interaction null using a two-sided 5% test.

The Wald statistic is \(z=\hat\gamma/SE(\hat\gamma)=-0.35/0.16=-2.188\). Its two-sided normal-reference p-value is approximately 0.0287. Since \(0.0287<0.05\), this is evidence against \(\gamma=0\) for the specified time function; it is not proof of a particular biological mechanism.

The fitted log hazard ratio for a one-unit increase in \(x\), at time \(t\), is \(0.80-0.35\log(t/12)\). At 12 months the estimated hazard ratio is \(\exp(0.80)=2.23\). At 48 months it is \(\exp\{0.80-0.35\log(4)\}=\exp(0.315)\approx1.37\). These are model-based estimates at the stated times, not fixed hazard ratios across all follow-up. The calculation illustrates why the PH interpretation is questionable when this interaction is supported.

Combining Evidence and Diagnosing Failure

The three approaches answer related but different questions. Schoenfeld residuals use event-level information and can highlight which fitted coefficients have time trends. Log-minus-log plots offer an accessible visual check for grouped predictors, but are generally unadjusted. A time interaction tests a particular pattern of changing effect in a specified model. Agreement among them strengthens concern, while disagreement can reflect different targets, limited power, or model details.

DiagnosticWhat to look forImportant limitation
Scaled Schoenfeld residualsSystematic trend against event time; coefficient-specific or global testTrend tests depend on time transformation and can be unstable with sparse events
Log-minus-log plotGroup curves that are not approximately parallelUsually unadjusted; categorizing a continuous omics variable loses information
Time interactionEvidence against the interaction coefficient being zeroOnly detects patterns represented by the selected time function

Before interpreting any diagnostic, check that time is coded consistently, the risk sets reflect entry and censoring correctly, and the covariate specification is plausible. A nonlinear relationship entered as a straight-line effect can create residual patterns that resemble time dependence. So can influential events, measurement error, or an omitted covariate. Inspect event counts and numbers at risk over time; a smooth curve in a region with very few events is weak evidence.

For an omics analysis, many tested features create an additional problem: some apparent violations will occur by chance. If feature-wise tests are used, account for multiplicity and distinguish exploratory findings from confirmatory ones. When features were selected using outcomes, diagnostics on the same data can also be optimistic; independent validation or an appropriately designed resampling procedure is needed for claims about generalization.

The main alternative to these diagnostics is to report a single Cox coefficient and proceed under PH. That is simpler and may be reasonable when checks show no meaningful concern, but it does not test the assumption. Conversely, replacing the model immediately after one flagged p-value can overreact to noise. Diagnostics are not model-selection rules: use plots, tests, subject-matter knowledge, and sensitivity analyses together. The specific modeling and interpretation of time-varying effects is a separate topic.

Common pitfall: “The PH test was not significant, so the proportional-hazards assumption is true.” This is incorrect. A nonsignificant test means the data did not provide sufficient evidence against its null under that test and its chosen time function. Limited events or an unsuitable transformation can hide departures. Report the diagnostic and its limits rather than treating a large p-value as proof.

Check Your Understanding

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

  1. In one sentence, state the PH assumption and explain what a constant hazard ratio means.
  2. At an event time, the risk set has binary covariates \(x=(0,1)\), and the fitted coefficient is zero. The subject with \(x=0\) has the event. Calculate the unscaled Schoenfeld residual.
  3. What pattern in a log-minus-log plot is expected under PH, and name one reason that a plot for an omics predictor could mislead.
  4. A time interaction has estimate 0.12 and standard error 0.10. Calculate the Wald statistic and state whether it rejects a two-sided 5% null using the approximate critical value 1.96.
  5. Why does a nonsignificant Schoenfeld time-trend test not prove PH?

Solutions. (1) PH requires each covariate coefficient to be constant over time, so the hazard ratio for a specified covariate contrast remains the same throughout follow-up. (2) With \(\hat\beta=0\), both risk-set weights are 1; the weighted mean is \((0+1)/2=0.5\), so \(r=0-0.5=-0.5\). (3) The transformed group curves should be approximately parallel. Categorizing a continuous omics predictor loses information and makes the plot depend on the cut point; unadjusted curves may also reflect other covariate differences. (4) \(z=0.12/0.10=1.20\). Since \(|1.20|<1.96\), this approximate test does not reject at 5%; that is not evidence that PH is established. (5) The test may have low power, sparse events, or a time transformation that does not capture the departure. It assesses evidence against a particular null, not proof that the assumption holds.

Clinical Biostats