As mentioned in
several previous posts, I strongly recommend using JAGS instead of BUGS, and I
have converted all the BUGS programs to JAGS versions. Here I provide
guidelines for how to make the conversion in case you want to convert your own
programs.

For a concrete
example, I will use the programs BernTwoBugs.R and BernTwoJags.R. I’ll proceed
section 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 the
graphics to work on non-Windows machines. This is just R, so it can work with
BUGS too:

if ( .Platform$OS.type != "windows" ) {

windows <- function(
... ) X11( ... )

}

**The model specification:**

*BUGS + BRugs version:*

modelstring = "

# BUGS model specification begins here...

model {

# Likelihood. Each
flip is Bernoulli.

for ( i in 1 : N1 ) {
y1[i] ~ dbern( theta1 ) }

for ( i in 1 : N2 ) {
y2[i] ~ dbern( theta2 ) }

# Prior. Independent
beta 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. Each
flip is Bernoulli.

for ( i in 1 : N1 ) {
y1[i] ~ dbern( theta1 ) }

for ( i in 1 : N2 ) {
y2[i] ~ dbern( theta2 ) }

# Prior. Independent
beta 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 the
model specification is the same in JAGS as in BUGS. Also, in both cases the
modelString gets written to a file called “model.txt”. The JAGS + rjags version
uses a streamlined version of writeLines that would also work in the BUGS
program, as it is just an R command. The only difference is in how the
specification gets communicated to BUGS or JAGS: BRugs uses the modelCheck
command, but there is no analogous command in rjags.

**The data:**

*BUGS + BRugs version:*

# Specify the data in a form that is compatible with BRugs
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
)

)

# 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
)

)

The
specification of the data is the same in JAGS as in BUGS. The only difference
is in how the specification gets communicated to BUGS or JAGS.: BRugs uses the modelData
command, 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 out
inits argument.

# Otherwise could be established as:

# initsList = list( theta1 =
sum(dataList$y1)/length(dataList$y1) ,

# theta2
= sum(dataList$y2)/length(dataList$y2) )

The BUGS
version has to compile the model (using the BRugs modelCompile command) and
then generate initial values (using the BRugs modelGenInits command). The JAGS
version does not need any explicit initialization at this point, as the
commented code explains.

**Run the chains:**

*BUGS + BRugs version:*

samplesSet( c( "theta1" , "theta2" ) ) #
Keep a record of sampled "theta" values

chainlength = 10000 # Arbitrary length of
chain 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 the
equivalent of BRugs modelCompile is rjags jags.model. For burning in, the rough
equivalent of BRugs modelUpdate before samplesSet is rjags update. Notice that the BUGS version
here did no burning in. For the final chain, the rough equivalent of BRugs
modelUpdate is rjags coda.samples. Notice that rjags specifies the parameters
to 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 sampled
values in a vector.

theta2Sample = samplesSample( "theta2" ) # Put sampled
values in a vector.

# Plot the trajectory of the last 500 sampled values.

…

*JAGS + rjags version:*

# Convert coda-object codaSamples to matrix object for easier
handling.

# But note that this concatenates the different chains into one
long chain.

# Result is mcmcChain[ stepIdx , paramIdx ]

mcmcChain = as.matrix( codaSamples )

theta1Sample = mcmcChain[,"theta1"] # Put sampled
values in a vector.

theta2Sample = mcmcChain[,"theta2"] # Put sampled
values in a vector.

# Plot the trajectory of the last 500 sampled values.

…

Once the chain
is put into variables in R, they can be plotted the same way.

One big
difference not shown above is how the chains can be examined for
autocorrelation and convergence. My homegrown plotChains function uses BRugs
commands and will not work with the JAGS output. Instead, JAGS+rjags returns
the chain as a coda-package object, named codaSamples above. There are many
useful functions in the coda package for displaying the chains, not shown
above. For example, summary( codaSamples ), plot( codaSamples ), and autocorr.plot(
codaSamples ).

This is a very useful post. Thank you!

ReplyDeleteThanks, this is helpful. Much appreciated.

ReplyDeleteThanks for this post Dr. Kruschke. What brought me here is exercise 9.1A, where the program BernBetaMuKappaBugs.R won't run (theta[bleh]'s value is too low). I tried to increase the censoring, to no avail.

ReplyDeleteThanks again for this post. I'll try to convert to JAGS sooner rather than later.

The program BernBetaMuKappaJags.R should work. Just try running it as is (after installing JAGS, of course). Converting to JAGS is really pretty easy. Good luck and thanks for your interest.

ReplyDeleteThanks very much for this. Stumbled across it after being informed I need to convert to JAGS from WinBUGS. I'm a beginner Bayesian and this was very clear so thank you!

ReplyDelete