The proportional hazards model proposed by Cox (1972) has been widely used in modeling censored survival data. It is both easy to implement and easy to interpret, usually making it the biostatistician’s first model to attempt when faced with time-to-event or survival data. In addition, its usage has extended to fields beyond biostatistics, such as predicting bank failures in finance (Lane et al., 1986), estimating customer attrition probability in insurance (Kang and Han, 2004), identifying determinants for duration of unemployment in labor market research (Kupets, 2006), and modeling time until a policy is adopted in political science (Jones and Branton, 2005). It has been deemed one of the “breakthroughs in statistics” (Kotz and Johnson, 1992), and has been cited over 46,256 times.

Due to its pervasive applicability, before taking the results from a fitted Cox model as valid, one should address a few questions: is the proportional hazards assumption satisfied? Are the functional forms of the variables appropriate? Are there any outliers or influential observations? To answer these questions, multiple methods have been proposed, many of which rely on different types of residuals of the model.

In this article, we review methods for assessing goodness-of-fit of the Cox proportional hazards model. The rest of this article is organized as follows. In Section 2, we present the notation and important quantities in estimating a Cox model. In Section 3, we discuss the potential problems, and review the literature on tests and graphical methods proposed to identify them. The diagnostic plots and tests for each problem are illustrated with an application regarding dental clinic visits using existing R software in Section 4. We conclude with a brief discussion and suggestions for possible remedies when nonproportionality is identified.

Let Ti* be the true time to event for individual i, i = 1 . . . , n. Define T_i to be the minimum of Ti* and censoring time C_i, which is a time beyond which individual i cannot be observed, i.e., Ti=min(Ti*,Ci). We assume C_i is independent of Ti* for each i = 1, . . . , n. Define δ_i as an indicator equal to 1 if Ti* is observed for individual i and 0 if the observation is censored, i.e., δi=I(Ti*<Ci). Additionally, let X_i be the p-dimensional vector of covariates for the i^th individual. Assume we observe (T_i, δ_i, X_i) independently for i = 1, . . . , n.

The Cox model as given in the 1972 manuscript specifies the hazard for individual i as

where λ₀ is an unspecified, nonnegative function called the baseline hazard, and β is a p-dimensional vector of coefficients. The model is also known as the proportional hazards model, since the hazard ratio for two subjects with covariate vectors X_i and X_j is

which is independent of time. In the case of a single binary predictor, β is the logarithm of the hazard ratio between the corresponding two subgroups of the data.

The Cox model has subsequently been extended to incorporate time-dependent covariates. We will use the notation X_i(t) to allow for the possibility of time-dependent covariates. In this case, the hazard ratio for two subjects with covariate vectors X_i(t) and X_j(t) becomes

such that the model assumes the covariates have proportionate effects on the hazard function over time (see, e.g., Fisher and Lin, 1999).

Estimation of β is based on the partial likelihood introduced by Cox (1972) and later formulated by Cox (1975):

where δ_i is the censoring indicator, and Y_j(t) = I(T _j > t) indicates whether the j^th subject is still at risk (alive) at time t.

Assume there are no ties in the failure time data. Let N_i(t) be the number of events for subject i at time t, which is either 0 or 1, and define

where Δt is chosen to be sufficiently small so that ∑i=1ndNi(t) is either 0 or 1 at any time t. Using the counting process formulation of Fleming and Harrington (1991), the partial likelihood (2.2) can be alternatively written as

where γ_i(β, t) is the risk score for subject i, γ_i(β, t) = exp{β′X_i(t)} ≡ γ_i(t). We will use this form for further derivations. Taking the logarithm of (2.3), we have the log partial likelihood:

Differentiating (2.4) with respect to β yields the p × 1 score vector U(β):

where χ̄(β, s) is a weighted average of X over those observations still at risk at time s with weights Y_i(s)γ_i(s):

The p × p information matrix is the negative second derivative of (2.4), given by

where V(β, s) is the weighted variance of X at time s:

We obtain the maximum partial likelihood estimator by solving the partial likelihood equation

The solution β̂ is consistent and asymptotically normal with mean β, the true parameter vector, and variance {ℰℐ(β)}⁻¹, the inverse of the expected information matrix. While the expected information matrix is unknown, the observed information matrix ℐ(β̂) is available as

We may use its inverse, ℐ⁻¹(β̂), as the estimated variance of β̂.

Most diagnostic procedures for the Cox model are based on its residuals. There are multiple types of residuals defined for the Cox model, and they can often serve different purposes in model diagnostics. In this section, we review methods to identify (1) violations of the proportional hazards assumption, (2) appropriate functional forms of covariates, (3) outlying observations, and (4) influential observations. We demonstrate the usage of many of these methods in Section 4, and provide implementation details in the Supplemental Materials.

