3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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

enter image description here

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.

$\endgroup$
  • $\begingroup$ I guess fitting 'weibull' distribution increases the variance close to one ( 0.8302648) holding the mean unity. $\endgroup$ – KDG Oct 8 '18 at 8:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy