I do a survival regression on some time-to-event data (vehicle breakdowns) with some covariates (essentially the age of the vehicle and some boolean variables for vehicle type). I must admit that I am not so sure about the quality of the data and whether they result from multiple processes instead of one.

For the implementation, I use the Python package lifelines and the CoxPH fitter in particular.

This is the distribution of survival times in my data (for both censored and non-censored points!):

Distribution of survival times

Here is the model summary:

<lifelines.CoxPHFitter: fitted with 2081 total observations, 354 right-censored observations>
             duration col = 'duration'
                event col = 'observed'
              cluster col = 'vehicle'
          robust variance = True
      baseline estimation = breslow
   number of observations = 2081
number of events observed = 1727
   partial log-likelihood = -11593.48
         time fit was run = 2023-03-18 06:32:51 UTC

                               coef  exp(coef)   se(coef)   coef lower 95%   coef upper 95%  exp(coef) lower 95%  exp(coef) upper 95%
year_of_manufacture            0.05       1.06       0.01             0.03             0.08                 1.03                 1.08
N7411   0.90       2.46       0.18             0.54             1.26                 1.72                 3.52
N7412   1.25       3.50       0.18             0.91             1.60                 2.48                 4.94
N7413   1.13       3.11       0.17             0.80             1.46                 2.24                 4.31
N7460  -0.47       0.63       0.33            -1.12             0.19                 0.33                 1.21
N7440   0.76       2.14       0.20             0.38             1.14                 1.46                 3.13

                               cmp to     z      p   -log2(p)
year_of_manufacture              0.00  4.15 <0.005      14.86
N7411     0.00  4.95 <0.005      20.33
N7412     0.00  7.11 <0.005      39.68
N7413     0.00  6.76 <0.005      36.06
N7460     0.00 -1.39   0.16       2.61
N7440     0.00  3.89 <0.005      13.28
Concordance = 0.56
Partial AIC = 23198.95
log-likelihood ratio test = 163.41 on 6 df
-log2(p) of ll-ratio test = 106.14

This is how the modeled baseline hazard looks:

baseline hazard

I now have some questions:

  1. When I plot the deviance residuals of the model, a very clear pattern is visible: Deviance residuals of CoxPH model

(Red dots indicate censored values, blue dots indicate events).

I assume that this pattern indicates that the model is missing some crucial feature which allows it to explain the visible non-linear relation between explanatory and target variable. And it seems to me that the model is underestimating hazard for low durations and overestimating it for large durations (please correct me if I am wrong).

  • But why can't the model adapt the hazard function to this slope, as the baseline hazard is non-parametric and thus not bound to any specific distribution of survival times?
  • Besides adding other covariables: What could I do to improve the model?
  1. Is my interpretation correct that the model says all features except N7460 do have a significant effect on the outcome and that a concordance value of 0.56 tells us that the model does slightly better than plainly guessing average values?

1 Answer 1


The plots that you show aren't very useful for evaluating a Cox model. Intuition from models like ordinary least squares often doesn't extend to survival modeling.

The baseline hazard over time is necessarily 0 except at observed event times, as the Cox-model hazard is exactly 0 between event times. The lines connecting the points in your first plot are thus, at best, misleading. If you want to show the baseline hazard then show the cumulative hazard, a step function. Remember, however, that the baseline hazard isn't even a part of the Cox model itself. It's something you can extract from the data after you fit the model. It doesn't evaluate the fit of the model.

The discussion of goodness-of-fit measures for survival models in lifelines doesn't include deviance residuals, for good reason. Those coming new to survival modeling often hope to use deviance residuals (or the martingale residuals to which they are related) in the ways that one uses residuals in least-squares modeling. For many purposes you can't, as Therneau and Grambsch explain in Chapter 4. For example, in Section 4.2.3, they note that "martingale residuals and the fitted values are negatively correlated, frequently strongly so," which is related to the plot you show.

With time-fixed covariates, the deviance residuals $d_i$ are close to a normalized difference between the observed $(N_i)$ and predicted $(\hat E_i)$ numbers of events for an individual at the event or censoring time $t_i$ (Therneau and Grambsch, Section 4.3):

