R/jagsUI.R

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

#' Call JAGS from R
#'
#' The \code{jagsUI} function is a basic user interface for running JAGS 
#'  analyses via package \code{rjags} inspired by similar packages like 
#'  \code{R2WinBUGS}, \code{R2OpenBUGS}, and \code{R2jags}. The user provides 
#'  a model file, data, initial values (optional), and parameters to save. 
#'  The function compiles the information and sends it to \code{JAGS}, 
#'  then consolidates and summarizes the MCMC output in an object of 
#'  class \code{jagsUI}.
#'
#' @param model_file Path to file containing the model written in 
#'    \code{BUGS} code.
#' @param data A named list of the data objects required by the model. 
#' @param inits A list with \code{n.chains} elements; each element of the
#'    list is itself a list of starting values for the \code{BUGS} model,
#'    \emph{or} a function creating (possibly random) initial values. If inits 
#'    is \code{NULL}, \code{JAGS} will generate initial values for parameters. 
#' @param params Character vector of the names of the parameters in the 
#'    model which should be monitored.
#' @param n_cores If > 1, the number of parallel CPU cores to use.
#'    MCMC chains are split between cores.
#' @param n_chains Number of Markov chains to run.
#' @param n_adapt Number of iterations to run in the adaptive phase.
#' @param n_iter Total number of iterations per chain (*including* warm-up).
#' @param n_warmup Number of iterations at the beginning of each chain to
#'    discard. Does not include the adaptive iterations.
#' @param n_thin Thinning rate. Must be a positive integer.
#' @param modules List of JAGS modules to load before analysis.
#' @param factories Optional character vector of factories to enable or 
#'    disable, in the format <factory> <type> <setting>. For example, to turn 
#'    \code{TemperedMix} on you would provide 
#'    \code{'mix::TemperedMix sampler TRUE'} (note spaces between parts). 
#'    Make sure you have the corresponding modules loaded as well.
#' @param quiet If \code{TRUE}, suppress console output.
#'
#' @return An object of class \code{jagsUI}. Notable elements in the output 
#'  object include:
#'  \item{samples}{The original output object from the \code{rjags} package, 
#'  as class \code{mcmc.list}.}
#'  \item{summary}{A summary of various statistics calculated based on model 
#'  output, in matrix form.}
#'  \item{model}{The \code{rjags} model object; this will contain multiple 
#'  elements if \code{n_cores > 1}.}
#'
#' @examples
#'  #Analyze Longley economic data in JAGS 
#'  #Package data
#'  data(longley)
#'  data <- list(gnp=longley$GNP, employed=longley$Employed, 
#'               n=length(longley$Employed))
#'
#'  #Write model file
#'  model_file <- tempfile()
#'  writeLines("
#'  model{
#'  #Likelihood
#'  for (i in 1:n){ 
#'    employed[i] ~ dnorm(mu[i], tau)     
#'    mu[i] <- alpha + beta*gnp[i]
#'  }
#'  #Priors
#'  alpha ~ dnorm(0, 0.00001)
#'  beta ~ dnorm(0, 0.00001)
#'  sigma ~ dunif(0,1000)
#'  tau <- pow(sigma,-2)
#'  }
#'  ", con=model_file)
#'
#' #Inits function
#'  inits <- function(){  
#'    list(alpha=rnorm(1,0,1),beta=rnorm(1,0,1),sigma=runif(1,0,3))  
#'  }
#'
#' #Parameters to save posteriors for
#' params <- c('alpha','beta','sigma')
#'
#' #Run analysis
#' out <- jagsUI(model_file, data, inits, params, 
#'             n_chains=3, n_iter=1000, n_warmup=500)
#'
#' #Take another 1000 posterior samples
#' out <- update(out, n_iter=1000)
#' 
#' @export

jagsUI <- function(model_file, data, inits=NULL, params, n_cores=1,
                 n_chains, n_adapt=100, n_iter, n_warmup=0, n_thin=1, 
                 modules=c("glm","dic"), factories=NULL, quiet=FALSE){

  #Process input
  #TODO: Don't combine everything like this
  run_dat <- mget(c("model_file", "data", "inits", "params", 
                    "n_cores", "n_chains", "n_adapt", "n_iter", "quiet", 
                    "n_warmup", "n_thin", "modules", "factories"))

  run_dat$inits <- process_inits(run_dat$inits, run_dat$n_chains, 
                                 run_dat$n_cores) 

  #Run MCMC
  if(n_cores == 1){    
    out <- run_model(run_dat)
  } else{   
    out <- run_parallel(run_dat)
  }

  #Calc summary stats
  out$summary <- calc_stats(out$samples)

  #Assign class
  class(out) <- 'jagsUI'

  return(out)

}