Checking a survival model

Diagnostics for survival models: the proportional hazards test, martingale residuals for functional form, and Cox-Snell residuals. None proves the model right.

Corrected 7 October 2026: the Cox-Snell check said it plotted the Nelson-Aalen cumulative hazard of the residuals but computed minus the log of their Kaplan-Meier curve; the code now uses the Nelson-Aalen estimate, which moves the two largest departures from the diagonal very slightly and leaves the comparison unchanged. The caption of the Cox-Snell figure said the Weibull drifts away from the line; it now says the log-logistic fit tracks the line a little more closely than the Weibull.

Corrected 27 September 2026: the residual figure’s caption, its text description and the paragraph on the test result called the fitted line for the constant covariate flat; it rises gently, and the text now says it shows a much weaker trend that the test does not flag.

The last four posts fitted survival models: the Cox model, parametric and accelerated failure time models, the cumulative incidence function for competing risks, and frailty and recurrent-event models. Each rests on assumptions, and this post is about checking them: does the hazard ratio really stay constant over time, did you get the shape of a continuous covariate right, does the parametric distribution fit?

There are good, specific diagnostics for each of those questions. There is also a hard ceiling on what any of them can tell you, and the honest reading of a clean diagnostic is “no evidence against”, never “assumption verified”. This is the recurring lesson of every checking post on the blog, and it is at its sharpest here.

library(ggplot2)
library(survival)
library(grid)

FOREST <- "#275139"; GOLD <- "#c9b458"; RUST <- "#b5534e"; SAGE <- "#93a87f"
INK <- "#16241d"; BODY <- "#2c3a31"; LINE <- "#dad9ca"

