R/S3_methods.R

In kenkellner/jagsUI: A Wrapper Around 'rjags' to Streamline 'JAGS' Analyses

# Summary method
summary.jagsUI <- function(object, ...){
  object$summary
}

#Plot method
plot.jagsUI <- function(x, parameters=NULL, per_plot=4, ask=NULL, ...){

  if(is.null(ask))
    ask <- grDevices::dev.interactive(orNone = TRUE)
  plot_info <- get_plot_info(x, parameters, NULL, ask)
  dims <- c(min(length(plot_info$params), per_plot), 2)
  if(length(plot_info$params) <= per_plot)
    ask <- FALSE
  new_par <- list(mfrow = dims, mar = c(4,4,2.5,1), oma=c(0,0,0,0), ask=ask)
   
  #Handle par()
  old_par <- graphics::par(new_par)
  on.exit(graphics::par(old_par))  
  

  #Make plot
  for (i in plot_info$params){
    param_trace(x, i)
    param_density(x, i)
  }
}

# Print method
print.jagsUI <- function(x,digits=3,...){

  mc <- x$mcmc.info

  # Header
  #bugs.format=TRUE prints a nearly exact replica of WinBUGS-style output
  if(x$bugs.format){
    cat('Inference for Bugs model at \'',x$modfile,'\', fit using JAGS,','\n',sep="")
    cat(mc$n.chains,'chains, each with',mc$n.iter,'iterations (first ',mc$n.burnin,'discarded), n.thin =',mc$n.thin) 
    cat('\nn.sims =',mc$n.samples,'iterations saved','\n')  
  } else {
    cat('JAGS output for model \'',x$modfile,'\', generated by jagsUI.','\n',sep="")
    cat('Estimates based on',mc$n.chains,'chains of',mc$n.iter,'iterations,\n') 
    if(all(mc$sufficient.adapt)){
      cat('adaptation =',mean(mc$n.adapt),'iterations (sufficient),\n')
    } else{
      cat('adaptation =',mean(mc$n.adapt),'iterations (possibly insufficient),\n')
    }
    cat('burn-in = ',mc$n.burnin,' iterations and thin rate = ',mc$n.thin,',','\n',sep="")
    cat('yielding',mc$n.samples,'total samples from the joint posterior.','\n')
    if(!x$parallel){
      cat('MCMC ran for ',mc$elapsed.mins,' minutes at time ',paste(x$run.date),'.\n','\n',sep="")
    } else{
      cat('MCMC ran in parallel for ',mc$elapsed.mins,' minutes at time ',paste(x$run.date),'.\n','\n',sep="")
    }
  }

  # Table
  col_idx <- 1:ncol(x$summary)
  if(x$bugs.format){
    col_idx <- 1:9
    if(mc$n.chains == 1) col_idx <- 1:7
  } else {
    col_idx <- c(1,2,3,5,7,10,11,8,9)
    if(mc$n.chains == 1) col_idx <- c(1,2,3,5,7,10,11)
  }
  tab <- x$summary[,col_idx,drop=FALSE]
  tab <- as.data.frame(round(tab, digits))
  if(!x$bugs.format) tab[,"overlap0"] <- tab[,"overlap0"] == 1
  print(tab)

  #Print Rhat/n.eff information if necessary
  if(mc$n.chains>1){
    if(x$bugs.format){
      cat('\nFor each parameter, n.eff is a crude measure of effective sample size,','\n')
      cat('and Rhat is the potential scale reduction factor (at convergence, Rhat=1).','\n')
    } else {
      rhats <- unlist(tab[,"Rhat"])
      if(any(is.na(rhats))){
        cat('\n**WARNING** Some Rhat values could not be calculated.')
      }
      if(all(is.na(rhats))){
        cat('\n')
      } else if (max(rhats, na.rm=TRUE) > 1.1){
        cat('\n**WARNING** Rhat values indicate convergence failure.','\n')
      } else {
        cat('\nSuccessful convergence based on Rhat values (all < 1.1).','\n')
      }
      cat('Rhat is the potential scale reduction factor (at convergence, Rhat=1).','\n')
      cat('For each parameter, n.eff is a crude measure of effective sample size.','\n')
    }
  }

  #Print overlap0/f statistic info
  if(!x$bugs.format){
    cat('\noverlap0 checks if 0 falls in the parameter\'s 95% credible interval.\n')
    cat('f is the proportion of the posterior with the same sign as the mean;\n')
    cat('i.e., our confidence that the parameter is positive or negative.\n')
  }
  
  #Print DIC info
  if(x$calc.DIC & !is.null(x$pD)){
    if(x$bugs.format){
      cat('\nDIC info (using the rule, pD = var(deviance)/2)','\npD =',round(x$pD,1),'and DIC =',round(x$DIC,digits),'\n')
      cat('DIC is an estimate of expected predictive error (lower deviance is better).\n')
    } else {
      cat('\nDIC info: (pD = var(deviance)/2)','\npD =',round(x$pD,1),'and DIC =',round(x$DIC,digits),'\n')
      cat('DIC is an estimate of expected predictive error (lower is better).\n')
    }  
  }
    
    
}