R/JMsummary.R

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

######################################
#' Summary of a joint model fit
#' @description Produce a summary of the fit, including density, trace and autocorrelation plots. 
#' @keywords Joint model, MCMC, Posterior summary
#' @param bayes.fit fit of the joint model using \code{jmreg.aft}.
#' @param trace.plot produce trace plots if TRUE (default is FALSE).
#' @param density.plot produce density plots if TRUE (default is FALSE).
#' @param autocorr.plot produce autocorrelation plots if TRUE (default is FALSE).
#' @details The \code{coda} package (Plummer et al., 2006) is used to summarize and plot the output from 
#'     MCMC simulations. It also includes Gelman-Rubin statistic for convergence diagnostic.
#' @return Returns posterior summaries of the parameters, including posterior mean, median,
#'     standard deviation, 0.025 and 0.975 quantiles, and Gelman-Rubin statistic values. It also returns
#'     some basic information about the model, computational time, and MCMC simulation setup.
#' @references Plummer M, Best N, Cowles K, and Vines K, CODA: Convergence diagnosis and output 
#'        analysis for MCMC, R News, 6: 7-11, 2006. 
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{jm.DIC}}, \code{\link{jmreg.aft}}, \code{\link{jm.WAIC}}
#' @examples Consider the pbc data from package JM.
#'   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)
#'   jm.summary(jmfit.ew,trace.plot=TRUE,density.plot=TRUE)
#' @export


