Schoenfeld Residual

Schoenfeld Residual in Cox Models and Survival Evidence for HTA

Schoenfeld residuals are the usual diagnostic for asking whether a covariate in a Cox proportional hazards model, most often the treatment arm of a trial, has the same effect early and late in follow-up. At every death or other event, the residual compares the covariate value of the patient who had the event with the value the fitted model expected among everyone still at risk, so a drift in the residuals over time points to a hazard ratio that changes. Health economists need this check because many cost-effectiveness models apply one hazard ratio to a baseline survival curve, which is valid only under the proportional hazards assumption, and the NICE manual (PMG36, section 4.6.22) states that this assumption should always be assessed. Schoenfeld introduced the residuals in 1982, and the scaled form from Grambsch and Therneau underlies the proportional hazards tests in R and Stata.

How each residual is built from a risk set

The Cox model estimates its coefficients from risk sets, the patients still under observation and event-free just before each event time. At an event, the fitted model implies an expected covariate value for the patient who fails, a weighted average over the risk set in which each patient counts in proportion to their risk score. The residual is the gap between the observed value and that expectation.

$$r_{ik} = x_{ik} - \bar{x}_{k}(t_i), \qquad \bar{x}_{k}(t_i) = \frac{\sum_{j \in R(t_i)} x_{jk} \exp(x_j^{\mathsf T}\hat{\beta})}{\sum_{j \in R(t_i)} \exp(x_j^{\mathsf T}\hat{\beta})}$$

where $r_{ik}$ is the residual for covariate $k$ at the event of patient $i$, $t_i$ is that event time, $x_{ik}$ is the patient's value of covariate $k$, $\bar{x}_{k}(t_i)$ is the risk-score weighted mean of covariate $k$ at $t_i$, $R(t_i)$ is the risk set, $x_j$ is the covariate vector of patient $j$ in it, $\hat{\beta}$ is the vector of estimated coefficients and $\exp(x_j^{\mathsf T}\hat{\beta})$ is patient $j$'s risk score.

Only patients who fail contribute a residual, the Stata manual notes, though censored patients still sit in the risk sets behind each expectation. The residuals are each event's contribution to the first derivative of the log partial likelihood, which is zero at the fitted coefficients, so for each covariate they sum to zero. With tied event times the expectation depends on the tie method, and the manual gives Peto-Breslow and Efron versions.

Schoenfeld's 1982 paper proposed plotting the residuals against time to test the proportional hazards assumption, and using histograms of them to examine fit and detect outlying covariate values.

Scaled residuals and a changing log hazard ratio

An unscaled residual for a binary treatment lies between minus 1 and 1 and cannot be read as a hazard ratio. Grambsch and Therneau framed non-proportional hazards as a Cox model with a time-varying coefficient, so that for one covariate the hazard is multiplied by $\exp\{\beta(t) x\}$, and gave a treatment effect that decreases with time as a common example. They showed that $\beta(t)$ can be seen directly by smoothing a suitable residual plot.

$$s^{*}_{i} = d\, \widehat{V}\, r_{i}, \qquad \mathrm{E}\left(s^{*}_{i}\right) + \hat{\beta} \approx \beta(t_i)$$

where $r_i$ is the vector of Schoenfeld residuals at the $i$-th event, $s^{*}_{i}$ is the vector of scaled Schoenfeld residuals, $d$ is the total number of events, $\widehat{V}$ is the estimated variance matrix of $\hat{\beta}$, and $\beta(t_i)$ is the coefficient that would apply at time $t_i$ if the effect were free to vary.

A plot of $\hat{\beta} + s^{*}_{i}$ against time is therefore a rough picture of the log hazard ratio over follow-up. Stata stores the scaled residuals with $\hat{\beta}$ already added and describes them as centred at $\hat{\beta}$, with zero slope against functions of time when proportional hazards holds. For a treatment with $\hat{\beta} < 0$, an upward trend means the log hazard ratio is moving towards zero, the pattern of a fading benefit.

Tests for a trend against time

A formal test asks whether the scaled residuals drift with time. The Stata manual notes that many popular tests for proportional hazards are in fact tests of non-zero slope in a regression of the scaled residuals on time, equivalent to testing that the log hazard ratio is constant. For one covariate the statistic is

$$\chi^{2} = \frac{\left[\sum_{i=1}^{d} \{g(t_i) - \bar{g}\}\, s^{*}_{i}\right]^{2}}{d\, \widehat{\mathrm{Var}}(\hat{\beta}) \sum_{i=1}^{d} \{g(t_i) - \bar{g}\}^{2}}$$

