Title: | Using R to Run 'JAGS' |
---|---|
Description: | Providing wrapper functions to implement Bayesian analysis in JAGS. Some major features include monitoring convergence of a MCMC model using Rubin and Gelman Rhat statistics, automatically running a MCMC model till it converges, and implementing parallel processing of a MCMC model for multiple chains. |
Authors: | Yu-Sung Su [aut, cre] , Masanao Yajima [aut], Gianluca Baio [ctb] |
Maintainer: | Yu-Sung Su <[email protected]> |
License: | GPL (>= 2) |
Version: | 0.8-9 |
Built: | 2024-11-08 03:24:30 UTC |
Source: | https://github.com/suyusung/r2jags |
These are wraper functions for attach.bugs
and detach.bugs
, which
attach or detach three-way-simulation array of bugs object to the search path. See attach.all
for details.
attach.jags(x, overwrite = NA) detach.jags()
attach.jags(x, overwrite = NA) detach.jags()
x |
An |
overwrite |
If |
See attach.bugs
for details
Yu-Sung Su [email protected],
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.
# See the example in ?jags for the usage.
# See the example in ?jags for the usage.
The autojags
takes a rjags
object as input.
autojags
will update the model until it converges.
## S3 method for class 'rjags' update(object, n.iter=1000, n.thin=1, refresh=n.iter/50, progress.bar = "text", ...) autojags(object, n.iter=1000, n.thin=1, Rhat=1.1, n.update=2, refresh=n.iter/50, progress.bar = "text", ...)
## S3 method for class 'rjags' update(object, n.iter=1000, n.thin=1, refresh=n.iter/50, progress.bar = "text", ...) autojags(object, n.iter=1000, n.thin=1, Rhat=1.1, n.update=2, refresh=n.iter/50, progress.bar = "text", ...)
object |
an object of |
n.iter |
number of total iterations per chain, default=1000 |
n.thin |
thinning rate. Must be a positive integer, default=1 |
... |
further arguments pass to or from other methods. |
Rhat |
converegence criterion, default=1.1. |
n.update |
the max number of updates, default=2. |
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 |
Yu-Sung Su [email protected]
Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B. (2003): Bayesian Data Analysis, 2nd edition, CRC Press.
# see ?jags for an example.
# see ?jags for an example.
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.
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)
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)
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 |
pD |
logical; if |
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.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 |
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.
Yu-Sung Su [email protected], Masanao Yajima [email protected], Gianluca Baio [email protected]
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.
# 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)
# 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)
This function reads Markov Chain Monte Carlo output in the
CODA format produced by jags and returns an object of class
mcmc.list
for further output analysis using the
coda package.
jags2bugs(path=getwd(), parameters.to.save, n.chains=3, n.iter=2000, n.burnin=1000, n.thin=2, DIC=TRUE)
jags2bugs(path=getwd(), parameters.to.save, n.chains=3, n.iter=2000, n.burnin=1000, n.thin=2, DIC=TRUE)
path |
sets working directory during execution of this function; This should be the directory where CODA files are. |
parameters.to.save |
character vector of the names of the parameters to save which should be monitored. |
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.thin |
thinning rate, default is 2 |
DIC |
logical; if |
Yu-Sung Su [email protected], Masanao Yajima [email protected]
The recompile
takes a rjags
object as input.
recompile
will re-compile the previous saved rjags
object.
recompile(object, n.iter, refresh, progress.bar) ## S3 method for class 'rjags' recompile(object, n.iter=100, refresh=n.iter/50, progress.bar = "text")
recompile(object, n.iter, refresh, progress.bar) ## S3 method for class 'rjags' recompile(object, n.iter=100, refresh=n.iter/50, progress.bar = "text")
object |
an object of |
n.iter |
number of iteration for adapting, default is 100 |
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 |
Yu-Sung Su [email protected]
# see ?jags for an example.
# see ?jags for an example.
Displays a plot of iterations vs. sampled values for each variable in the chain, with a separate plot per variable.
traceplot(x, ...) ## S4 method for signature 'rjags' traceplot(x, mfrow = c(1, 1), varname = NULL, match.head = TRUE, ask = TRUE, col = rainbow( x$n.chains ), lty = 1, lwd = 1, ...)
traceplot(x, ...) ## S4 method for signature 'rjags' traceplot(x, mfrow = c(1, 1), varname = NULL, match.head = TRUE, ask = TRUE, col = rainbow( x$n.chains ), lty = 1, lwd = 1, ...)
x |
A bugs object |
mfrow |
graphical parameter (see |
varname |
vector of variable names to plot |
match.head |
matches the variable names by the beginning of the variable names in bugs object |
ask |
logical; if |
col |
graphical parameter (see |
lty |
graphical parameter (see |
lwd |
graphical parameter (see |
... |
further graphical parameters |
Masanao Yajima [email protected].
densplot
, plot.mcmc
,
traceplot