3.1. Proportional hazards assumption

3.2. Functional forms

Martingale residuals, defined in Equation (3.6), play an important role in functional form diagnostics. Barlow and Prentice (1988) provided more detailed discussion and illustrated that plots of such residuals may provide insight to the choice of model form. Therneau et al. (1990) discussed the usage of martingale residuals in investigating the functional form of covariates. To examine a particular covariate, they suggest fitting a proportional hazards model omitting that covariate and computing the martingale residuals M̂_i as given in Equation (3.6). Then a smoothed plot of M̂_i versus the omitted covariate often gives approximately the correct functional form of the covariate (e.g., linear, quadratic) to place in the exponent of a Cox model. They also pointed out, however, that this plot did not work well when dealing with large covariate effects, and that it requires the covariate of interest to be uncorrelated with other covariates in the model.

Henderson and Milner (1991) noticed that plots of the martingale residuals against time, although useful, can exhibit systematic patterns which are not a priori predictable even when the model fails. They suggested two amendment approaches and gave an example for illustration. One approach was to superimpose the estimated mean when plotting residuals, which enables comparison between the observed patterns and the expected patterns. The other approach was to subtract the conditional expected value from each observed residual and scale it using its standard deviation, which could be consistently estimated from the data according to Barlow and Prentice (1988). Then the standardized residuals, when plotted, should be randomly scattered if the model is appropriate.

Grambsch (1995) proposed two aspects from which the martingale residual plot in Therneau et al. (1990) can be improved. One aspect is to modify the martingale plot for counting process data because of the close relationship between counting process models and Poisson regression. Suppose Z is the variable of interest. If a monotonic relationship between Z and the hazard λ_i(t) is expected, a log-linear form is often adequate. The model

is fitted, and the expected count for the i^th individual is

where λ̂₀(t) is the estimated baseline hazard and α̂ is the estimated parameter for Z. The martingale residual in this case would be M̂_i = δ_i − Ê_i, and the generalized linear model (GLM) partial residual is given by

McCullagh and Nelder (1983) recommended plotting the partial residual against Z as an informal check for the correctness of the guess for functional form.

The other aspect mentioned by Grambsch (1995) comes from the penalized likelihood approach of Hastie and Tibshirani (1993). They assumed that the functional form of covariate X_j is an unknown, smooth function f _j, and proposed the alternative formulation

which enables estimation of all functional forms at the same time. To avoid overfitting, they maximized the penalized partial likelihood with penalty ∑j=1pνj∫fj″(s)2ds, where ν_j ≥ 0, j = 1, . . . , p, are smoothing parameters that can be tuned. Both approaches lead to approximately the same solution, but the latter is computationally more complex since the optimization is done within the kernel of the partial likelihood.

3.3. Outlying observations

A plot of martingale residuals against the linear prediction β̂′X(t) or the risk score γ_i(t) defined in Section 2 often helps to identify the observations who have died too soon or lived too long, based on the assumed model. Nevertheless, having a range of (−∞, 1], the martingale residual is often heavily skewed, and may be misleading. Therneau et al. (1990) used a liver disease data set to demonstrate these scenarios, where the martingale residual plot indicated that some observations died too soon while in actuality they were not outliers at all. They pointed out that it is a favorable practice to transform the residuals to a more normal shaped distribution to help assess the prediction accuracies for individual subjects.

Inspired by the deviance residuals for GLM in McCullagh and Nelder (1983), Therneau et al. (1990) introduced the deviance residual for a Cox model:

where δ_i is again the censoring indicator for subject i. From the functional form it is apparent that the deviance residual is essentially a transformation of the martingale residual. Therneau et al. (1990) concluded that with less than 25% of censoring, the deviance residual is approximately normally distributed. With censoring rates greater than 40%, too many points will lie near 0 and make the distribution not normal, but the set of residuals is still symmetrized. Plotting d_i against β̂′X_i or exp(β̂′X_i) will help identify potential outliers which have deviance residuals with too large absolute values.

