R/JMwaic.R

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

##################################
# WAIC
#####################################
#' Computes WAIC for Bayesian fit of the joint model
#' @description Returns WAIC.  
#' @keywords WAIC, Joint model, MCMC
#' @param bayes.fit a fit of the joint model using \code{jmreg.aft}.
#' @details WAIC is computed using -2*(lppd - pwaic), where
#'      lppd =  log pointwise predictive density, and
#'      pwaic = effective number of parameters. 
#'      The effective number of parameters is estimated using the 
#'      variance of individual terms in the log predictive density 
#'      summed over the data points. See  Gelman et al. (2014) for details.
#' @return Returns lppd, pwaic (effective number of parameters), and WAIC.
#' @references Gelman A, Hwang J, and Vehtari A, Understanding predictive information
#' criteria for bayesian models, Statistics and Computing 24: 997-1016, 2014.
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{jmreg.aft}}, \code{\link{jm.DIC}}, \code{\link{jm.summary}}
#' @examples
#'   data(pbc.long)
#'   data(pbc.surv)
#'   agec<-pbc.surv$age-mean(pbc.surv$age)
#'   pbc.surv0<-data.frame(pbc.surv,agec=agec)
#' # use natural splines
#'   lme.fit<-lme(log(bilirubin)~drug+ns(futime,2),data=pbc.long,
#'          random=~ns(futime,2)|id)
#'   surv.fit<-coxph(Surv(st,status2)~drug*agec,data=pbc.surv0,x=TRUE)
#' # use rand.model="ns"
#'   jmfit.w<-jmreg.aft(surv.fit=surv.fit,lme.fit=lme.fit,
#'           surv.model="weibull",rand.model="ns",timevar="futime",
#'           n.chain=2,n.adapt = 5000, n.burn = 15000, 
#'           n.iter = 150000, n.thin = 5)
#'   jm.WAIC(jmfit.w)
#' @export

jm.WAIC<-function(bayes.fit){
pwaic<-bayes.fit[["MCMC output"]]$pwaic
lppd<-bayes.fit[["MCMC output"]]$lppd
waic<-bayes.fit[["MCMC output"]]$waic
waic0<-cbind(lppd=lppd,pwaic=pwaic,WAIC=waic)
rownames(waic0)<-""
return(waic0)
}