jm.summary<-function(bayes.fit,trace.plot=FALSE,density.plot=FALSE,autocorr.plot=FALSE){
surv.data<-bayes.fit[["MCMC output"]]$surv.data
surv.covariates<-colnames(surv.data)[-(1:3)]
p2<-length(surv.covariates)
long.data<-bayes.fit[["MCMC output"]]$long.data
long.covariates<-colnames(long.data)[-(1:2)]
p1<-length(long.covariates)
n.chain<-length(bayes.fit[["MCMC output"]]$MCMC)
chains<-list()
for(i in 1:n.chain){
chains[[i]]<-bayes.fit[["MCMC output"]]$MCMC[[i]]
}
sd.cor<-function(vv,col2){
v.mat<-solve(matrix(vv,ncol=sqrt(col2)))
v.cor0<-cov2cor(v.mat)
v.cor0[lower.tri(v.cor0,diag=T)] <- NA
v.cor1<-c(na.omit(c(t(v.cor0))))
v.sd<-sqrt(diag(v.mat))
return(c(v.cor1,v.sd))
}
L.coef<-list()
L.corr<-list()
L.sd<-list()
S.dist<-list()
S.coef<-list()
for(i in 1:n.chain){
L.coef0<-chains[[i]][,bayes.fit[[4]]$lcoef.loc]
colnames(L.coef0)<-c(long.covariates)
L.coef[[i]]<-L.coef0
L.sig1<-chains[[i]][,bayes.fit[[4]]$lcov.loc]
col2<-ncol(L.sig1)
L.cor.sd<-t(apply(L.sig1,1,sd.cor,col2=col2))
sd.ran<-L.cor.sd[,(ncol(L.cor.sd)-sqrt(col2)+1):ncol(L.cor.sd)]
L.cor2<-L.cor.sd[,1:(ncol(L.cor.sd)-sqrt(col2))]
if(is.vector(L.cor2)){L.cor2<-matrix(L.cor2,ncol=1)}
colnames(sd.ran)<-paste(rep("sd(b_",ncol(sd.ran)),c(0:(ncol(sd.ran)-1)),rep(")",ncol(sd.ran)),sep="")
sd.ran1<-cbind(sd.ran,sqrt(1/chains[[i]][,bayes.fit[[4]]$lres.loc]))
colnames(sd.ran1)[ncol(sd.ran1)]<-"sd(resid)"
L.sd[[i]]<-sd.ran1
vm<-sqrt(ncol(L.sig1))
mat <- matrix(0, vm, vm)
mat1<-which(mat ==0, arr.ind = T)
mat1[,1:2]<-mat1[,2:1]
mat2<-matrix(mat1[mat1[,"row"]<mat1[,"col"],],ncol=2)-1
gg<-NULL
for(k in 1:(vm*(vm-1)/2)){
gg<-c(gg,paste("corr",paste(mat2[k,], collapse = ''),collapse = ''))
}
colnames(L.cor2)<-gg
L.corr[[i]]<-L.cor2
S.dist0<-chains[[i]][,bayes.fit[[4]]$sdpar.loc]
S.coef0<-chains[[i]][,bayes.fit[[4]]$scoef.loc]
colnames(S.coef0)<-c(surv.covariates,"association")
S.dist[[i]]<-S.dist0
S.coef[[i]]<-S.coef0
}
S.coef.mcmc<-as.mcmc.list(lapply(S.coef,mcmc))
S.coef.sum<-summary(S.coef.mcmc)
R.surv.coef<-round(gelman.diag(S.coef.mcmc,multivariate=FALSE,autoburnin=FALSE)[[1]],digits=2)
stat.surv.coef<-round(cbind(mean=S.coef.sum[1]$statistics[,1],
       median=S.coef.sum[2]$quantiles[,3],
       SD=S.coef.sum[1]$statistics[,2],
       quantile.025=S.coef.sum[2]$quantiles[,1],
       quantile.975=S.coef.sum[2]$quantiles[,5],
       Rhat=R.surv.coef[,1],Rhat.upper.CI=R.surv.coef[,2]),digits=4)
S.dist.mcmc<-as.mcmc.list(lapply(S.dist,mcmc))
S.dist.sum<-summary(S.dist.mcmc)
R.surv.dist<-round(gelman.diag(S.dist.mcmc,multivariate=FALSE,autoburnin=FALSE)[[1]],digits=2)
stat.surv.dist<-round(cbind(mean=S.dist.sum[1]$statistics[,1],
       median=S.dist.sum[2]$quantiles[,3],
       SD=S.dist.sum[1]$statistics[,2],
       quantile.025=S.dist.sum[2]$quantiles[,1],
       quantile.975=S.dist.sum[2]$quantiles[,5],
       Rhat=R.surv.dist[,1],Rhat.upper.CI=R.surv.dist[,2]),digits=4)
L.coef.mcmc<-as.mcmc.list(lapply(L.coef,mcmc))
L.coef.sum<-summary(L.coef.mcmc)
R.long.coef<-round(gelman.diag(L.coef.mcmc,multivariate=FALSE,autoburnin=FALSE)[[1]],digits=2)
stat.long.coef<-round(cbind(mean=L.coef.sum[1]$statistics[,1],
       median=L.coef.sum[2]$quantiles[,3],
       SD=L.coef.sum[1]$statistics[,2],
       quantile.025=L.coef.sum[2]$quantiles[,1],
       quantile.975=L.coef.sum[2]$quantiles[,5],
       Rhat=R.long.coef[,1],Rhat.upper.CI=R.long.coef[,2]),digits=4)
L.corr.mcmc<-as.mcmc.list(lapply(L.corr,mcmc))
L.corr.sum<-summary(L.corr.mcmc)
R.long.corr<-round(gelman.diag(L.corr.mcmc,multivariate=FALSE,
  autoburnin=FALSE)[[1]],digits=2)
if(vm>2){
stat.long.corr<-round(cbind(mean=L.corr.sum[1]$statistics[,1],
       median=L.corr.sum[2]$quantiles[,3],
       SD=L.corr.sum[1]$statistics[,2],
       lower.95=L.corr.sum[2]$quantiles[,1],
       upper.95=L.corr.sum[2]$quantiles[,5],
       Rhat=R.long.corr[,1],Rhat.upper.CI=R.long.corr[,2]),digits=4)
} else{
stat.long.corr<-round(cbind(mean=L.corr.sum[1]$statistics[1],
       median=L.corr.sum[2]$quantiles[3],
       SD=L.corr.sum[1]$statistics[2],
       quantile.025=L.corr.sum[2]$quantiles[1],
       quantile.975=L.corr.sum[2]$quantiles[5],
       Rhat=R.long.corr[1],Rhat.upper.CI=R.long.corr[2]),digits=4)
}
rownames(stat.long.corr)<-gg
L.sd.mcmc<-as.mcmc.list(lapply(L.sd,mcmc))
L.sd.sum<-summary(L.sd.mcmc)
R.long.sd<-round(gelman.diag(L.sd.mcmc,multivariate=FALSE,autoburnin=FALSE)[[1]],digits=2)
stat.long.sd<-round(cbind(mean=L.sd.sum[1]$statistics[,1],
       median=L.sd.sum[2]$quantiles[,3],
       SD=L.sd.sum[1]$statistics[,2],
       quantile.025=L.sd.sum[2]$quantiles[,1],
       quantile.975=L.sd.sum[2]$quantiles[,5],
       Rhat=R.long.sd[,1],Rhat.upper.CI=R.long.sd[,2]),digits=4)
if(trace.plot){
plot(S.coef.mcmc,trace=T,density=F,ask = dev.interactive())
x11()
plot(S.dist.mcmc,trace=T,density=F,ask = dev.interactive())
x11()
plot(L.coef.mcmc,trace=T,density=F,ask = dev.interactive())
x11()
plot(L.sd.mcmc,trace=T,density=F,ask = dev.interactive())
x11()
plot(L.corr.mcmc,trace=T,density=F,ask = dev.interactive())
}
if(density.plot){
x11()
plot(S.coef.mcmc,trace=F,density=T,ask = dev.interactive())
x11()
plot(S.dist.mcmc,trace=F,density=T,ask = dev.interactive())
x11()
plot(L.coef.mcmc,trace=F,density=T,ask = dev.interactive())
x11()
plot(L.sd.mcmc,trace=F,density=T,ask = dev.interactive())
x11()
plot(L.corr.mcmc,trace=F,density=T,ask = dev.interactive())
}
if(autocorr.plot){
for(i in 1:n.chain){
x11()
autocorr.plot(S.coef.mcmc[[i]],ask = dev.interactive())
x11()
autocorr.plot(S.dist.mcmc[[i]],ask = dev.interactive())
x11()
autocorr.plot(L.coef.mcmc[[i]],ask = dev.interactive())
x11()
autocorr.plot(L.sd.mcmc[[i]],ask = dev.interactive())
x11()
autocorr.plot(L.corr.mcmc[[i]],ask = dev.interactive())
}
}
surv.dist<-bayes.fit[["MCMC output"]]$surv.model.full
if(bayes.fit[["MCMC output"]]$surv.model=="eweibull" || bayes.fit[["MCMC output"]]$surv.model=="ggamma"){
par.des1<-"rho = rate parameter"
par.des2<-"kappa and gamma = shape parameters"
} else{
par.des1<-"rho = rate parameter"
par.des2<-"kappa = shape parameter"
}
lme.fit<-bayes.fit[["MCMC output"]]$lme.fit
id <- as.vector(unclass(lme.fit$groups[[1]]))
surv.data<-bayes.fit[["MCMC output"]]$surv.data
cens<-sum(1-surv.data[,"status"])
cens.p<-round(cens*100/nrow(surv.data),digits=1)
cat(paste("Survival data:"),
    paste("          ","number of observations = ",nrow(surv.data)),
    paste("          ","censoring = ",cens,"(",cens.p,"%)"),
    paste("Longitudinal data:"), 
    paste("          ","number of observations = ",length(id)),
    paste("          ","number of individuals = ",length(unique(id))),paste(" "),sep="\n")
if(bayes.fit[["MCMC output"]]$lme.model=="ns"){
cat(paste("Survival model:",surv.dist), 
    paste("          ",par.des1),
    paste("          ",par.des2),
    paste("Longitudinal model:","LME"),
    paste("Association:",bayes.fit[[2]]),
    paste("Numerical integration:",bayes.fit[[6]]),paste(" "),sep="\n")
}
if(bayes.fit[["MCMC output"]]$lme.model=="simple"){
cat(paste("Survival model:",surv.dist), 
    paste("          ",par.des1),
    paste("          ",par.des2),
    paste("Longitudinal model:","LME"),
    paste("Association:",bayes.fit[[2]]),paste(" "),sep="\n")
}
chain.size<-nrow(bayes.fit[["MCMC output"]]$MCMC[[1]])
cat(paste("Number of chains =",bayes.fit[["MCMC output"]][6]), 
   paste("Iterations =",bayes.fit[["MCMC output"]][3],":",bayes.fit[["MCMC output"]][4]),
   paste("Thinning =",bayes.fit[["MCMC output"]][5]),
   paste("Sample size per chain =",chain.size),
   paste("Time =",bayes.fit$time),paste(""),sep="\n")
return(list("Distributional parameters (survival model)"=stat.surv.dist,
  "Regression coefficients (survival model)"=stat.surv.coef,
  "Regression coefficients (longitudinal model)"=stat.long.coef,
  "Variance components (longitudinal model)"=stat.long.sd,
  "Correlations (longitudinal model)"=stat.long.corr))
}