$$d_i \approx \frac{N_i(t_i)-\hat E_i(t_i)}{\sqrt {\hat{E_i}(t_i)}}. $$

The numerator is the martingale residual process evaluated at the event or censoring time. The estimated expected number of events for individual $i$ with time-fixed covariate values $X_i$ and corresponding coefficient estimates $\hat \beta$ is a non-decreasing function of time (see Chapter 4 of Therneau and Grambsch):

$$\hat{E_i}(t) = \int_0^t Y_i(s) e^{X_i \hat\beta} d \hat{\Lambda}_0(s),$$

where $Y_i(s)$ is 1 while the individual is at risk (0 otherwise) and $\hat{\Lambda}_0(s)$ is the estimated baseline cumulative hazard over time.

$\hat{E}_i(t)$ is thus related to the estimated baseline cumulative hazard, which has a step increase at each event time over the entire data sample. That has two implications.

First, $\hat{E}_i(t)$ increases by $e^{X_i \hat\beta} d \hat{\Lambda}_0(s)$ at each event time $s$ while at risk. For any set of individuals with the same set of covariate values, the deviance residuals thus necessarily decrease with increasing times to events (or to right censoring).

Second, as the baseline cumulative hazard estimate extends out until the last event in the entire data set, it's quite possible for $\hat{E}_i(t)$ to exceed 1 at a late enough time, if $e^{X_i \hat\beta}$ is large enough.

Put another way, a survival model is of a distribution of event times as a function of covariate values, typically a pretty wide distribution. Some individuals with the same covariate values are going to have event times earlier than "expected" and others will have later event times. That (plus the fact that the deviance residual can't exceed 0 if the event time is censored) is all that your plot of deviance residuals against observation time demonstrates.

Plots of martingale residuals against values of a covariate can be useful to evaluate a model's functional form for a continuous covariate; see this page. What can work more directly than starting from martingale residuals is to use a flexible regression spline to fit a continuous predictor like year_of_manufacture and let the data tell you an approximate functional form.

Finally, the concordance doesn't evaluate "average values" in the way that you seem to think. It's the fraction of pairs of cases for which the observed and model-predicted event orders agree (among case pairs that can be evaluated). It's a measure of discrimination among cases, not calibration. In your case, only 0.56 of pairs had agreement between observed and predicted event orders, close to the 0.5 you would expect by chance.

  • $\begingroup$ Thank you so much for your very helpful explanations! I understand that the "lines" I see in my deviance residuals plot come from different parameter combinations. 1. Can you explain once more why the residuals are necessarily higher for smaller event times? 2. The plot deviance residuals vs event time is then only useful to detect outliers? I understood that it might also make sense to plot deviance residuals against the predictors to see how those predictors might be related to the target variable. I'm going to try this now. $\endgroup$
    – Requin
    Commented Mar 20, 2023 at 9:31
  • 1
    $\begingroup$ @Requin Therneau and Grambsch say in Section 4.3 about the deviance residual: "In practice it has not been as useful as anticipated" even in plots. For plots against predictor values, martingale residuals have a stronger theoretical basis. In the approximation shown, $\hat E_i$ is the expected number of events for individual $i$ at the individual's event/censoring time. For a given set of covariate values, an early event ($N_i=1$) might have $\hat E_i=0.1$ and a late event time have $\hat E_i=0.9$. Corresponding approximate $d_i$ are 2.8 and 0.1, respectively. $\endgroup$
    – EdM
    Commented Mar 20, 2023 at 13:34
  • $\begingroup$ Thanks! I guess I didn't understand the notion "expected number of events" for an individual. Intuitively, I'd say this is something like "1 minus the survival curve of an individual", but with this, E^i>0 is not possible. So, where do I get the "expected number of events" for an individual from? $\endgroup$
    – Requin
    Commented Mar 20, 2023 at 19:15
  • 1
    $\begingroup$ @Requin I've added the formal definition of the expected number of events to the answer. It increases at each event time in the data set while an individual is at risk, and can exceed 1 if there's a high estimated hazard ratio and the data set has enough events. $\endgroup$
    – EdM
    Commented Mar 20, 2023 at 20:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.