From BUGS with BRugs to JAGS with rjags

As mentioned inseveral previous posts, I strongly recommend using JAGS instead of BUGS, and Ihave converted all the BUGS programs to JAGS versions. Here I provideguidelines for how to make the conversion in case you want to convert your ownprograms.
For a concreteexample, I will use the programs BernTwoBugs.R and BernTwoJags.R. I’ll proceedsection by section through the programs.
The header:
BUGS + BRugs version:
library(BRugs)
JAGS + rjags version:
require(rjags)
Instead of “require”it could say “library”. Also in the JAGS + rjags version I added a way for thegraphics to work on non-Windows machines. This is just R, so it can work withBUGS too:
if ( .Platform$OS.type != "windows" ) {   windows <- function(... ) X11( ... ) }

The model specification:
BUGS + BRugs version:
modelstring = "# BUGS model specification begins here...model {    # Likelihood. Eachflip is Bernoulli.     for ( i in 1 : N1 ) {y1[i] ~ dbern( theta1 ) }    for ( i in 1 : N2 ) {y2[i] ~ dbern( theta2 ) }    # Prior. Independentbeta distributions.    theta1 ~ dbeta( 3 , 3)    theta2 ~ dbeta( 3 , 3)}# ... end BUGS model specification" # close quote for modelstring# Write model to a file:.temp = file("model.txt","w") ;writeLines(modelstring,con=.temp) ; close(.temp) # Load model file into BRugs and check its syntax:modelCheck( "model.txt" )
JAGS + rjags version:
modelString = "# JAGS model specification begins here...model {    # Likelihood. Eachflip is Bernoulli.     for ( i in 1 : N1 ) {y1[i] ~ dbern( theta1 ) }    for ( i in 1 : N2 ) {y2[i] ~ dbern( theta2 ) }    # Prior. Independentbeta distributions.    theta1 ~ dbeta( 3 , 3)    theta2 ~ dbeta( 3 , 3)}# ... end JAGS model specification" # close quote for modelstring# Write the modelString to a file, using R commands:writeLines(modelString,con="model.txt")
Notice that themodel specification is the same in JAGS as in BUGS. Also, in both cases themodelString gets written to a file called “model.txt”. The JAGS + rjags versionuses a streamlined version of writeLines that would also work in the BUGSprogram, as it is just an R command. The only difference is in how thespecification gets communicated to BUGS or JAGS: BRugs uses the modelCheckcommand, but there is no analogous command in rjags.
The data:
BUGS + BRugs version:
# Specify the data in a form that is compatible with BRugsmodel, as a list:datalist = list(    N1 = 7 ,    y1 = c( 1,1,1,1,1,0,0) ,    N2 = 7 ,    y2 = c( 1,1,0,0,0,0,0))# Get the data into BRugs:modelData( bugsData( datalist ) )
JAGS + rjags version:
# Specify the data in a form that is compatible with JAGS model,as a list:dataList = list(    N1 = 7 ,    y1 = c( 1,1,1,1,1,0,0) ,    N2 = 7 ,    y2 = c( 1,1,0,0,0,0,0))
Thespecification of the data is the same in JAGS as in BUGS. The only differenceis in how the specification gets communicated to BUGS or JAGS.: BRugs uses the modelDatacommand, but there is no analogous command in rjags.
Initialize the chains:
BUGS + BRugs version:
modelCompile()modelGenInits()
JAGS + rjags version:
# Can be done automatically in jags.model() by commenting outinits argument.# Otherwise could be established as:# initsList = list( theta1 =sum(dataList$y1)/length(dataList$y1) , #                   theta2= sum(dataList$y2)/length(dataList$y2) )
The BUGSversion has to compile the model (using the BRugs modelCompile command) andthen generate initial values (using the BRugs modelGenInits command). The JAGSversion does not need any explicit initialization at this point, as thecommented code explains.

Run the chains:
BUGS + BRugs version:
samplesSet( c( "theta1" , "theta2" ) ) #Keep a record of sampled "theta" valueschainlength = 10000                     # Arbitrary length ofchain to generate.modelUpdate( chainlength )             # Actually generate the chain.
JAGS + rjags version:
parameters = c( "theta1" , "theta2" )     # The parameter(s) to be monitored.adaptSteps = 500             # Number of steps to "tune" the samplers.burnInSteps = 1000           # Number of steps to "burn-in" the samplers.nChains = 3                  # Number of chains to run.numSavedSteps=50000          # Total number of steps in chains to save.thinSteps=1                  # Number of steps to "thin" (1=keep every step).nIter = ceiling( ( numSavedSteps * thinSteps ) / nChains ) #Steps per chain.# Create, initialize, and adapt the model:jagsModel = jags.model( "model.txt" , data=dataList ,# inits=initsList ,                        n.chains=nChains , n.adapt=adaptSteps )# Burn-in:cat( "Burning in the MCMC chain...\n" )update( jagsModel , n.iter=burnInSteps )# The saved MCMC chain:cat( "Sampling final MCMC chain...\n" )codaSamples = coda.samples( jagsModel ,variable.names=parameters ,                            n.iter=nIter , thin=thinSteps )# resulting codaSamples object has these indices: #   codaSamples[[ chainIdx]][ stepIdx , paramIdx ]
Roughly theequivalent of BRugs modelCompile is rjags jags.model. For burning in, the roughequivalent of BRugs modelUpdate before samplesSet  is rjags update. Notice that the BUGS versionhere did no burning in. For the final chain, the rough equivalent of BRugsmodelUpdate is rjags coda.samples. Notice that rjags specifies the parametersto be stored with the variable.names argument in the coda.samples command,whereas BRugs specifies the parameters in its samplesSet command.
Examine the results:
BUGS + BRugs version:
theta1Sample = samplesSample( "theta1" ) # Put sampledvalues in a vector.theta2Sample = samplesSample( "theta2" ) # Put sampledvalues in a vector.
# Plot the trajectory of the last 500 sampled values.…
JAGS + rjags version:
# Convert coda-object codaSamples to matrix object for easierhandling.# But note that this concatenates the different chains into onelong chain.# Result is mcmcChain[ stepIdx , paramIdx ]mcmcChain = as.matrix( codaSamples )
theta1Sample = mcmcChain[,"theta1"] # Put sampledvalues in a vector.theta2Sample = mcmcChain[,"theta2"] # Put sampledvalues in a vector.
# Plot the trajectory of the last 500 sampled values.…
Once the chainis put into variables in R, they can be plotted the same way.
One bigdifference not shown above is how the chains can be examined forautocorrelation and convergence. My homegrown plotChains function uses BRugscommands and will not work with the JAGS output. Instead, JAGS+rjags returnsthe chain as a coda-package object, named codaSamples above. There are manyuseful functions in the coda package for displaying the chains, not shownabove. For example, summary( codaSamples ), plot( codaSamples ), and autocorr.plot(codaSamples ).