Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Cox Proportional Hazards › Visualizing Cox Proportional Hazards

Michael BrenndoerferOctober 26, 20268 min read

Machine Learning from ScratchAll chapters

Part 6 of 15Survival Analysis

Chapter 50Cox Proportional Hazards › Visualizing Cox Proportional Hazards

Reading section 7 of 10 · Visualizing Cox Proportional Hazards · 8 minSection 7 of 10 · 8 min

Section 7 of 10 · 8 min read

Cox Proportional Hazards › Visualizing Cox Proportional HazardsSection view · Sections in this chapter7 / 10 ↓

Sections in this chapter10

Reading tools Definitions, links & notes

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Notes & comments

Notes & commentsLeave notes and join the discussionSelection toolbarShow comment actions when selecting text. Turn off for the browser menu only.Comment indicatorsShow small counts beside annotated passagesPassage highlightsGently shade annotated text selections

The first visualization shows what the fitted Cox model implies about survival. We plot the baseline-implied survival curves for maintained and unmaintained equipment on a common time axis. The cumulative baseline is evaluated as a right-continuous step function, with each event time contributing an increment and the value staying constant between events. The comparison holds the fitted coefficient and baseline fixed while changing only the maintenance covariate. This is a model-implied association in a constructed dataset, not a causal experiment conclusion.

The residual panels ask a separate question: whether the fixed-coefficient model leaves a pattern over event time. Individual residuals are noisy, and a fitted trend is a descriptive aid rather than a formal proportional-hazards test. We compare the small equipment example with a separate seeded generator whose two groups have nonproportional hazards by construction. The absence of a visible trend in a single small sample cannot validate the assumption.

VisualizationOut [28]
CollapseExpand

Two declining step curves on a shared time axis. The maintained curve lies at or above the unmaintained curve, with a changing gap between them. — Survival curves fitted by the Cox model for maintained and unmaintained equipment. The maintained curve stays at or above the unmaintained curve, and the vertical gap changes with time even though the fitted hazard ratio is constant. This model-implied contrast comes from a constructed teaching dataset and is descriptive, not causal.

Both fitted curves are nonincreasing and start at one before the first event; the first plot cell shows a visible title and uses the shared style. A negative fitted maintenance coefficient places its curve at or above the reference curve. The absolute survival gap changes with cumulative hazard, and it need not reach its maximum or return to zero within the observed follow-up. The model's constant hazard multiplier is therefore a different quantity from a constant survival-probability difference or ratio.

The fitted survival curves imply a survival-probability ratio that changes with time. That ratio is a different quantity from the constant fitted hazard ratio, even though both are derived from the same Cox model.

The accompanying plot contrasts the two quantities directly: the survival-probability ratio rises in steps with accumulated baseline hazard, while the hazard ratio remains flat. This compares a rate ratio with a ratio of survival probabilities; neither plotted ratio is itself a probability.

PythonIn [29]
CollapseExpandsurvival_ratio = surv_maintained / surv_unmaintained
fitted_hazard_ratio = float(np.exp(beta_hat[0]))

VisualizationOut [30]
CollapseExpand

A right-continuous rising step curve for the survival-probability ratio, compared with a flat hazard ratio line. — Survival-probability ratio between maintained and unmaintained equipment over time, compared with the constant fitted hazard ratio. The survival ratio rises in steps as the cumulative baseline hazard grows, while the hazard ratio stays fixed under the model. This model-implied contrast comes from a constructed teaching dataset and is descriptive, not causal.

The next visualization is a diagnostic plot. For each event, the Schoenfeld residual is a single number per covariate. Under proportional hazards, the true-parameter conditional residual mean is zero. Estimated residuals remain noisy; a trend is a reason to investigate, not proof of a changing effect.

VisualizationOut [31]
CollapseExpand

A scatter plot of Schoenfeld residuals against event time, with a dashed horizontal line at zero and a dotted straight-line trend. The residuals scatter around zero. — Schoenfeld residuals for the maintenance covariate plotted against event time. The residuals scatter around zero with no strong linear trend, which is consistent with proportional hazards. With only thirteen event residuals, this descriptive plot leaves substantial uncertainty about time-varying effects; an apparently flat trend does not validate proportional hazards.

In our equipment data the residuals are sparse (thirteen events) and the trend line is nearly flat, which is consistent with proportional hazards. Thirteen event residuals give only a sparse view of changes over time. A nearly flat fitted line can miss departures that occur outside the observed event times or have little linear association with time. We will discuss how to interpret residuals rigorously, and how to test the proportional-hazards assumption, in the Practical Implications section.

