# Cox Snell residuals in R

I'm trying to fit parametric survival models for a data that I have, and I don't know how to get the Cox Snell residuals in R.

An example of dataset with Exponential model fitted

install.packages("survival")
library(survival)
data(lung)
fit<-survreg(Surv(time,status)~age+as.factor(sex)+ph.karno,na.action = "na.omit",
dist = "exponential",data=lung)


I want to use the Cox Snell residuals to check if the model distribution is good for the data. I will fit a Exponential, Log-Normal, Weibull and Log-logistic model and do these residuals for all them.

So I want to find the Cox Snell residuals and check if these residuals follow a exponential distribution with parameter 1.

A model fits the data well if the Cox-Snell residuals follow an exponential distribution of parameter 1; the Komologorov-Smirnov Goodness of Fit Test (KS-test) is used to assess whether this is the case.

I will do the Komolgorov test later, but I want to check the residuals first.

Based on the book's formulas, I wrote this code :

#rC : Cox-Snell residuals
#rM : Martingale residuals
#rD : Deviance residuals

rC<-exp(((fit$$y[,1])-log(predict(fit,lung,na.action = "na.omit")))/fit$$scale)
rM<-fit$$y[,2]-rC rD<-sign(rM)*sqrt(-2*(rM+fit$$y[,2]*log(rC)))  # -residuals(fit,type='deviance')

mean(rC)
var(rC)

qqplot((qexp(ppoints(length(rC)))),(rC));qqline(rC, distribution=qexp,col="red", lty=2)


(The exponential qqplot reference is here: https://stackoverflow.com/a/37031433/10042541)

And this is the result :

> mean(rC)
[1] 0.722467
> var(rC)
[1] 0.3411307


This graph seems not bad but the mean and the variance isn't close to 1 and deletion of the 6th row which looks like outlier lower the variance. I guess the residuals are more spread out than the exponential distribution with parameter 1 and some of them are smaller than they supposed to be.

According to the reference book, we need to modify rC for censored observations because they will be too small.

rC<-rC+(1-fit\$y[,2])*1  # or *log (2) instead of *1


This makes the mean to be unity though it makes the qqplot seem worse.

I haven't studied the survival analysis thoroughly so I am not sure how to improve this model. Any improvement on this answer should be welcomed. Thank you very much.

• I guess fitting 'weibull' distribution increases the variance close to one ( 0.8302648) holding the mean unity.
– KDG
Commented Oct 8, 2018 at 8:40

The Cox-Snell residual for case $$j$$ in a survival model is $$r_j=\hat H(T_j|X_j)$$, where $$\hat H(T|X)$$ is the estimated cumulative hazard function, $$T_j$$ is the event/censoring time for the case and $$X_j$$ its vector of covariate values. That's a convenient form for proportional hazard models, for which $$\hat H(T_j|X_j)= \hat H_0(T) \exp(X_j' \hat\beta)$$, the product of an estimated baseline hazard $$H_0(T)$$ and the hazard ratio for the case calculated from the vector $$\hat \beta$$ of regression coefficients. That simple multiplicative form for $$\hat H$$, however, doesn't hold without proportional hazards.

As $$S(T|X)=\exp(-H(T|X))$$ for continuous-time survival models, once you have a parametric equation for $$\hat S(T|X)$$ you can calculate a Cox-Snell residual directly from that equation: $$r_j =\hat H(T_j|X_j) = -\ln(\hat S(T_j|X_j))$$. The Survival() function in the R rms package, for example, provides an easy way to get survival probability estimates from a psm parametric model fit for specified times and linear-predictor values, which then can be transformed into Cox-Snell residuals.

There is a simple relationship between Cox-Snell residuals and standardized residuals for models that can be expressed in a location-scale form:

$$f(T) \sim X' \beta + \sigma W,$$

where $$\sigma$$ is a scale factor, $$W$$ is the standard probability distribution corresponding to the parametric survival model, and $$f(T)$$ is a link function, typically $$\ln(T)$$. With a log link, $$W$$ is standard minimum extreme value for Weibull and exponential models ($$\sigma = 1$$ for exponential), standard logistic for log-logistic, and standard normal for lognormal. These notes provide a concise introduction to such parametric modeling.

The survival function $$S(T)$$ is the complement of the cumulative distribution function (CDF) of the survival times, so the Cox-Snell residual can be written $$r_j = -\ln(1-\widehat {\text{CDF}}(T_j|X_j))$$. For a location-scale model with distribution $$W$$, $$\widehat {\text{CDF}}(T_j|X_j)$$ can be calculated from the standardized residuals

$$s_j=\frac {f(T_j)-X_j' \hat \beta}{\hat \sigma},$$

as those $$s_j$$ should be distributed according to $$W$$. For example, if you want a Cox-Snell residual for a lognormal model, $$r_j = -\ln(1-\Phi(s_j))$$ where $$\Phi$$ is the standard normal CDF. This 1:1 relationship between the residual types means that evaluation of a model via Cox-Snell residuals $$r_j$$ is equivalent to evaluation via standardized residuals $$s_j$$, as Klein and Moeschberger note on page 415.

The answer from @KDG shows problems with using Cox-Snell residuals to evaluate a model. First is how to deal with censoring. Second is the nature of the plots typically used to evaluate the agreement with the theoretical exponential form, which tend to overemphasize the larger values at the upper right while squishing together most of the cases.

Both of these problems are solved by working directly with the (potentially censored) standardized residuals $$s_j$$ instead. You examine the survival function of the $$s_j$$ (the complement of the CDF) with a Kaplan-Meier plot that incorporates the censoring and compare that against the survival function of the standard distribution $$W$$. That both takes care of the censoring problem and spaces cases more evenly.

Harrell's rms package provides a simple implementation of standardized residual analysis for standard parametric families. Starting with a psm object from an rms parametric survival model fit, a residuals() function provides standardized residuals along with censoring indicators in a Surv object. The survplot() function applied to that Surv object of (censored) standardized residuals displays the Kaplan-Meier (KM) curve of the $$s_j$$ along with the theoretical CDF based on the assumed distribution $$W$$. For the example in the OP:

where the thin lines are the KM plots and the thick lines are the theoretical forms for the indicated parametric families. The inadequacy of the exponential model and the quality of the Weibull model for these data are apparent. I don't see any reason to use Cox-Snell residuals for such parametric models.

Code for plots:

library(survival)
library(rms)
lungCC <- lung[complete.cases(lung[,c("age","sex","ph.karno")]),] ## instead of using na.action
psmE <- psm(Surv(time,status)~age+sex+ph.karno,dist="exponential",data=lungCC)
residE <- residuals(psmE)
psmW <- psm(Surv(time,status)~age+sex+ph.karno,dist="weibull",data=lungCC)
residW <- residuals(psmW)
psmLN <- psm(Surv(time,status)~age+sex+ph.karno,dist="lognormal",data=lungCC)
residLN <- residuals(psmLN)
psmLL <- psm(Surv(time,status)~age+sex+ph.karno,dist="loglogistic",data=lungCC)
residLL <- residuals(psmLL)
par(mfrow=c(2,2))
survplot(residE,main="Exponential",ylab="Complement of residual CDF")
survplot(residW,main="Weibull",ylab="Complement of residual CDF")
survplot(residLN,main="Lognormal",ylab="Complement of residual CDF")
survplot(residLL,main="Log Logistic",ylab="Complement of residual CDF")