R/JMresiduals.R

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

###########################################################
#' Cox-Snell residual plot for the event process of the joint model
#' @description Returns Cox-Snell residual plot.  
#' @keywords Joint model, Cox-Snell residual
#' @param bayes.fit a \code{list} of \code{jmreg.aft} fits.
#' @param summary \code{"posterior.mean"} or \code{"posterior.median"} of the Cox-Snell residuals are used 
#'        (default is \code{"posterior.mean"}).
#' @param col (optional) line colors for multiple fits (see "Color Specification" in \code{par})
#' @details In MCMC simulations, log(survival probabilities) are monitored to obtain posterior means/medians of the Cox-Snell residuals. 
#' @return Returns Cox-Snell residual plot. The plot should be roughly a straight line 
#'     with unit slope when the survival model is adequate. For a list of multiple fits, the plots are superimposed
#'     on the same graph to facilitate comparison.
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{jmreg.aft}}, \code{\link{jm.resid}}, \code{\link{jm.summary}}
#' @examples 
#'   library(JM)
#'   data(aids.id)
#'   data(aids)
#'   surv.fit<-coxph(Surv(Time,death)~drug+gender+prevOI+AZT,
#'           data=aids.id,x=TRUE)
#'   lme.fit<-lme(CD4~obstime+obstime:drug+gender+prevOI+AZT,
#'          random=~obstime|patient,data=aids)
#'   jmfit.ew<-jmreg.aft(surv.fit=surv.fit,lme.fit=lme.fit,
#'          surv.model="eweibull",rand.model="simple",timevar="obstime",
#'          n.chain=2,n.adapt = 5000, n.burn = 15000, 
#'          n.iter = 150000, n.thin = 5)
#'   jmfit.ln<-jmreg.aft(surv.fit=surv.fit,lme.fit=lme.fit,
#'          surv.model="lnormal",rand.model="simple",timevar="obstime",
#'          n.chain=2,n.adapt = 5000, n.burn = 15000, 
#'          n.iter = 150000, n.thin = 5)
#'   par(mfrow=c(1,2))
#'   jm.resid.plot(list(jmfit.ew))
#'   jm.resid.plot(list(jmfit.ln))
#'   jm.resid.plot(list(jmfit.ew,jmfit.ln))
#' @export


jm.resid.plot<-function(bayes.fit,summary="posterior.mean",col=NULL){
if(length(bayes.fit)==1){
resid<-bayes.fit[[1]][["MCMC output"]][["cox.snell.residuals"]][,summary]
delta<-bayes.fit[[1]][["MCMC output"]][["surv.data"]][,"status"]
km.Rhat<-survfit(Surv(resid,delta)~1,type="kaplan-meier")
H<-cbind(km.Rhat$time,-log(km.Rhat$surv))
H[which(!is.finite(H))]<-NA
H<-na.omit(H)
xlim<-ylim<-c(min(c(H)),max(c(H)))
plot(H[,1],H[,2],type="p",pch=20,xlab="Residual",
  ylab="Estimated Cumulative Hazard",
  xlim=xlim,ylim=ylim,main=bayes.fit[["MCMC output"]]$surv.model.full,
  cex.lab=1.25,cex.main=1.5)
abline(a=0,b=1)
}
if(length(bayes.fit)>1){
H1<-NULL
H2<-NULL
for(i in 1:length(bayes.fit)){
resid<-bayes.fit[[i]][["MCMC output"]][["cox.snell.residuals"]][,summary]
ordered<-order(resid)
resid<-resid[ordered]
delta<-bayes.fit[[i]][["MCMC output"]][["surv.data"]][,"status"][ordered]
km.Rhat<-survfit(Surv(resid,delta)~1,type="kaplan-meier")
H1<-cbind(H1,km.Rhat$time)
H2<-cbind(H2,-log(km.Rhat$surv))
}
H<-cbind(H1,H2)
H[which(!is.finite(H))]<-NA
H<-na.omit(H)
xlim<-ylim<-c(min(c(H)),max(c(H)))
HH1<-H[,1:length(bayes.fit)]
HH2<-H[,(length(bayes.fit)+1):ncol(H)]
plot(HH1[,1],HH2[,1],type="l",xlab="Residual",ylab="Estimated Cumulative Hazard",
  lty=1,xlim=xlim,ylim=ylim,main="",
  cex.lab=1.25,cex.main=1.5,col=ifelse(is.null(col),1,col[1]))
abline(a=0,b=1)
s.model<-bayes.fit[[1]][["MCMC output"]]$surv.model.full
for(i in 2:length(bayes.fit)){
lines(HH1[,i],HH2[,i],type="l",lty=i,col=ifelse(is.null(col),i,col[i]))
s.model<-c(s.model,bayes.fit[[i]][["MCMC output"]]$surv.model.full)
}
ccol<-if(is.null(col)) (1:length(bayes.fit)) else col
legend("topleft",s.model,lty=1:length(bayes.fit),col=ccol)
}
}
###################################################################
# Cox-Snell Residuals
############################################################
#' Cox-Snell residuals for the event process of the joint model
#' @description Returns Cox-Snell residuals.  
#' @keywords Joint model, Cox-Snell residual
#' @param bayes.fit a fit of the joint model using \code{jmreg.aft}.
#' @param posterior.mean returns posterior means if TRUE, and posterior medians if FALSE.
#' @details In MCMC simulations, log(survival probabilities) are monitored to obtain posterior means/medians of the Cox-Snell residuals. 
#' @return Returns posterior means/medians of the Cox-Snell residuals.
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{jmreg.aft}}, \code{\link{jm.resid.plot}}, \code{\link{jm.summary}},
#' @examples
#'   library(JM)
#'   data(aids.id)
#'   data(aids)
#'   surv.fit<-coxph(Surv(Time,death)~drug+gender+prevOI+AZT,
#'           data=aids.id,x=TRUE)
#'   lme.fit<-lme(CD4~obstime+obstime:drug+gender+prevOI+AZT,
#'          random=~obstime|patient,data=aids)
#' # exponentiated weibull time-to-event process
#'   jmfit.ew<-jmreg.aft(surv.fit=surv.fit,lme.fit=lme.fit,
#'          surv.model="eweibull",rand.model="simple",timevar="obstime",
#'          n.chain=2,n.adapt = 5000, n.burn = 15000, 
#'          n.iter = 150000, n.thin = 5)
#' jm.resid(jmfit.ew)
#' @export

jm.resid<-function(bayes.fit,posterior.mean=TRUE){
if(posterior.mean){
return(bayes.fit[["MCMC output"]][["cox.snell.residuals"]][,"posterior.mean"])
} else{
return(bayes.fit[["MCMC output"]][["cox.snell.residuals"]][,"posterior.median"])
}
}