jags: Run 'JAGS' from R

In R2jags: Using R to Run 'JAGS'

Run ‘JAGS’ from R

Description

The jags function takes data and starting values as input. It automatically writes a jags script, calls the model, and saves the simulations for easy access in R.

Usage

jags(data, inits, parameters.to.save, model.file="model.bug",
  n.chains=3, n.iter=2000, n.burnin=floor(n.iter/2),
  n.thin=max(1, floor((n.iter - n.burnin) / 1000)),
  DIC=TRUE, pD = FALSE, n.iter.pd = NULL, n.adapt = 100,
  working.directory=NULL, jags.seed = 123,
  refresh = n.iter/50, progress.bar = "text", digits=5,
  RNGname = c("Wichmann-Hill", "Marsaglia-Multicarry",
              "Super-Duper", "Mersenne-Twister"),
  jags.module = c("glm","dic"), quiet = FALSE,
  checkMissing = FALSE
  )

jags.parallel(data, inits, parameters.to.save, model.file = "model.bug",
             n.chains = 2, n.iter = 2000, n.burnin = floor(n.iter/2),
             n.thin = max(1, floor((n.iter - n.burnin)/1000)),
             n.cluster= n.chains, DIC = TRUE,
             working.directory = NULL, jags.seed = 123, digits=5,
             RNGname = c("Wichmann-Hill", "Marsaglia-Multicarry",
              "Super-Duper", "Mersenne-Twister"),
             jags.module = c("glm","dic"),
             export_obj_names=NULL,
             envir = .GlobalEnv
             )

jags2(data, inits, parameters.to.save, model.file="model.bug",
  n.chains=3, n.iter=2000, n.burnin=floor(n.iter/2),
  n.thin=max(1, floor((n.iter - n.burnin) / 1000)),
  DIC=TRUE, jags.path="",
  working.directory=NULL, clearWD=TRUE,
  refresh = n.iter/50)

Arguments

data	(1) a vector or list of the names of the data objects used by the model, (2) a (named) list of the data objects themselves, or (3) the name of a "dump" format file containing the data objects, which must end in ".txt", see example below for details.
inits	a list with n.chains elements; each element of the list is itself a list of starting values for the BUGS model, or a function creating (possibly random) initial values. If inits is NULL, JAGS will generate initial values for parameters.
parameters.to.save	character vector of the names of the parameters to save which should be monitored.
model.file	file containing the model written in BUGS code. Alternatively, as in R2WinBUGS, model.file can be an R function that contains a BUGS model that is written to a temporary model file (see tempfile) using write.model
n.chains	number of Markov chains (default: 3)
n.iter	number of total iterations per chain (including burn in; default: 2000)
n.burnin	length of burn in, i.e. number of iterations to discard at the beginning. Default is n.iter/2, that is, discarding the first half of the simulations. If n.burnin is 0, jags() will run 100 iterations for adaption.
n.cluster	number of clusters to use to run parallel chains. Default equals n.chains.
n.thin	thinning rate. Must be a positive integer. Set n.thin > 1 to save memory and computation time if n.iter is large. Default is max(1, floor(n.chains * (n.iter-n.burnin) / 1000)) which will only thin if there are at least 2000 simulations.
DIC	logical; if TRUE (default), compute deviance, pD, and DIC. The rule pD=var(deviance) / 2 is used.
pD	logical; if TRUE and DIC is also TRUE, then adds the computation of 'pD', using 'rjags::dic.samples()'. Defaults to FALSE.
n.iter.pd	number of iterations to feed 'rjags::dic.samples()' to compute 'pD'. Defaults at 1000.
n.adapt	number of iterations for which to run the adaptation, when creating the model object. Defaults at 100.
working.directory	sets working directory during execution of this function; This should be the directory where model file is.
jags.seed	random seed for JAGS, default is 123. This function is used for jags.parallell() and does not work for jags(). Use set.seed() instead if you want to produce identical result with jags()

.

jags.path	directory that contains the JAGS executable. The default is “”.
clearWD	indicating whether the files ‘data.txt’, ‘inits[1:n.chains].txt’, ‘codaIndex.txt’, ‘jagsscript.txt’, and ‘CODAchain[1:nchains].txt’ should be removed after jags has finished, default=TRUE.
refresh	refresh frequency for progress bar, default is n.iter/50
progress.bar	type of progress bar. Possible values are “text”, “gui”, and “none”. Type “text” is displayed on the R console. Type “gui” is a graphical progress bar in a new window. The progress bar is suppressed if progress.bar is “none”
digits	as in write.model in the R2WinBUGS package: number of significant digits used for BUGS input, see formatC. Only used if specifying a BUGS model as an R function.
RNGname	the name for random number generator used in JAGS. There are four RNGS supplied by the base moduale in JAGS: Wichmann-Hill, Marsaglia-Multicarry, Super-Duper, Mersenne-Twister
jags.module	the vector of jags modules to be loaded. Default are “glm” and “dic”. Input NULL if you don't want to load any jags module.
export_obj_names	character vector of objects to export to the clusters.
envir	default is .GlobalEnv
quiet	Logical, whether to suppress stdout in jags.model().
checkMissing	Default: FALSE. When TRUE, checks for missing data in categorical parameters and returns a sim.list with NA values if detected. It's recommended to supply jags() with complete data.

Details

To run:

1. Write a JAGS model in an ASCII file.

2. Go into R.

3. Prepare the inputs for the jags function and run it (see Example section).

4. The model will now run in JAGS. It might take awhile. You will see things happening in the R console.

Author(s)

Yu-Sung Su [e-mail], Masanao Yajima [e-mail], Gianluca Baio [e-mail]

References

Plummer, Martyn (2003) “JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling.” https://www.r-project.org/conferences/DSC-2003/Proceedings/Plummer.pdf.

Gelman, A., Carlin, J. B., Stern, H.S., Rubin, D.B. (2003) Bayesian Data Analysis, 2nd edition, CRC Press.

Sibylle Sturtz and Uwe Ligges and Andrew Gelman. (2005). “R2WinBUGS: A Package for Running WinBUGS from R.” Journal of Statistical Software 3 (12): 1–6.

Examples

  # An example model file is given in:
  model.file <- system.file(package="R2jags", "model", "schools.txt")
  # Let's take a look:
  file.show(model.file)
  # you can also write BUGS model as a R function, see below:

#=================#
# initialization  #
#=================#

  # data
  J <- 8.0
  y <- c(28.4,7.9,-2.8,6.8,-0.6,0.6,18.0,12.2)
  sd <- c(14.9,10.2,16.3,11.0,9.4,11.4,10.4,17.6)


  jags.data <- list("y","sd","J")
  jags.params <- c("mu","sigma","theta")
  jags.inits <- function(){
    list("mu"=rnorm(1),"sigma"=runif(1),"theta"=rnorm(J))
  }

  ## You can input data in 4 ways
  ## 1) data as list of character
  jagsfit <- jags(data=list("y","sd","J"), inits=jags.inits, jags.params,
                n.iter=10, model.file=model.file)

  ## 2) data as character vector of names
  jagsfit <- jags(data=c("y","sd","J"), inits=jags.inits, jags.params,
                n.iter=10, model.file=model.file)

  ## 3) data as named list
  jagsfit <- jags(data=list(y=y,sd=sd,J=J), inits=jags.inits, jags.params,
                n.iter=10, model.file=model.file)

  ## 4) data as a file
  fn <- "tmpbugsdata.txt"
  dump(c("y","sd","J"), file=fn)
  jagsfit <- jags(data=fn, inits=jags.inits, jags.params,
                  n.iter=10, model.file=model.file)
  unlink("tmpbugsdata.txt")

  ## You can write bugs model in R as a function

  schoolsmodel <- function() {
    for (j in 1:J){                     # J=8, the number of schools
      y[j] ~ dnorm (theta[j], tau.y[j]) # data model:  the likelihood
      tau.y[j] <- pow(sd[j], -2)        # tau = 1/sigma^2
    }
    for (j in 1:J){
      theta[j] ~ dnorm (mu, tau)        # hierarchical model for theta
    }
    tau <- pow(sigma, -2)               # tau = 1/sigma^2
    mu ~ dnorm (0.0, 1.0E-6)            # noninformative prior on mu
    sigma ~ dunif (0, 1000)             # noninformative prior on sigma
  }

  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
                n.iter=10, model.file=schoolsmodel)


