JAGS (Just Another Gibbs Sampler) is a general purpose MCMC engine similar to WinBUGS and OpenBUGS. I have a slight preference for JAGS as it is free and portable, works well on Linux, and interfaces well with R. It is tempting to write a tutorial introduction to JAGS and the corresponding R package, `rjags`, but there is a lot of material freely available on-line already, so it isn’t really necessary. If you are new to JAGS, I suggest starting with Getting Started with JAGS, rjags, and Bayesian Modelling. In this post I want to focus specifically on the problem of inlining JAGS models in R scripts as it can be very useful, and is usually skipped in introductory material.

#### JAGS and rjags on Ubuntu Linux

On recent versions of Ubuntu, assuming that R is already installed, the simplest way to install JAGS and `rjags` is using the command

sudo apt-get install jags r-cran-rjags

Now rjags is a CRAN package, so it can be installed in the usual way with `install.packages("rjags")`. However, taking JAGS and `rjags` direct from the Ubuntu repos should help to ensure that the versions of JAGS and `rjags` are in sync, which is a good thing.

#### Toy model

For this post, I will use a trivial toy example of inference for the mean and precision of a normal random sample. That is, we will assume data

with priors on and of the form

#### Separate model file

The usual way to fit this model in R using `rjags` is to first create a separate file containing the model

model { for (i in 1:n) { x[i]~dnorm(mu,tau) } mu~dnorm(cc,d) tau~dgamma(a,b) }

Then, supposing that this file is called `jags1.jags`, an R session to fit the model could be constructed as follows:

require(rjags) x=rnorm(15,25,2) data=list(x=x,n=length(x)) hyper=list(a=3,b=11,cc=10,d=1/100) init=list(mu=0,tau=1) model=jags.model("jags1.jags",data=append(data,hyper), inits=init) update(model,n.iter=100) output=coda.samples(model=model,variable.names=c("mu", "tau"), n.iter=10000, thin=1) print(summary(output)) plot(output)

This is all fine, and it can be very useful to have the model declared in a separate file, especially if the model is large and complex, and you might want to use it from outside R. However, very often for simple models it can be quite inconvenient to have the model separate from the R script which runs it. In particular, people often have issues with naming files correctly, making sure R is looking in the correct directory, moving the model with the R script, etc. So it would be nice to be able to just *inline* the JAGS model within an R script, to keep the model, the data, and the analysis all together in one place.

#### Using a temporary file

What we want to do is declare the JAGS model within a text string inside an R script and then somehow pass this into the call to `jags.model()`. The obvious way to do this is to write the string to a text file, and then pass the name of that text file into `jags.model()`. This works fine, but some care needs to be taken to make sure this works in a generic platform independent way. For example, you need to write to a file that you know doesn’t exist in a directory that is writable using a filename that is valid on the OS on which the script is being run. For this purpose R has an excellent little function called `tempfile()` which solves exactly this naming problem. It should always return the name of a file which does not exist in a writable directly within the standard temporary file location on the OS on which R is being run. This function is exceedingly useful for all kinds of things, but doesn’t seem to be very well known by newcomers to R. Using this we can construct a stand-alone R script to fit the model as follows:

require(rjags) x=rnorm(15,25,2) data=list(x=x,n=length(x)) hyper=list(a=3,b=11,cc=10,d=1/100) init=list(mu=0,tau=1) modelstring=" model { for (i in 1:n) { x[i]~dnorm(mu,tau) } mu~dnorm(cc,d) tau~dgamma(a,b) } " tmpf=tempfile() tmps=file(tmpf,"w") cat(modelstring,file=tmps) close(tmps) model=jags.model(tmpf,data=append(data,hyper), inits=init) update(model,n.iter=100) output=coda.samples(model=model,variable.names=c("mu", "tau"), n.iter=10000, thin=1) print(summary(output)) plot(output)

Now, although there is a file containing the model temporarily involved, the script is stand-alone and portable.

#### Using a text connection

The solution above works fine, but still involves writing a file to disk and reading it back in again, which is a bit pointless in this case. We can solve this by using another under-appreciated R function, `textConnection()`. Many R functions which take a file as an argument will work fine if instead passed a `textConnection` object, and the `rjags` function `jags.model()` is no exception. Here, instead of writing the model string to disk, we can turn it into a `textConnection` object and then pass that directly into `jags.model()` without ever actually writing the model file to disk. This is faster, neater and cleaner. An R session which takes this approach is given below.

require(rjags) x=rnorm(15,25,2) data=list(x=x,n=length(x)) hyper=list(a=3,b=11,cc=10,d=1/100) init=list(mu=0,tau=1) modelstring=" model { for (i in 1:n) { x[i]~dnorm(mu,tau) } mu~dnorm(cc,d) tau~dgamma(a,b) } " model=jags.model(textConnection(modelstring), data=append(data,hyper), inits=init) update(model,n.iter=100) output=coda.samples(model=model,variable.names=c("mu", "tau"), n.iter=10000, thin=1) print(summary(output)) plot(output)

This is my preferred way to use `rjags`. Note again that `textConnection` objects have many and varied uses and applications that have nothing to do with `rjags`.

Tags: Bayes, Bayesian, Gibbs, inline, jags, MCMC, model, R, rjags, rstats, sampler, tutorial

15/11/2012 at 21:48 |

An even nicer way (in my opinion) is to copy the write.model() function from R2WinBUGS …

20/11/2012 at 23:09 |

[...] a previous post I gave a quick introduction to using the rjags R package to access the JAGS Bayesian inference from [...]