where $g(t_i)$ is the chosen function of the $i$-th event time, $\bar{g}$ is its mean over the $d$ events, $s^{*}_{i}$ is the scaled residual for the covariate, and $\widehat{\mathrm{Var}}(\hat{\beta})$ is the estimated variance of its coefficient. The statistic is compared with a chi-squared distribution on 1 degree of freedom, and a global test across $m$ covariates has $m$ degrees of freedom.

The choice of $g$ is part of the test. Stata's estat phtest uses analysis time by default, with log time, one minus the Kaplan-Meier estimate and ranks as options, while the R function cox.zph defaults to a Kaplan-Meier based transform and, from survival version 3.0, computes the actual score test rather than a fast approximation. The two can report different p-values on the same data, so the transform and software should be stated.

Bradburn and colleagues cite a simulation study in which the scaled Schoenfeld residuals test was among the most powerful checks of proportionality, and advise that a trend calls for a time-dependent Cox model or a model without proportional hazards.

Worked example with a fading treatment effect

All numbers in this example are illustrative. Ten trial patients are followed for overall survival, five on a new treatment ($x = 1$) and five on the comparator ($x = 0$). Comparator patients die at months 1, 2, 4 and 6 and one is censored at month 10, while treated patients die at months 3, 7, 8 and 9 and one is censored at month 12.

1. Fit the Cox model. Maximising the partial likelihood gives $\hat{\beta} = -0.609$, a hazard ratio of $\exp(-0.609) = 0.544$. The information, the sum over the eight risk sets of the weighted variance of $x$, is 1.913, so $\widehat{\mathrm{Var}}(\hat{\beta}) = 1/1.913 = 0.523$.

2. Expected value and residual at each event. With $n_0$ comparator and $n_1$ treated patients at risk, the expected value of $x$ is $\hat{\psi} n_1 / (\hat{\psi} n_1 + n_0)$, where $\hat{\psi} = 0.544$ is the fitted hazard ratio. At month 1 the five treated patients contribute 0.544 × 5 = 2.720 and the five comparator patients 5, so the expected value is 2.720 / 7.720 = 0.352 and the residual for the comparator death is minus 0.352.

3. Scale. Since $d\,\widehat{\mathrm{Var}}(\hat{\beta}) = 8 \times 0.52285 = 4.1828$, the plotted value at month 1 is $-0.6090 + 4.1828 \times (-0.35229) = -2.083$.

4. Check the sum. The negative residuals add to minus 1.698 and the positive ones to 1.699, zero apart from rounding.

5. Look for a trend. The mean of $\hat{\beta} + s^{*}$ over the events in months 1 to 4 is minus 1.292, a hazard ratio of about $\exp(-1.292) = 0.27$, and over months 6 to 9 it is 0.074, about $\exp(0.074) = 1.08$.

6. Test. With $g(t) = t$ the event months have mean 5, and the squared deviations sum to 60. The weighted sum of residuals is $\sum (t_i-5)\, r_i = 5.159$, and since $s^{*}_i = 4.182\, r_i$ the statistic is 4.182 × 5.159 × 5.159 / 60 = 1.86, with $p = 0.17$. On log time it is 1.56, with $p = 0.21$.

7. Interpretation. The plot suggests an early benefit that fades over the later months, yet with eight events the test is far from significant. One hazard ratio of 0.54 applied throughout would understate the early benefit and overstate the later one, and with so few events the plot, read with clinical reasoning, says more than the p-value.

Schoenfeld residuals in NICE appraisals

PMG36 does not name Schoenfeld residuals. Section 4.6.22 prefers log-cumulative hazard plots as advised in NICE DSU Technical Support Document 14 (TSD14), inspection of hazard plots or hazard ratios over time, and interpretation of tests reported in the original trial publications, where a residual-based test would sit. TSD14 names tests of Cox-Snell, martingale or Schoenfeld residuals as possibly useful, and notes that with patient-level data the assumption need not be relied on.

In the NICE appraisal of abemaciclib with endocrine therapy for high-risk early breast cancer (TA810), the evidence review group reports that the company checked proportional hazards with a log-cumulative plot, the Grambsch-Therneau test and Schoenfeld residual plots, based on TSD14. Finding no violation for invasive disease-free survival, the company fitted one model with a hazard ratio adjustment, while time to treatment discontinuation was modelled separately for each arm. The base case still set a fixed treatment effect lasting 8 years, followed by waning, an assumption resting on other analyses and the literature rather than on the checks.