Noticing that deviance residuals do not have a reference distribution and the normal approximation can sometimes be unsatisfactory (Fleming and Harrington, 1991), Nardi and Schemper (1999) proposed two new types of residuals: (i) the log-odds residual L_i = log [S _i(T_i)/{1 − S _i(T_i)}] and (ii) normal deviate residual η_i = Φ⁻¹{S _i(T_i)}, i = 1, . . . , n, where Φ⁻¹ is the inverse normal cumulative distribution function. Assuming S _i(·) is known, the sampling distribution for L_i is logistic with E(L_i) = 0 and var(L_i) = π²/3, and standard normal for η_i. In practice, we use the predicted survival for observation i, Ŝ _i(T_i), to calculate L̂_i and η̂_i, which converge in probability to L_i and η_i, respectively. Based on simulations, they concluded the performances of these two residuals when identifying outliers are better than that of the deviance residual since they are both unimodal, and the empirical distribution of deviance residual often becomes bimodal because of censoring. They suggested that one can use the quantiles of the normal distribution, ±1.64 and ±1.96, and of the logistic distribution ±2.94 and ±3.66, to help identify potential outliers.

3.4. Influential observations

The score vector U defined in Equation (2.5) is of great importance in influential diagnostics. Again, using the counting process formulation, the score residual for the i^th individual is defined to be

where χ̄(β̂, s) is the χ̄(β, s) defined in Equation (2.6) evaluated at β = β̂, and M̂_i(s) is defined in Equation (3.5).

In studying the influence of one observation, a general practice is to delete that observation, fit the model again, and compare the parameter estimates with those of the model fit on the complete data. Nevertheless, the Cox model is conceptually different from linear or generalized linear models in that it involves both parametric and nonparametric estimation. Therefore, an observation could be influential in terms of more than just regression coefficients. We review measures of both in this section.

The dental restoration longevity data, provided by the University of Iowa College of Dentistry’s Geriatric and Special Needs (SPEC) Clinic (see Caplan et al., 2017), is used as a case study to demonstrate the diagnostic methods of the Cox model. For this analysis, electronic data was obtained during the 5-year period from 1/1/1995–12/31/1999. The health record numbers were scrambled by IT personnel to ensure that no Personal Health Information was included. Subsequently, the Institutional Review Board at the University of Iowa declared that this project is exempt from Human Subjects Review, due to the anonymous nature of the data.

We identified 697 unique patients who went to the SPEC Clinic to treat their molars upon their first visit and received restoration in amalgam, composite, or glass ionomer. The follow-up of their visits began on the date of restoration. Any restoration that was replaced with another intracoronal or extracoronal restoration, accessed for endodontic therapy, or extracted was deemed to have undergone an event. If the restoration results in an event, the event date would become the end of follow-up. Restorations that did not incur an event are considered censored up to the date of the patient’s second-to-last visit to any College of Dentistry’s clinic. Among the 697 patients, 228 experienced an event during the follow-up, giving a censoring rate of 67.3%.

We considered the following covariates: Gender, Age at receiving restoration (centered and scaled), Occupation (Faculty, Non-faculty) and Size (Small, Medium, Large). Analysis was performed using the survival package in R, and figures were produced using the survminer (Kassambara and Kosinski, 2017), ggplot2 (Wickham, 2009) and ggfortify (Tang et al., 2016) packages. R code for analysis is available in our supplementary material.

4.1. Functional form

As suggested in Section 3.2, we should determine the appropriate form of covariates to include in the model before testing for proportionality. Age is the only continuous covariate whose form needs to be assessed. We use the methods of Therneau et al. (1990): fit a model excluding Age and obtain its martingale residuals. The martingale residuals are plotted against Age in Figure 1. We also superimpose the loess pointwise confidence band. The curvy behavior of the loess fit indicates that we should consider higher orders of Age.

4.2. Proportional hazards

As suggested in Section 4.1, we consider including the square of Age (Age²) in the model. We fit two models: one with only linear Age effects (Model 1) and another model with linear and quadratic Age effects (Model 2) to assess the improvement to the model when correcting the functional form. To assess the proportional hazards assumption, we used the cox.zph function from the survival package to obtain both the individual χ12 statistics for each covariate and the global χp2 statistic for each model. The test results are summarized in Table 2. While both models passed the global test, Age in Model 1 did not pass the individual test at the 0.05 level. In Model 2, however, both Age and Age² pass the individual test of proportionality at the 0.05 level.

The parameter estimates of Model 2 are summarized in Table 3. Restorations for males tend to fail later than for females, while restorations for older patients tend to fail sooner. Compared to large restorations, medium and small restorations are less likely to fail.

As mentioned in Section 3.1.1, when the proportional hazards assumption holds, the Schoenfeld residuals will be close to zero. Therefore a plot of the Schoenfeld residuals against survival times would be informative (Schoenfeld, 1982). Figure 2 shows the plots for Model 2. For all six covariates, the smoothed pointwise confidence bands are all around 0, which again confirms that there is no obvious evidence against the proportional hazards assumption.