#===============================#
# RUN jags and postprocessing   #
#===============================#
  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,
    n.iter=5000, model.file=model.file)

  # Can also compute the DIC using pD (=Dbar-Dhat), via dic.samples(), which
  # is a closer approximation to the original formulation of Spiegelhalter et
  # al (2002), instead of pV (=var(deviance)/2), which is the default in JAGS
   jagsfit.pD <- jags(data=jags.data, inits=jags.inits, jags.params,
     n.iter=5000, model.file=model.file, pD=TRUE)

  # Run jags parallely, no progress bar. R may be frozen for a while,
  # Be patient. Currenlty update afterward does not run parallelly
  #
   jagsfit.p <- jags.parallel(data=jags.data, inits=jags.inits, jags.params,
     n.iter=5000, model.file=model.file)
  
  # display the output
  print(jagsfit)
  plot(jagsfit)

  # traceplot
  traceplot(jagsfit.p)
  traceplot(jagsfit)

  # or to use some plots in coda
  # use as.mcmmc to convert rjags object into mcmc.list
  jagsfit.mcmc <- as.mcmc(jagsfit.p)
  jagsfit.mcmc <- as.mcmc(jagsfit)
  ## now we can use the plotting methods from coda
  #require(lattice)
  #xyplot(jagsfit.mcmc)
  #densityplot(jagsfit.mcmc)

  # if the model does not converge, update it!
  jagsfit.upd <- update(jagsfit, n.iter=100)
  print(jagsfit.upd)
  print(jagsfit.upd, intervals=c(0.025, 0.5, 0.975))
  plot(jagsfit.upd)

  # before update parallel jags object, do recompile it
  recompile(jagsfit.p)
  jagsfit.upd <- update(jagsfit.p, n.iter=100)



  # or auto update it until it converges! see ?autojags for details
  # recompile(jagsfit.p)
  jagsfit.upd <- autojags(jagsfit.p)
  jagsfit.upd <- autojags(jagsfit)

  # to get DIC or specify DIC=TRUE in jags() or do the following#
  dic.samples(jagsfit.upd$model, n.iter=1000, type="pD")

  # attach jags object into search path see "attach.bugs" for details
  attach.jags(jagsfit.upd)

  # this will show a 3-way array of the bugs.sim object, for example:
  mu

  # detach jags object into search path see "attach.bugs" for details
  detach.jags()

  # to pick up the last save session
  # for example, load("RWorkspace.Rdata")
  recompile(jagsfit)
  jagsfit.upd <- update(jagsfit, n.iter=100)

  recompile(jagsfit.p)
  jagsfit.upd <- update(jagsfit, n.iter=100)

#=============#
# using jags2 #
#=============#
  ## jags can be run and produces coda files, but cannot be updated once it's done
  ## You may need to edit "jags.path" to make this work,
  ## also you need a write access in the working directory:
  ## e.g. setwd("d:/")

  ## NOT RUN HERE
  ## Not run: 
    jagsfit <- jags2(data=jags.data, inits=jags.inits, jags.params,
      n.iter=5000, model.file=model.file)
    print(jagsfit)
    plot(jagsfit)
    # or to use some plots in coda
    # use as.mcmmc to convert rjags object into mcmc.list
    jagsfit.mcmc <- as.mcmc.list(jagsfit)
    traceplot(jagsfit.mcmc)
    #require(lattice)
    #xyplot(jagsfit.mcmc)
    #densityplot(jagsfit.mcmc)
  
## End(Not run)

R2jags documentation built on Oct. 13, 2024, 9:07 a.m.