Latimer and colleagues re-examined the 12-month data cut of the CheckMate 649 trial of nivolumab plus chemotherapy, the main evidence in NICE TA857. The data cut's Schoenfeld residual and hazard plots and the Grambsch-Therneau test were not definitive, so the authors turned to external evidence that immunotherapies have delayed and durable effects and used independently fitted survival curves. PMG36 section 4.6.23 makes the same distinction, allowing hazard ratios to be pooled only if proportional hazards holds within the trial and is clinically plausible during extrapolation.

Assumptions and choices that change the result

The test uses one pooled variance matrix, which the Stata manual says assumes homogeneity of variance across risk sets, and it recommends checking each stratum separately in a stratified Cox analysis. The R documentation notes that random effects terms, such as a frailty model term, are not checked and are treated as a fixed offset. The residuals exist for any covariate, continuous ones included, whereas Stata's log-log plot needs a nominal or ordinal covariate.

How Schoenfeld residuals differ from other Cox model residuals

Each Cox model residual answers a different question about the fit. The Stata manual lists the traditional uses.

A time-varying covariate, whose value changes during follow-up, differs from the time-varying coefficient the test looks for.

Common misreadings

Most errors treat the test as a verdict on the economic model, when it describes only the trial data. Four recur.

1. Reading a non-significant test as proof. With few events a clear trend can give a large p-value.
2. Treating a significant test as the end of Cox modelling. Bradburn and colleagues kept a Cox model in their ovarian cancer example, adding interactions with log time for the covariates whose effects drifted.
3. Carrying the within-trial result into the extrapolated period. PMG36 asks separately for clinical plausibility during extrapolation, and waning assumptions break proportionality by design.
4. Reading the unscaled residual plot as a hazard ratio curve. Only the scaled residuals plus $\hat{\beta}$ estimate the log hazard ratio over time.

Sources

• Bradburn MJ, Clark TG, Love SB, Altman DG. Survival analysis part III: multivariate data analysis, choosing a model and assessing its adequacy and fit. British Journal of Cancer. 2003;89(4):605-611. https://doi.org/10.1038/sj.bjc.6601120
• Grambsch PM, Therneau TM. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika. 1994;81(3):515-526. https://doi.org/10.1093/biomet/81.3.515
• Latimer NR. NICE DSU Technical Support Document 14: Survival analysis for economic evaluations alongside clinical trials, extrapolation with patient-level data. Decision Support Unit; June 2011 (last updated March 2013). https://www.sheffield.ac.uk/nice-dsu/tsds/survival-analysis
• Latimer NR, Taylor K, Hatswell AJ, Ho S, Okorogheye G, Chen C, et al. An evaluation of an algorithm for the selection of flexible survival models for cancer immunotherapies: pass or fail? PharmacoEconomics. 2024;42(12):1395-1412. https://doi.org/10.1007/s40273-024-01429-0
• National Institute for Health and Care Excellence. NICE technology appraisal and highly specialised technologies guidance: the manual (PMG36), sections 4.6.22 and 4.6.23. London: NICE; 2022 (last updated 13 October 2026). https://www.nice.org.uk/process/pmg36/chapter/economic-evaluation-2
• Orozco Leal G, Armstrong N, Kernohan A, Ahmadu C, Coughlan D, McDermott K, et al. Abemaciclib in combination with endocrine therapy for adjuvant treatment of hormone receptor-positive, HER2-negative, node-positive early breast cancer: an evidence review group perspective of a NICE single technology appraisal. PharmacoEconomics. 2023;41(7):741-750. https://doi.org/10.1007/s40273-023-01259-6
• Schoenfeld D. Partial residuals for the proportional hazards regression model. Biometrika. 1982;69(1):239-241. https://doi.org/10.1093/biomet/69.1.239
• StataCorp. PH plots (right-censored). Stata 19 Survival Analysis Reference Manual (accessed 28 October 2026). https://www.stata.com/manuals/stphplotsright-censored.pdf
• StataCorp. stcox postestimation. Stata 19 Survival Analysis Reference Manual (accessed 28 October 2026). https://www.stata.com/manuals/ststcoxpostestimation.pdf
• Therneau TM. cox.zph: Test the proportional hazards assumption of a Cox regression. survival: Survival Analysis, R package documentation (accessed 28 October 2026). https://rdrr.io/cran/survival/man/cox.zph.html