Cox-Snell Residual

Cox-Snell Residual for Checking Survival Curve Fit Before Extrapolation

Economic models of cancer and other chronic conditions often rest on survival curves fitted to trial data and then projected far beyond the end of follow-up. The Cox-Snell residual turns a fitted parametric survival model into one number per patient, and the plot built from those numbers shows whether the model reproduces the pattern of events actually observed. When the model is right the residuals should resemble a censored unit exponential sample, and the plot, built with the Nelson-Aalen estimator, is read through departures from its reference line. The plot belongs within NICE Decision Support Unit guidance on choosing survival curves, and it has clear limits for the extrapolated part of a curve.

From a fitted survival curve to one residual per patient

Observed minus fitted values do not work as residuals for censored, skewed survival times. Cox and Snell published a general definition of residuals in 1968 that was not limited to linear models, and Kay later used their approach to derive residuals for proportional hazards regression. For a time-to-event model the Cox-Snell residual of patient $i$ is the fitted cumulative hazard evaluated at that patient's own observed time, as follows.

$$r_i = \hat{H}(t_i \mid x_i) = -\log \hat{S}(t_i \mid x_i)$$

where $r_i$ is the Cox-Snell residual for patient $i$, $t_i$ is the observed time (an event time or a censoring time), $x_i$ are the patient's covariate values, such as treatment arm, $\hat{H}$ is the cumulative hazard implied by the fitted model and $\hat{S}$ is the corresponding fitted survival function. For a Cox proportional hazards model the same quantity is written $r_i = \hat{H}_0(t_i)\exp(\hat{\beta}' x_i)$, with $\hat{H}_0$ the estimated cumulative baseline hazard and $\hat{\beta}$ the estimated coefficients.

Each residual measures how much hazard the model says a patient had accumulated by the time that patient was last seen. It is never negative and is not symmetric around zero, so it does not behave like an ordinary residual in a regression plot.

Why a correct model yields unit exponential residuals

The reference distribution comes from a property that holds for any continuous survival time. If a random survival time $T$ has cumulative hazard $H$, the transformed value $H(T)$ follows an exponential distribution with mean 1. This follows because survival and cumulative hazard are linked by $S(t) = \exp\{-H(t)\}$, so that the following holds.

$$P\{H(T) > u\} = P\{T > H^{-1}(u)\} = \exp(-u)$$

where $u$ is any non-negative value, $H^{-1}$ is the inverse of the cumulative hazard and $\exp(-u)$ is the survival function of the unit exponential distribution. When the fitted model is correct, $\hat{H}$ is close to the true $H$, so in the absence of censoring the residuals behave like draws from that distribution.

Censoring changes how the check has to be done. A censored patient's residual is only a lower bound on the value it would have reached, so the residuals form a censored sample from the unit exponential, and a histogram, mean or quantile plot of the raw values misleads when follow-up is incomplete. The usual remedy treats the residuals as survival times, keeps each patient's original event or censoring flag, and estimates their cumulative hazard nonparametrically. The unit exponential has a constant hazard of 1, so its cumulative hazard is a straight line through the origin with slope 1, and that line is the benchmark for the plot.

Building and reading the residual plot

The plot takes four steps once a candidate model has been fitted, by maximum likelihood estimation for a parametric model or by partial likelihood for a Cox model. The flexsurv package for R returns the residuals, and Stata's estat gofplot command draws the whole plot after streg or stcox.

1. Compute the residuals. Each patient's fitted cumulative hazard at the observed time is calculated from the model, including the patient's covariates.

2. Treat the residuals as survival times. The residuals are sorted and each one keeps its original event or censoring indicator.

3. Estimate their cumulative hazard. The Nelson-Aalen estimator is the usual choice, with minus the log of the Kaplan-Meier estimator as an alternative. The estimate is as follows.

$$\hat{H}_R(r) = \sum_{j:\ r_{(j)} \le r} \frac{d_j}{n_j}$$

where $\hat{H}_R(r)$ is the estimated cumulative hazard of the residuals at value $r$, $r_{(j)}$ are the distinct residual values at which events occurred, $d_j$ is the number of events at $r_{(j)}$ and $n_j$ is the number of patients whose residual is at least $r_{(j)}$.

4. Plot against the 45 degree line. $\hat{H}_R(r)$ is plotted against $r$ together with the line of slope 1 through the origin.