theme_te <- function() {
  theme_minimal(base_size = 12) +
    theme(plot.title = element_text(colour = INK, face = "bold", size = 12),
          plot.subtitle = element_text(colour = BODY, size = 10),
          axis.title = element_text(colour = BODY),
          axis.text = element_text(colour = BODY),
          legend.position = "none",
          panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = LINE, linewidth = 0.3),
          plot.background = element_rect(fill = "#f5f4ee", colour = NA),
          panel.background = element_rect(fill = "#f5f4ee", colour = NA))
}
two_panel <- function(pa, pb) {
  grid.newpage(); pushViewport(viewport(layout = grid.layout(1, 2)))
  print(pa, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
  print(pb, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
}

Does the hazard ratio hold over time

A Cox model claims one covariate multiplies the hazard by the same factor at every age. When the true effect fades or grows over time, that claim is false, and the standard tool to catch it is the scaled Schoenfeld residual and its test, cox.zph. A residual that trends against time is the signature of a time-varying effect.

To see it work, we build one covariate whose effect genuinely converges over time (data from a log-logistic accelerated failure time model, which gives non-proportional hazards) and one whose effect is genuinely constant (data from a Weibull proportional hazards model).

set.seed(4260)
n <- 400
x_np <- rbinom(n, 1, 0.5)
t_np <- exp(log(2) - 0.6 * x_np + 0.6 * qlogis(runif(n)))
o_np <- pmin(t_np, 3.5); d_np <- as.integer(t_np <= 3.5)
zph_np <- cox.zph(coxph(Surv(o_np, d_np) ~ x_np))

x_p <- rbinom(n, 1, 0.5); a <- 1.4; lam0 <- 2
t_p <- lam0 * exp(-log(2) * x_p / a) * (-log(runif(n)))^(1 / a)
o_p <- pmin(t_p, 4); d_p <- as.integer(t_p <= 4)
zph_p <- cox.zph(coxph(Surv(o_p, d_p) ~ x_p))

p_np <- zph_np$table["x_np", "p"]
p_p <- zph_p$table["x_p", "p"]

The test flags the converging covariate, with a proportional hazards p-value of 0.016, and clears the constant one, at 0.31. The residual plot below tells the same story visually: a downward slope for the non-proportional covariate, a much weaker upward trend for the proportional one, too weak for the test to flag. When cox.zph fires, the usual fixes are stratifying on the offending variable or letting its coefficient vary with time.

Did you get the shape right

The second question is functional form: when a continuous covariate enters the hazard, a Cox model assumes it does so linearly on the log scale. Martingale residuals from a model with no covariates, plotted against the covariate, trace out its true shape.

set.seed(4261)
nb <- 500; z <- runif(nb, 0, 4)
t_b <- -log(runif(nb)) / (0.4 * exp(0.9 * (z - 2)^2))
o_b <- pmin(t_b, 3); d_b <- as.integer(t_b <= 3)
mart_null <- residuals(coxph(Surv(o_b, d_b) ~ 1), type = "martingale")
sm <- lowess(z, mart_null, f = 0.5)
z_min <- sm$x[which.min(sm$y)]
aic_lin <- AIC(coxph(Surv(o_b, d_b) ~ z))
aic_poly <- AIC(coxph(Surv(o_b, d_b) ~ poly(z, 2)))

Here the covariate acts through a U-shape, low risk in the middle and high at both ends, and the smoother through the residuals recovers it, bottoming out at 2.11 against a true minimum at 2. Fitting the covariate as linear is badly wrong; adding a quadratic term drops the AIC from 4854.5 to 4552.8, a fall of 302.

Figure 1: Left: scaled Schoenfeld residuals against time. The converging covariate trends downward (its effect fades); the constant one shows only a much weaker upward trend, which the test does not flag. Right: martingale residuals against a continuous covariate, with a smoother recovering the true U-shaped functional form.

Does the distribution fit

For a parametric model, the question is whether the assumed distribution matches the data. Cox-Snell residuals answer it: if the model fits, the residuals behave like a sample from a unit-exponential distribution, so their Nelson-Aalen cumulative hazard, plotted against the residual values, should follow the forty-five degree line. We reuse the log-logistic data from the parametric survival post and fit both the wrong distribution (Weibull) and the right one (log-logistic).

set.seed(4259)
nc <- 500; mu_c <- log(2.2); s_c <- 0.55
t_ll <- exp(mu_c + s_c * qlogis(runif(nc)))
cc <- 3.5; o_c <- pmin(t_ll, cc); d_c <- as.integer(t_ll <= cc)
fit_w  <- survreg(Surv(o_c, d_c) ~ 1, dist = "weibull")
fit_ll <- survreg(Surv(o_c, d_c) ~ 1, dist = "loglogistic")
cs_w  <- { a <- 1 / fit_w$scale;  l <- exp(coef(fit_w)[1]); (o_c / l)^a }
cs_ll <- { b <- 1 / fit_ll$scale; al <- exp(coef(fit_ll)[1]); log(1 + (o_c / al)^b) }
na_dev <- function(cs, d) {
  sf <- survfit(Surv(cs, d) ~ 1); H <- sf$cumhaz
  max(abs(H - sf$time))
}
dev_w <- na_dev(cs_w, d_c); dev_ll <- na_dev(cs_ll, d_c)

The correct model tracks the line more closely: the largest departure from forty-five degrees is 0.059 for the log-logistic against 0.080 for the Weibull. The gap is real but modest, and that caveat matters: with censoring the Cox-Snell check is a fairly blunt instrument; it will catch a badly wrong distribution but not split hairs between two plausible ones.

The two things a proportional hazards test cannot do

Now the ceiling. The proportional hazards test is neither a proof nor a disproof, and two short simulations show why.

First, a clean test does not prove proportionality; it may simply lack power. We generate genuinely non-proportional data and ask how often cox.zph detects it at various sample sizes.

power_zph <- function(nn, B = 200, seed = 7) {
  set.seed(seed); k <- 0
  for (b in 1:B) {
    x <- rbinom(nn, 1, 0.5); tt <- exp(log(2) - 0.6 * x + 0.6 * qlogis(runif(nn)))
    oo <- pmin(tt, 3.5); dd <- as.integer(tt <= 3.5)
    if (cox.zph(coxph(Surv(oo, dd) ~ x))$table["x", "p"] < 0.05) k <- k + 1
  }
  k / B
}
ns <- c(80, 200, 500, 1000)
pw <- sapply(ns, power_zph)

At 200 animals, about 160 of which die before the study ends, the test finds the violation only 33 per cent of the time, and even at 500 it misses roughly a third. A non-significant cox.zph on a modest data set is very weak evidence that hazards are proportional.

Second, a failed test does not prove the effect is time-varying; an omitted covariate can manufacture the appearance of non-proportionality. We generate data where two covariates both have strictly constant effects, then fit a model leaving one out.

set.seed(4262); B <- 250; flag_omit <- flag_full <- 0
for (b in 1:B) {
  nd <- 350; x1 <- rbinom(nd, 1, 0.5); x2 <- rnorm(nd)
  td <- lam0 * exp(-(log(2) * x1 + 1.5 * x2) / a) * (-log(runif(nd)))^(1 / a)
  od <- pmin(td, 4); dd <- as.integer(td <= 4)
  if (cox.zph(coxph(Surv(od, dd) ~ x1))$table["x1", "p"] < 0.05) flag_omit <- flag_omit + 1
  if (cox.zph(coxph(Surv(od, dd) ~ x1 + x2))$table["x1", "p"] < 0.05) flag_full <- flag_full + 1
}
rate_omit <- flag_omit / B; rate_full <- flag_full / B

With the strong covariate omitted, cox.zph flags the first covariate as non-proportional 16 per cent of the time, far above the nominal five, even though its true effect never changes. Put the covariate back and the flag rate drops to 5 per cent, where it belongs. The apparent time-varying effect was an artefact of unmodelled heterogeneity, not a real feature of the hazard.

Figure 2: Left: Cox-Snell residual cumulative hazard against the forty-five degree line; the log-logistic fit tracks it a little more closely than the Weibull. Right: the power of the proportional hazards test to detect a genuine violation, low at small samples, so a clean test is weak evidence of proportionality.

What to take away

Survival models come with sharp, specific diagnostics, and you should run all of them. Use scaled Schoenfeld residuals and cox.zph for proportional hazards, martingale residuals for the functional form of a continuous covariate, and Cox-Snell residuals for a parametric distribution; for the frailty and recurrent-event models, the gap between the model-based and cluster-robust standard error is itself a check for unmodelled correlation, and for competing risks you would test the Fine-Gray proportional subdistribution hazards the same way. But read every result with the ceiling in mind. A clean proportional hazards test can simply lack power, and a failed one can be an omitted covariate in disguise. The diagnostics test the model you fit, not the question you care about, and the deepest assumptions, that competing events are handled correctly, that no important covariate is missing, that the sampled animals represent the population, are not checkable from the data at all. Run the checks, act on what they flag, and treat silence as the absence of evidence rather than evidence of correctness.

References

Cox DR, Snell EJ 1968 Journal of the Royal Statistical Society: Series B 30:248-265 (10.1111/j.2517-6161.1968.tb00724.x)

Grambsch PM, Therneau TM 1994 Biometrika 81:515-526 (10.1093/biomet/81.3.515)

Schoenfeld D 1982 Biometrika 69:239-241 (10.1093/biomet/69.1.239)

Therneau TM, Grambsch PM 2000 Modeling Survival Data: Extending the Cox Model. Springer. ISBN 978-0-387-98784-2

Therneau TM, Grambsch PM, Fleming TR 1990 Biometrika 77:147-160 (10.1093/biomet/77.1.147)

Related tutorials

• Latency tests when some animals never respond
• Parametric survival and the AFT model
• Frailty and recurrent event models
• Time-varying covariates in a survival model

Reproducibility

• Code last run: 2026-10-15
• R version 4.5.3 (2026-03-11), aarch64-apple-darwin20
• Packages: ggplot2 4.0.3, survival 3.8-12