Another reason to use JAGS instead of BUGS

BUGS is the pioneering software that made MCMC available to so many of us, but it has some problems with robustness that are not suffered by the subsequent software JAGS. Readers of DBDA tell me of some new problems running models in BUGS, which I have not been able to solve in BUGS, but the same model runs fine in JAGS. The purpose of this post is to point out the particular problem and the complete solution in JAGS. (If you have a solution for BUGS, please comment.)

I converted all the BUGS programs to JAGS versions back on 01 January 2012. Please see the previous post for information about JAGS.

A reader reports this problem:

I [the reader] have received this error several times:
"update error for node <theta[66]> algorithm conjugate beta updater error argument two  has too small a value "
...  In particular, I [the reader] noticed for exercise 12.3, if I remove the 0's from the data, I do not get the error. 

When I try running the program from Exercise 12.3 with BUGS, I also get the error, although I did not get any error a few years ago when I created the program and reported its results in the solutions manual. I have tried to solve the problem by arbitrarily censoring/truncating the distributions in BUGS (as I did in some other programs in the book), but without success in this application.

I quickly re-wrote the BUGS+BRugs program for JAGS+rjags, going step-by-step through the guidelines from this previous post. The resulting program, using the identical model statement, runs fine in JAGS and produces results identical to those in the solutions manual. Here is the complete program and results:


graphics.off()
rm(list=ls(all=TRUE))
fileNameRoot="Exercise12.3.A.JAGS" # for constructing output filenames

require(rjags)

if ( .Platform$OS.type != "windows" ) {
  windows <- function( ... ) X11( ... )
}


#------------------------------------------------------------------------------
# THE MODEL.

modelstring = "
model {
   for ( subjIdx in 1:nSubj ) {
      # Likelihood:
      z[subjIdx] ~ dbin( theta[subjIdx] , N[subjIdx] )
      # Prior on theta: Notice nested indexing.
      theta[subjIdx] ~ dbeta( a[cond[subjIdx]] , b[cond[subjIdx]] )
   }
   for ( condIdx in 1:nCond ) {
      a[condIdx] <- mu[condIdx] * kappa[condIdx]
      b[condIdx] <- (1-mu[condIdx]) * kappa[condIdx]
      # Hyperprior on mu and kappa:
      mu[condIdx] ~ dbeta( Amu , Bmu )
      kappa[condIdx] ~ dgamma( Skappa , Rkappa )
   }
   # Constants for hyperprior:
   Amu <- 0.35 * 2
   Bmu <- 0.65 * 2
   Skappa <- pow(meanGamma,2)/pow(sdGamma,2)
   Rkappa <- meanGamma/pow(sdGamma,2)
   meanGamma <- 10
   sdGamma <- 10
}
" # close quote for modelstring
# Write model to a file:
writeLines(modelstring,con="model.txt")

#------------------------------------------------------------------------------
# THE DATA.

nCorrOfSubj = c( 8, 4, 6, 3, 1, 4, 4, 6, 4, 2, 2, 1, 1, 4, 3 ,3 ,2 ,6, 3, 4,
                 2, 1, 1, 3, 2, 7, 2, 1, 3, 1, 0, 2, 4, 2, 3, 3, 0, 1, 2, 2 )
CondOfSubj = c( rep(1,20) , rep(2,20) )
nTrlOfSubj = rep(10,40)
nSubj = length(CondOfSubj)
nCond = length(unique(CondOfSubj))

# Specify the data in a form that is compatible with BRugs model, as a list:
dataList = list(
 nCond = nCond ,
 nSubj = nSubj ,
 cond = CondOfSubj ,
 N = nTrlOfSubj ,
 z = nCorrOfSubj
)

#------------------------------------------------------------------------------
# INTIALIZE THE CHAINS.


#------------------------------------------------------------------------------
# RUN THE CHAINS.

parameters = c("mu","kappa","theta")     # 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 ]


#------------------------------------------------------------------------------
# EXAMINE THE RESULTS.

source("plotPost.R")

mcmcChain = as.matrix( codaSamples )

# Extract parameter values and save them.
mu = NULL
kappa = NULL
for ( condIdx in 1:nCond ) {
   mu = rbind( mu , mcmcChain[, paste("mu[",condIdx,"]",sep="") ] )
   kappa = rbind( kappa , mcmcChain[, paste("kappa[",condIdx,"]",sep="") ] )
}
save( mu , kappa , file=paste(fileNameRoot,"MuKappa.Rdata",sep="") )
chainLength = NCOL(mu)

# Histograms of mu differences:
windows(width=7,height=7)
plotPost( mu[1,]-mu[2,] , xlab=expression(mu[1]-mu[2]) , main="" ,
          breaks=20 , compVal=0 , ROPE=c(-.1,.1) )
savePlot( file=paste(fileNameRoot,"MuDiffs.jpg",sep="") , type="jpg" )

# Scatterplot of mu,kappa in each condition:
windows(width=7,height=7)
muLim = c(min(mu),max(mu))
kappaLim = c(min(kappa),max(kappa))
mainLab="Posterior"
thindex = round( seq( 1 , chainLength , len=300 ) )
plot( mu[1,thindex] , kappa[1,thindex] , main=mainLab ,
      xlab=expression(mu[c]) , ylab=expression(kappa[c]) , cex.lab=1.75 ,
      xlim=muLim , ylim=kappaLim , log="y" , col="red" , pch="1" )
points( mu[2,thindex] , kappa[2,thindex] , col="blue" , pch="2" )
savePlot( file=paste(fileNameRoot,"Scatter.jpg",sep="") , type="jpg" )