jags {R2jags} | R Documentation |
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, 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 |
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 |
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.cluster |
number of clusters to use to run parallel chains. Default equals n.chains. |
n.thin |
thinning rate. Must be a positive integer. Set
|
DIC |
logical; if |
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.path |
directory that contains the |
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 |
refresh |
refresh frequency for progress bar, default is |
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 |
digits |
as in |
RNGname |
the name for random number generator used in JAGS. There are four RNGS
supplied by the base moduale in JAGS: |
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 |
checkMissing |
Default: FALSE. When TRUE, checks for missing data in categorical parameters
and returns a |
Details
To run:
Write a
JAGS
model in an ASCII file.Go into R.
Prepare the inputs for the
jags
function and run it (see Example section).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 suyusung@tsinghua.edu.cn, Masanao Yajima yajima@bu.edu
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)
# 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)