R/JMsurv.R

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

################################################################
########################################################################
#' Dynamic predictions of survival probabilites 
#' @description Dynamic predictions of survival probabilites \code{P(T >= t | T >= s)} for \code{t >= s},
#'         given longitudinal measurements, covariate information, and survival up to time \code{s}. 
#' @keywords AFT model, Joint moldeing, MCMC, Prediction, Survival probabilities
#' @param bayes.fit fit of the joint model using \code{jmreg.aft}.
#' @param newsurvdata a data frame that contains covariate information for the time-to-event process.
#'         Names of the variables in \code{newsurvdata} must be the same as in the data frame that was used 
#'         to fit the survival model using \code{coxph}.
#' @param newlongdata a data frame that contains repeated measurements and covariate information for the longitudinal process.
#'      Names of the variables in \code{newlongdata} must be the same as in the data frame that was used to fit 
#'      the linear mixed effects model using \code{lme}.
#' @param st a numeric vector of times at which survival probabilities \code{P(T >= t | T >= s)} are to be computed.
#'        The first element must be \code{s}. If \code{NULL}, \code{s} is taken to be the last time at which the 
#'        subject provided a longitudinal measurement, and \code{t} is \code{seq(s,max.st,length=15)}, where
#'        \code{max.st} is \code{floor(max(event times))}.
#' @param control a list of control values (see Details).
#' @details For estimation, we consider a Monte Carlo approach similar to the algorithm proposed by Rizopoulos (2016)
#'   to construct a subject's profile (i.e., drawing a sample from the posterior distribution of \eqn{b_i}),
#'   given a realization of the population parameters from the MCMC samples of the fitted joint model and 
#'   survival up to time \code{s}. This is implemented in JAGS (Plummer, 2017).
#'
#' \code{control}: a list of control values with components:
#'     \itemize{ 
#'      \item \code{n.sample}: size of the MCMC sample to be used to derive the posterior summeries of the conditional 
#'           probabilities (default is 200).
#'      \item \code{n.adapt}: number of steps for adaptation to draw \eqn{b_i} 
#'              (random effects), given a realization of the population parameters 
#'              from the MCMC samples of the fitted joint model and urvival up to time \code{s} (default is 1000).
#'      \item \code{n.update}: burn-in iterations to draw \eqn{b_i} 
#'              (random effects), given a realization of the population parameters 
#'              from the MCMC samples of the fitted joint model and urvival up to time \code{s} (default is 5000).
#'          }
#' @return Posterior summeries of the survival probabilities.
#' @references Plummer M, JAGS Version 4.3.0 user manual, 2017, 
#'           \url{http://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf}. 
#' @references Rizopoulos D, The R package JMbayes for fitting joint models for longitudinal and time-to-event data Using MCMC, 
#'           Journal of Statistical Software, 72(7): 1-45, 2016.
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{jmreg.aft}}, \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.gg<-jmreg.aft(surv.fit=surv.fit,lme.fit=lme.fit,
#'          surv.model="ggamma",rand.model="simple",timevar="obstime",
#'          n.chain=2,n.adapt = 5000, n.burn = 15000, 
#'          n.iter = 150000, n.thin = 5)
#'   newsurvdata<-data.frame(aids.id[aids.id$patient==409,])
#'   newlongdata<-data.frame(aids[aids$patient==409,])
#'   surv.id0<-jm.surv(bayes.fit=jmfit.gg,newsurvdata=newsurvdata,newlongdata=newlongdata)
#'   newsurvdata<-data.frame(aids.id[aids.id$patient==404,])
#'   newlongdata<-data.frame(aids[aids$patient==404,])
#'   surv.id1<-jm.surv(bayes.fit=jmfit.gg,newsurvdata=newsurvdata,newlongdata=newlongdata)
#'   plot(surv.id1[[1]][,1],surv.id1[[1]][,2],ylim=c(0.38,1),
#'        type="l",lty=2,xlab="Time (t)", ylab="P(T > t | T > 12)")
#'   lines(surv.id0[[1]][,1],surv.id0[[1]][,2],lty=1)
#'   legend("topright",legend=c("Patient id # 409","Patient id # 404"),lty=c(1,2))
#' @export

jm.surv<-function(bayes.fit,newsurvdata,newlongdata,st=NULL,control=list()){
con<-list(n.sample=200,seed=sample(123456,1),n.adapt=1000,n.update =5000)
nmsC <- names(con)
con[(namc <- names(control))] <- control
if (length(noNms <- namc[!namc %in% nmsC])) 
  stop("unknown names in control: ", paste(noNms, collapse = ", "))
if(is.null(bayes.fit[["MCMC output"]][["quad.points"]])){
if(bayes.fit[["MCMC output"]]$surv.model=="eweibull"){
sfit<-jm.surv.ew.simple(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
if(bayes.fit[["MCMC output"]]$surv.model=="ggamma"){
sfit<-jm.surv.gg.simple(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
if(bayes.fit[["MCMC output"]]$surv.model=="weibull"){
sfit<-jm.surv.w.simple(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
if(bayes.fit[["MCMC output"]]$surv.model=="llogistic"){
sfit<-jm.surv.ll.simple(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
if(bayes.fit[["MCMC output"]]$surv.model=="lnormal"){
sfit<-jm.surv.ln.simple(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
} else{
if(bayes.fit[["MCMC output"]]$surv.model=="eweibull"){
sfit<-jm.surv.ew(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
if(bayes.fit[["MCMC output"]]$surv.model=="ggamma"){
sfit<-jm.surv.gg(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
if(bayes.fit[["MCMC output"]]$surv.model=="weibull"){
sfit<-jm.surv.w(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
if(bayes.fit[["MCMC output"]]$surv.model=="llogistic"){
sfit<-jm.surv.ll(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
if(bayes.fit[["MCMC output"]]$surv.model=="lnormal"){
sfit<-jm.surv.ln(bayes.fit=bayes.fit,newsurvdata=newsurvdata,newlongdata=newlongdata,st=st,
   n.sample=con$n.sample,n.adapt=con$n.adapt,
   n.update=con$n.update,seed=con$seed)
}
}
return(sfit)
}