JAGS with Stata

WinBUGS and OpenBUGS are just two of a growing number of blackbox programs for performing Bayesian analysis. Others include, JAGS (http://mcmc-jags.sourceforge.net/), Stan (http://mc-stan.org/) and BiiPS(https://alea.bordeaux.inria.fr/biips/doku.php). Of these, the program that is closest in style to WinBUGS and OpenBUGS is JAGS; it has a similar structure and it uses very similar samplers. So it should be easy to modify the wbs family of ado files so that JAGS can be called from Stata.

JAGS stands for Just Another Gibbs Sampler; this is a appropriate name because the program follows very much the same approach as the BUGS family of programs. The rationale for creating a new program was to provide open source code that could be used across many platforms and that would be easy for the user to extend, for example by adding their own samplers. JAGS is written in C++ and it can be downloaded from the web address mentioned above and it is free. The original version was written by Martyn Plummer and was released in 2007.

JAGS has several design features that make it particularly easy to use with R and there is  an R package called rjags that takes advantage of those features. My aim is to provide similar facilities for Stata users, as far as possible I will integrate them within the programs that already link Stata to WinBUGS and OpenBUGS.

So far, my experience is with using JAGS under Windows so I will restrict this posting to that implementation and return to discuss JAGS under Linux on another occasion. The cross-platform design should mean that the extension presents few problems. For the present, I will also limit myself to the base version of JAGS and not talk about writing one’s own extensions.

JAGS is designed for running via a script file in just the same way as one runs WinBUGS from Stata, because of this, JAGS has only a very basic/non-existent front-end and if the program is run interactively, one merely obtains a Windows terminal. However, this is not a problem for us as we will be running the program through Stata.

In the book ‘Bayesian analysis with Stata’ I describe how WinBUGS and OpenBUGS are run by providing four text files; a model file, a data file, an initial values file and a script file. JAGS is very similar except that these files have a slightly different format.

Since anyone who uses JAGS via Stata is likely to already have tried using WinBUGS or OpenBUGS with Stata, I will concentrate on the differences. More details can be found in the JAGS User Manual that is provided with the program when it is installed.

The Script File

The scripting language for JAGS is different to either WinBUGS or OpenBUGS but my objective is to hide that complexity from the Stata user by adding a jags option to wbsscript so that it will automatically writes the jags script file. The JAGS user manual describes the scripting language for anyone who is interested.

The Data & Initial Values File

One of the ways in which JAGS is geared to work with R is that JAGS formats its data files in the layout provided by R’s dump() function. Here is a piece of R code that demonstrates the format

> x <- c(1,2,3,4)

> p <- 5

> m <- matrix(c(1,2,3,4,5,6),2,3)

> m
         [,1] [,2] [,3]
         [1,] 1 3 5
         [2,] 2 4 6
> dump(c("x","p","m"),file="text.txt")
> file.show("text.txt")
         x <-
         c(1, 2, 3, 4)
         p <-
         5
         m <-
         structure(c(1, 2, 3, 4, 5, 6), .Dim = 2:3)

Results are indented and coloured green to distinguish them from the commands.

The first three lines define a vector, x, a scalar, p, and a matrix, m. The matrix is displayed and we see that, unlike WinBUGS, an R matrix is filled by column and not by row.

The dump() command sends the three data objects to a text file and R’s file.show() command displays the file in a text window. As this is not the format used by WinBUGS, wbslist needs to be re-written and unfortunately we cannot use the WinBUGS array format as implemented in wbsarray. JAGS cannot parse an R data frame.

In these dumped files, the variable names may optionally be included in quotes.

The Model File

There are small differences between the syntax of JAGS and OpenBUGS but they are too minor to bother us here. The most important of these differences are in the handling of censored data, otherwise an OpenBUGS model file is likely to run in JAGS without alteration.

Updating wbs

If you want to use JAGS with Stata then you will need to download new versions of some of the wbs programs. Eventually I will make these the versions that are supplied when you net install wbs into Stata, but for the moment I am treating them as experimental and under test. I have included the code for the new ado files in a pdf from where they can be copied and saved as ado files in the PLUS folder. As they are still being tested I have given them names such as wbsrun_v2.ado. In this way they do not replace the current release of the wbs programs but can sit alongside them. As yet there are no help files. The new ado files are listed in a downloadable pdf called  v2_adofiles . The security restrictions on this blog do not allow me to attach the ado files themselves.

Executables.txt

Once you have downloaded and installed JAGS, you will need to tell Stata were to find it on your computer. I do this by creating a text file called executables.txt in my PERSONAL ado folder that contains a line such as
JAGS,”c:/jags/jags-3.4.0/x64/bin/jags.bat”
The quotes contain the path to the Windows batch file jags.bat that is supplied when JAGS is installed.

An Example

Let us take a extremely simple example. Here is some Stata code that calls JAGS and estimates the mean, mu, of some simulated normally distributed data.

cd "C:\Temp"
 *-----------------------------------
 * Simulate the data & write to file
 *-----------------------------------
 clear
 set seed 390173
 set obs 5
 gen x = round(rnormal(),0.01)
 scalar n = 5
 wbslist_v2 (scalar n) (vect x) using data.txt , replace jags
 *-----------------------------------
 * Set initial value & write to file
 *-----------------------------------
 scalar mu = 1
 wbslist_v2 (scalar mu) using init.txt , replace jags
 *-----------------------------------
 * Write model file
 *-----------------------------------
 wbsmodel thisfile.do model.txt
 /*
 model {
    for( i in 1:n ) {
       x[i] ~ dnorm(mu,1)
    }
    mu ~ dnorm(1,1)
 }
 */
 *-----------------------------------
 * Write script file
 *-----------------------------------
 wbsscript_v2 using script.txt, replace ///
 model(model.txt) data(data.txt) init(init.txt) ///
 set(mu) coda(cd) burn(500) update(1000) jags
 *-----------------------------------
 * Run script & inspect log
 *-----------------------------------
 wbsrun_v2 using script.txt , jags
 type jagslog.txt
 *-----------------------------------
 * Read the data and plot
 *-----------------------------------
 wbscoda_v2 using cd, clear jags
 mcmcstats mu
 mcmctrace mu

Anyone who is familiar with running WinBUGS from Stata will recognise this pattern of code; the only substantial difference is that the wbs*_v2 commands allow a jags option.

However the text files that are created are different. Let’s look at each of them

data.txt
“n” <- 5.000
“x” <- c( 0.730, -0.050, 0.890, -0.860, 1.010)

init.txt
“mu” <- 1.000

model.txt
model {
     for( i in 1:n ) {
         x[i] ~ dnorm(mu,1)
     }
    mu ~ dnorm(1,1)
}

script.txt
model in ‘C:/Temp/model.txt’
data in ‘C:/Temp/data.txt’
compile
parameters in ‘C:/Temp/init.txt’
initialize
update 500
monitor mu
update 1000
coda * , stem(‘C:/Temp/cd’)
exit

As currently written the wbsrun_v2.ado file writes the JAGS log to the file jagslog.txt. Eventually I will add an option that allows the filename to be specified by the user. For our simple example the log reads,

Welcome to JAGS 3.4.0 on Mon Nov 10 15:27:22 2014
 JAGS is free software and comes with ABSOLUTELY NO WARRANTY
 Loading module: basemod: ok
 Loading module: bugs: ok
 Reading data file C:/Temp/data.txt
 Compiling model graph
 Resolving undeclared variables
 Allocating nodes
 Graph Size: 8
 Reading parameter file C:/Temp/init.txt
 Initializing model
 Updating 500
 Updating 1000

Once we have read the results into Stata, we can inspect them as we would results from WinBUGS by using the mcmc commands. Here is the trace plot of mu. The example is very simple but everything works as one would expect.

Next time I will take a slightly more complex model and try fitting it in WinBUGS, OpenBUGS and JAGS so that we can see if there are any noticeable differences in performance.