R/AFTresiduals.R

In sa4khan/AFTjmr: Recurrent event analysis and joint modeling using accelerated failure time (AFT) models

#################################################
#' Cox-Snell residuals
#' @description Cox-Snell residuals of the \code{survreg.aft} fit.
#' @keywords Cox-Snell residual
#' @param fit \code{survreg.aft} fit.
#' @param plot returns residual plot if \code{TRUE} (default is \code{FALSE}).
#' @param xlim x limits (x1, x2) of the plot (optional).
#' @param ylim y limits of the plot (optional).
#' @param main main title of the plot (optional).
#' @details Plot of the Cox-Snell (hazard-based) residuals should be roughly a straight line 
#'     with unit slope when the model is adequate. 
#'     See Section 6.2 of Lawless and Section 3.7.3 of Cook and Lawless for details.
#' @return \code{residuals} Cox-Snell residuals.
#' @return \code{plot} residual plot if \code{plot=TRUE}.
#' @references Cook RJ and Lawless J, The statistical analysis of recurrent events, Springer, 2007.
#' @references Lawless J, Statistical models and methods for lifetime data, Wiley, 2003.
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{survreg.aft}}, \code{\link{LR.test}}.
#' @examples 
#' # Example 1: Recurrent Event Analysis
#'   library(frailtypack)
#'   data(readmission)
#' # Recoding sex and chemo
#'   sex<-factor(2-as.numeric(readmission$sex),
#'        labels=c("Female","Male"))
#'   chemo<-factor(as.numeric(readmission$chemo)-1,
#'          labels=c("NonTreated","Treated"))
#'   Data<-readmission
#'   Data$sex<-sex
#'   Data$chemo<-chemo
#'   fit.gg<-survreg.aft(c(t.start,t.stop,event)~sex+chemo+charlson,
#'       Data=Data,model="ggamma")
#'   fit.w<-survreg.aft(c(t.start,t.stop,event)~sex+chemo+charlson,
#'       Data=Data,model="weibull")
#'   surv.resid(fit.w)
#'   par(mfrow=c(1,2))
#'   surv.resid(fit.w,plot=TRUE,xlim=c(0,4.75),ylim=c(0,4.75),main="Weibull")
#'   surv.resid(fit.gg,plot=TRUE,xlim=c(0,4.75),ylim=c(0,4.75),main="Generalized Gamma")
#'
#' # Example 2: Non-recurrent Data (Right Censored)
#'   library(survival)
#'   fit.ll<-survreg.aft(c(time,status)~karno+diagtime+age+prior+celltype+trt,
#'       Data=veteran,model="llogistic")
#'   surv.resid(fit.ll,plot=TRUE,main="Log-logistic")
#' @export

surv.resid<-function(fit,plot=FALSE,xlim=NULL,ylim=NULL,main=NULL){
if(fit$surv.model=="ewaft") {resid<-ewaft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
if(fit$surv.model=="Rewaft") {resid<-Rewaft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
if(fit$surv.model=="waft") {resid<-waft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
if(fit$surv.model=="Rwaft") {resid<-Rwaft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
if(fit$surv.model=="ggaft") {resid<-ggaft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
if(fit$surv.model=="Rggaft") {resid<-Rggaft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
if(fit$surv.model=="llaft") {resid<-llaft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
if(fit$surv.model=="Rllaft") {resid<-Rllaft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
if(fit$surv.model=="lnaft") {resid<-lnaft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
if(fit$surv.model=="Rlnaft") {resid<-Rlnaft.resid(fit,plot=plot,xlim=xlim,ylim=ylim,main=main)}
return(resid)
}