A departure from the line can be read in terms of the model's error. For patients who share the same covariate values, the cumulative hazard of the residuals is the true cumulative hazard re-expressed on the scale of the fitted one, as follows.

$$H_R(x) = H\{\hat{H}^{-1}(x)\}$$

where $H_R(x)$ is the cumulative hazard of the residuals at value $x$, $H$ is the true cumulative hazard and $\hat{H}^{-1}$ is the inverse of the fitted cumulative hazard. Here a point above the line means that the model understates the cumulative hazard up to the corresponding time, and so overstates survival at that time; a point below the line means the reverse. The local slope of the plot reflects the ratio of the true to the fitted hazard at that time. A curve that bends steadily away from the line suggests the wrong hazard shape, for example a constant hazard fitted where the observed hazard rises, which the log-cumulative hazard plots recommended in NICE DSU guidance can then help to diagnose.

Worked example of an exponential curve fitted to eight patients

The figures in this example are illustrative. Eight patients in one trial arm are followed for up to 14 months; six die and two are censored. An exponential model, which assumes a constant hazard rate, is fitted. Its maximum likelihood estimate of the hazard is the number of events divided by the total time at risk, as follows.

$$\hat{\lambda} = \frac{6}{60} = 0.1$$

where $\hat{\lambda}$ is the estimated monthly hazard, 6 is the number of deaths and 60 is the sum of the eight follow-up times in months. Under this model $\hat{H}(t) = \hat{\lambda} t$, so each residual is 0.1 times the observed time. The Nelson-Aalen column adds $d_j / n_j$ at each event residual and stays flat at censored ones.

The value at the sixth residual, for example, is as follows.

$$\hat{H}_R(1.0) = \frac{1}{8} + \frac{1}{7} + \frac{1}{6} + \frac{1}{4} + \frac{1}{3} = 1.018$$

where each fraction is one death divided by the number still at risk. The first five event points lie close to the line of slope 1 (0.125 against 0.2, 0.268 against 0.3, 0.435 against 0.5, 0.685 against 0.8 and 1.018 against 1.0), slightly below it early on, which would mean the constant hazard slightly overstates early deaths. The final point jumps to 2.018 against 1.4 because only one patient is still at risk, so a single death adds a full unit; the jump is consistent with right-tail noise and says little about the model.

The residuals also sum to 6, the number of deaths, which is a property of the exponential maximum likelihood fit, because the log-likelihood $d \log \lambda - \lambda \sum t_i$ is maximised at $\hat{\lambda} = d / \sum t_i$, so $\sum r_i = \hat{\lambda} \sum t_i = d$. Under this fit the raw mean of the residuals always equals the proportion of patients with an event, here 6/8 = 0.75, so the mean carries no information on fit. With eight patients, the plot can do little more than rule out gross misfit, which is why such plots are read alongside the other checks set out below.

Where the residual plot sits in survival model selection for HTA

Survival curves set the state occupancy of a partitioned survival model, and the extrapolated tail often drives mean survival and the QALY gain. The NICE Decision Support Unit's Technical Support Document 14 (TSD 14), written by Latimer in 2011 and last updated in March 2013, sets out how fitted models should be justified when patient-level trial data are available. It names tests of residuals such as Cox-Snell, martingale or Schoenfeld residuals as possibly useful, while describing visual inspection, log-cumulative hazard plots and information criteria as particularly relevant, and it gives separate weight to external data and clinical validity for judging the extrapolated portion. Its review of past NICE appraisals found that tests of residuals were also occasionally undertaken.

Each check answers its own question, and they work together. The Akaike information criterion and the Bayesian information criterion rank candidate distributions against each other, so the best-ranked model can still fit poorly if every candidate does; the AIC page sets out how TSD 14 and the NICE manual treat these criteria. A Cox-Snell plot looks at one model at a time and asks whether its goodness of fit is adequate in absolute terms across the whole observed period.

TSD 14 states that visual inspection, log-cumulative hazard plots and AIC and BIC address internal validity only, and that when substantial extrapolation is needed the predictions should be validated by other means. A residual plot is confined to the observed data in the same way.

Choices that change what the plot shows

Several analysis choices change the look of the plot without any change in the underlying data. They matter most when a plot is used to separate noise from misfit.

The fitted parameters are themselves estimates. The documentation for Stata's Cox model notes that the unit exponential result holds only when the true parameters are used, so deviations from the line can partly reflect estimation uncertainty, especially in small samples and in the right tail, where prior events and censoring have reduced the effective sample. The estimator of the residuals' cumulative hazard is a second choice. Stata uses the Nelson-Aalen estimator by default after streg and stcox and offers minus the log of the Kaplan-Meier estimator as an option, so a report should say which was used.

The model structure matters as well. When arms are fitted separately, as TSD 14 allows when hazards are not proportional, each arm has its own residuals and its own plot. When treatment enters as a covariate, the residuals use each patient's covariate values, and Stata's estat gofplot can draw a separate curve for each group, so a misfit confined to one arm is not averaged away in a pooled plot. For flexible parametric models the residual is still the fitted cumulative hazard, and the flexsurv residual function accepts spline models as well as standard distributions.

What the plot cannot show

The limits of a Cox-Snell plot matter most where economic models lean on it. It describes only the observed follow-up. TSD 21 on flexible methods, by Rutherford and colleagues, states that a good fit within the range of the data may still lead to implausible extrapolations, and a plot of the observed residuals can say nothing about the tail of the curve beyond the trial.

The plot also gives an overall verdict without a diagnosis. A bend away from the line shows that something is wrong but not whether the cause is the distribution, a violated proportional hazards assumption or a missing covariate. It offers no formal test statistic either, so judging how close is close enough remains a matter of judgement. Each event adds one step to the plot, so the plot carries most information where events are numerous, while its right-hand end, closest to the period that matters for extrapolation, rests on the few patients still at risk and is the noisiest part.

One misreading recurs in practice, namely comparing the raw residuals with an exponential distribution without allowing for censoring, which can make a well-fitting model look poor, as the residual mean of 0.75 in the worked example shows.

Cox-Snell, martingale, deviance and Schoenfeld residuals

The Cox-Snell residual is one of a small family of survival residuals, and the martingale and deviance residuals are built from it. Each member suits a different diagnostic task.

The martingale residual is the event indicator minus the Cox-Snell residual, $M_i = d_i - r_i$, where $d_i$ is 1 for an event and 0 for censoring; it compares observed with expected events for each patient and is used to examine the functional form of a covariate. The deviance residual rescales the martingale residual so that it is symmetric about zero, which makes it the preferred choice for identifying outlying patients. The Schoenfeld residual, defined separately at each event time, is used in Cox models to check and test the proportional hazards assumption. Of the four, only the Cox-Snell plot compares the whole fitted model with a fixed reference distribution, which is why it is the residual suited to judging the overall fit of each candidate parametric curve.

Sources

• Collett D. Modelling Survival Data in Medical Research. 3rd ed. Boca Raton: CRC Press; 2015.
• Cox DR, Snell EJ. A general definition of residuals. Journal of the Royal Statistical Society: Series B (Methodological). 1968;30(2):248-265. https://doi.org/10.1111/j.2517-6161.1968.tb00724.x
• Jackson C. flexsurv: Flexible parametric survival and multi-state models. R package version 2.3.2. 2024. Help page coxsnell_flexsurvreg. https://search.r-project.org/CRAN/refmans/flexsurv/html/coxsnell_flexsurvreg.html
• Latimer N. NICE DSU Technical Support Document 14: Survival analysis for economic evaluations alongside clinical trials, extrapolation with patient-level data. Sheffield: Decision Support Unit, ScHARR, University of Sheffield; 2011, last updated March 2013. Sections 3, 3.1, 3.2, 3.5, 3.6, 4.4 and 6.1. https://doi.org/10.15131/shef.data.33745843
• Rutherford MJ, Lambert PC, Sweeting MJ, Pennington B, Crowther MJ, Abrams KR, Latimer NR. NICE DSU Technical Support Document 21: Flexible methods for survival analysis. Sheffield: Decision Support Unit, ScHARR, University of Sheffield; 2020 (updated March 2022). Section 3.2.2. https://doi.org/10.15131/shef.data.33849118
• StataCorp. Stata Survival Analysis Reference Manual, Release 19. Entries streg postestimation, stcox postestimation and estat gofplot. College Station, TX: Stata Press; 2025. https://www.stata.com/manuals/ststregpostestimation.pdf, https://www.stata.com/manuals/ststcoxpostestimation.pdf, https://www.stata.com/manuals/stestatgofplot.pdf
• Xue Y, Schifano ED. Diagnostics for the Cox model. Communications for Statistical Applications and Methods. 2017;24(6):583-604. https://doi.org/10.29220/CSAM.2017.24.6.583