Proportional-hazards diagnostics and limitationsLink Copied

The Schoenfeld residual is the building block of the standard proportional-hazards diagnostic. Formally, for the $p$ -th covariate and event time $t_{(k)}$ where individual $i_k$ fails, we defined

r_{i_k, p} = x_{i_k, p} - \bar{x}_p(t_{(k)}), \qquad \bar{x}_p(t_{(k)}) = \frac{\sum_{j \in R(t_{(k)})} x_{j, p} \exp(\hat\beta^\top x_j)}{\sum_{j \in R(t_{(k)})} \exp(\hat\beta^\top x_j)}.

At the true coefficient under the proportional-hazards model, conditional on a risk set and one failure, the event's covariate has the weighted risk-set mean. Its unscaled Schoenfeld residual therefore has conditional mean zero. Residuals evaluated at an estimated coefficient are noisy, so a visible trend suggests investigation rather than proving a changing effect. Grambsch and Therneau (1994) develop tests based on scaled Schoenfeld residuals and transformed time. This chapter plots descriptive trends; it does not execute that formal test.

This diagnostic has two limitations. First, sensitivity depends on event count and timing, covariate variation in the risk sets, and the departure being tested. Sparse events can leave little information about time-varying effects, while a precisely estimated departure can be statistically detectable even when its practical effect is small. Interpret effect size and uncertainty alongside the p-value. Second, a test based on a chosen transformation of time can have low power against departures with little association with that transformation. Nonmonotone patterns can therefore be missed, which is why residual plots are useful alongside a formal test.

To make the failure mode concrete, we simulate a dataset in which the hazard ratio is explicitly increasing over time. The covariate is binary; the $x = 0$ group follows an exponential with constant hazard, and the $x = 1$ group follows a Weibull with shape 2, which has an increasing hazard. The hazard ratio between the groups therefore grows over time, violating proportional hazards by construction.

PythonIn [32]
CollapseExpandrng = np.random.default_rng(42)
n_sim = 300
x_sim = rng.integers(0, 2, size=n_sim).astype(float)
t_event_sim = np.where(
    x_sim == 0, rng.exponential(1.0, n_sim), 0.7 * rng.weibull(2.0, n_sim)
)
censor_sim = rng.uniform(0.0, 2.5, n_sim)
time_sim = np.minimum(t_event_sim, censor_sim)
event_sim = (t_event_sim <= censor_sim).astype(int)

X_sim = x_sim.reshape(-1, 1)
fit_sim = optimize.minimize(
    lambda b: -cox_loglike_untied(b, X_sim, time_sim, event_sim),
    x0=np.zeros(1),
    jac=lambda b: -cox_gradient(b, X_sim, time_sim, event_sim),
    method="BFGS",
    options={"gtol": 1e-7},
)
assert (
    fit_sim.success
    and np.max(np.abs(cox_gradient(fit_sim.x, X_sim, time_sim, event_sim)))
    < 1e-5
)
time_sim_events, resid_sim = schoenfeld_residuals(
    fit_sim.x, X_sim, time_sim, event_sim
)
trend_sim = np.polyfit(time_sim_events, resid_sim[:, 0], 1)
trend_sim_times = np.sort(time_sim_events)
trend_sim_values = np.polyval(trend_sim, trend_sim_times)

The following plot displays the individual residuals and their fitted trend for this known nonproportional generator.

VisualizationOut [33]
CollapseExpand

Two bands of binary-covariate Schoenfeld residuals against event time, with a dashed fitted trend rising from negative to positive. — Schoenfeld residuals from a simulated dataset whose hazard ratio increases with time by construction. Individual binary-covariate residuals form two bands, while the fitted descriptive trend rises from negative to positive. A similar trend in real data would prompt investigation rather than establish a violation on its own.

In this constructed example, the fitted residual trend rises over event time, while individual binary-covariate residuals form two bands. A similar trend in a real dataset is a reason to investigate proportional hazards and model specification. If a changing effect is supported, possible responses include a time-varying coefficient, stratification for a categorical covariate, or a spline interaction between the covariate and time.

Comments

Add to the conversation

A question, a correction, or an example from your own work.

No comments yet. Be the first to share your thoughts!

Proportional-hazards diagnostics and limitationsLink Copied

Comments

Add to the conversation

Stay up to date