24  Cox Model Diagnostics

24.1 Overview

After fitting a Cox model, we should check two key assumptions:

  1. Proportional hazards: the hazard ratio between any two individuals is constant over time.
  2. Linearity: continuous explanatory variables enter the log-hazard linearly.

Several diagnostic tools are available.

24.2 Graphical check: log-cumulative hazard plots

The simplest PH check for a categorical explanatory variable is the log-cumulative hazard plot (introduced in Chapter 21): plot \(\log\widehat{H}(t)\) from the KM estimator against \(\log t\) for each group. Under PH, the curves should be parallel.

Code 
library(survival)

km_g <- survfit(Surv(t, delta) ~ group, data = mydata)
plot(km_g, fun = "cloglog",          # log H(t) vs log(t)
     col = c("blue", "red"), lwd = 2,
     xlab = "log(t)", ylab = "log H(t)")

24.3 Schoenfeld residuals

Definition: Schoenfeld residuals

At each event time \(t_{(j)}\), the Schoenfeld residual for explanatory variable \(k\) is the observed value of that variable minus the risk-set-weighted expected value: \[ r_{jk} = x_{(j)k} - \frac{\sum_{\ell \in \mathcal{R}(t_{(j)})} x_{\ell k} \exp(\hat{\bbeta}^T\mathbf{x}_\ell)} {\sum_{\ell \in \mathcal{R}(t_{(j)})} \exp(\hat{\bbeta}^T\mathbf{x}_\ell)}. \]

Under the PH assumption, the Schoenfeld residuals should have mean zero and no trend with time. A plot of \(r_{jk}\) against event time \(t_{(j)}\) that shows a systematic trend (increasing, decreasing, or non-linear) suggests that the PH assumption is violated for explanatory variable \(k\).

The cox.zph() function in R tests formally for a time trend in each explanatory variable’s Schoenfeld residuals.

Code 
library(survival)

cox_fit <- coxph(Surv(t, delta) ~ x1 + x2, data = mydata,
                 ties = "breslow")

# Schoenfeld residual test for PH
ph_test <- cox.zph(cox_fit)
print(ph_test)   # p-values for each explanatory variable + global test

# Plot Schoenfeld residuals
plot(ph_test)    # one plot per explanatory variable + global

A small \(p\)-value from cox.zph() suggests the PH assumption is violated for that explanatory variable.

24.4 Martingale residuals

Definition: Martingale residuals

\[ M_i = \delta_i - \widehat{H}_0(t_i)\exp(\hat{\bbeta}^T\mathbf{x}_i), \] where \(\widehat{H}_0(t)\) is the Breslow estimate of the baseline cumulative hazard.

Martingale residuals lie in \((-\infty, 1]\) and sum to zero. Plotting them against a continuous explanatory variable that was assumed to enter linearly can reveal non-linearity: a non-linear trend in the residuals suggests that a transformation of the explanatory variable (or a smooth term) would improve the fit.

Code 
# Martingale residuals from a null (intercept-only) model
null_fit <- coxph(Surv(t, delta) ~ 1, data = mydata, ties = "breslow")
mart_res  <- residuals(null_fit, type = "martingale")

# Plot vs continuous explanatory variable to check linearity
plot(mydata$age, mart_res,
     xlab = "Age", ylab = "Martingale residual")
lines(lowess(mydata$age, mart_res), col = "red", lwd = 2)
abline(h = 0, lty = 2)

24.5 Deviance residuals

Deviance residuals are a transformation of the martingale residuals designed to be more symmetric around zero for a well-fitting model: \[ d_i = \text{sign}(M_i)\sqrt{-2\bigl(M_i + \delta_i\log(\delta_i - M_i)\bigr)}. \]

Large absolute deviance residuals (say, \(|d_i| > 3\)) may indicate influential observations or outliers.

Code 
dev_res <- residuals(cox_fit, type = "deviance")
plot(predict(cox_fit), dev_res,
     xlab = "Linear predictor", ylab = "Deviance residual")
abline(h = c(-2, 0, 2), lty = c(2, 1, 2))

24.6 Complete diagnostics workflow

Code 
library(survival)

h <- read.table("https://www.richardpmann.com/MATH3701/Datasets/hypernephroma.dat", header = TRUE)
h$age         <- as.factor(h$age)
levels(h$age)  <- c("<60", "60-70", ">70")
h$nephrectomy <- as.factor(h$nephrectomy)
levels(h$nephrectomy) <- c("no", "yes")

# Fit final model
fit_final <- coxph(Surv(time, status) ~ age + nephrectomy,
                   data = h, ties = "breslow")

# 1. PH test
ph_test <- cox.zph(fit_final)
print(ph_test)
plot(ph_test)

# 2. Martingale residuals (all explanatory variables are categorical here)
mart <- residuals(fit_final, type = "martingale")
plot(fitted(fit_final), mart,
     xlab = "Fitted log-hazard", ylab = "Martingale residual")

# 3. Deviance residuals
dev <- residuals(fit_final, type = "deviance")
plot(seq_along(dev), dev, type = "h",
     xlab = "Observation index", ylab = "Deviance residual")
abline(h = c(-2, 2), lty = 2)

Exercises

24.1. After fitting the Cox proportional hazards model from Exercise 23.1 to the Drug A/B data, a researcher wishes to check the model assumptions. Without using R, describe the diagnostic approach for each of the following.

Think about log–log survival plots (\(\log(-\log\hat{S}(t))\) vs \(\log t\) for each group), Schoenfeld residuals plotted against time, and dfbeta or deviance residuals for influence assessment.

  1. Checking the proportional hazards assumption. What plot(s) would you produce, what would you look for, and what would departures from PH look like?

  2. Identifying influential observations. What residual-based diagnostics would you use, and what pattern in the plot would flag an observation as influential?

24.2. Using the robot surgery data from Exercise 20.1, perform a full set of Cox model diagnostics in R.

Use coxph() to fit the model, cox.zph() to test PH, plot(cox.zph(fit)) to visualise Schoenfeld residuals, and residuals(fit, type="deviance") for deviance residuals.

  1. Fit a Cox model with group as the explanatory variable. Run cox.zph() and interpret the output.

  2. Plot the Schoenfeld residuals against time. What pattern would confirm that the PH assumption is satisfied?

  3. Plot the deviance residuals against the linear predictor. Identify any potential outliers.