For the three categorical covariates (Gender, Occupation, and Size), we also utilize the graphical methods in Section 3.1.4 to check the proportional hazards assumption. For each covariate, we plot the estimated survival curves (Ŝ (t)), the cumulative hazards (− log Ŝ (t)) and the log-log transformed survival (log(− log Ŝ (t))) against survival times in Figure 3. The three plots for Gender indicate that the hazards of the two gender strata are proportional, but the lack of large discrepancy indicates that this proportional effect is not significant. Similarly, the ignorable discrepancy between the two occupation strata tells the same story. The three plots for restoration size strata, however, are more informative, in that although the proportionality effect is small between SizeLarge and SizeMedium, it is highly significant between SizeLarge and SizeSmall.

As mentioned in Section 3.1, the martingale residual can be used to graphically assess the proportional hazards assumption as well. We plot the cumulative sum of martingale residuals ordered by Age in Figure 4. The curve fluctuates around zero as expected.

4.3. Outlying observations

As suggested in Section 3.3, both the martingale residual and the deviance residual are useful for identifying outlying observations, but the deviance residual is less skewed and therefore more useful. We plot both residuals against the linear predictions, β̂′X, in Figure 5. In Figure 5(a), the martingale residuals do not vary much against the linear predictions, and fail to identify any outlying observations. Using ±1.96 as thresholds, the deviance residuals plotted in Figure 5(b) identify 34 potential outliers. Upon further investigation, these subjects turned out to be much younger than other subjects (46.6 vs. 55.1) but their restorations failed very soon. Due to the high censoring rate, however, the normal-approximation-based thresholds may not be appropriate.

We also use the log-odds residual and the normal deviate residual discussed in Section 3.3 to look for potential outliers. Both the log-odds residual in Figure 6(a) and the normal deviate residual in Figure 6(b) identify the same set of 67 potential outliers, which is bigger than the set of outliers identified by the deviance residual. This set, however, still consists of younger individuals (51.7 vs. 55.1) whose restorations failed very soon.

4.4. Influential observations

We use the methods in Section 3.4 to perform influential diagnostics. We first look at influence of observations on parameter estimates and plot the dfbetas residuals in Figure 7. As illustrated, no observation caused any parameter change of more than 15% of that parameter’s standard error. Considering that there are 697 observations, we can conclude there are no significantly influential observations.

We also present the likelihood displacement approach in Figure 8. The absence of particularly large likelihood displacements further confirms our conclusion from Figure 7.

With such wide usage across a variety of disciplines, the importance of Cox regression for modeling time-to-event data cannot be overstated. As a consequence, one must consider the validity of the results from such an analysis before making any conclusions. This article presents a review of diagnostic methods for the Cox proportional hazards model, including methods for identifying violations of the proportional hazards assumption, appropriate functional forms of continuous covariates, outlying observations, and influential observations. Using a non-linear functional form of covariates can often improve model fit, while detected outlying or influential observations should be investigated further before taking any action.

Violations of the proportional hazards assumption can be addressed in several ways, the most common include the use of time-varying coefficients and stratified models. Time-varying coefficients, which are left out in this article, make a major source of nonproportionality and are often the alternatives to test against when assessing the proportional hazards assumption. More flexible models that incorporate this effect have been studied by Murphy and Sen (1991), Hastie and Tibshirani (1993), Verweij and van Houwelingen (1995), Sargent (1997), Marzec and Marzec (1997b), Cai and Sun (2003), Tian et al. (2005), Fan et al. (2006), and more recently by Chen et al. (2012). In practice, the graphical tools in the survival package enable us to check if there are any time-varying coefficients (Therneau et al., 2017). Another popular approach for addressing non-proportionality is using a stratified Cox model. In this case, it is assumed that individuals in different strata have different baseline hazard functions, but all other predictor variables satisfy the proportional hazards assumption within each stratum. Related chapters can be found in Therneau and Grambsch (2000), Kalbfleisch and Prentice (2002), Lawless (2003), and Collett (2015).

The Cox model has also been extended to the analysis of interval-censored survival data. Such models have been studied by Finkelstein (1986), Farrington (2000), Goggins and Finkelstein (2000) and recently Heller (2010). In particular, Farrington (2000) proposed the counterparts to the Cox-Snell, martingale, deviance, and Schoenfeld residuals and illustrated their usage in model diagnostics under the interval-censored framework.

The authors are very grateful to Professor Jun Yan for his helpful discussions and comments on an earlier version of this manuscript.