From BUGS with BRugs to JAGS with rjags

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

Posterior predictive check can and should be Bayesian

Perception apparently employs prior knowledge that illumination comes from above, and that the surface itself has constant color.

Here's a new article, in press. I'm a bit reluctant about it because it got zero feedback and is unchanged from the initial submission, but, damn the torpedos ..., fools rush in ..., and all that sort of thing.

Abstract: Bayesian inference is conditional on the space of models assumed by the analyst. The posterior distribution indicates only which of the available parameter values are less bad than the others, without indicating whether the best available parameter values really fit the data well. A posterior predictive check is important to assess whether the posterior predictions of the least bad parameters are discrepant from the actual data in systematic ways. Gelman and Shalizi (2012a) assert that the posterior predictive check, whether done qualitatively or quantitatively, is non-Bayesian. I suggest that the qualitative posterior predictive check might be Bayesian, and the quantitative posterior predictive check should be Bayesian. In particular, I show that the “Bayesian p value,” from which an analyst attempts to reject a model without recourse to an alternative model, is ambiguous and inconclusive. Instead, the posterior predictive check, whether qualitative or quantitative, should be consummated with Bayesian estimation of an expanded model. The conclusion agrees with Gelman and Shalizi (2012a) regarding the importance of the posterior predictive check for breaking out of an initially assumed space of models. Philosophically, the conclusion allows the liberation to be completely Bayesian instead of relying on a non-Bayesian deus ex machina. Practically, the conclusion cautions against use of the Bayesian p value in favor of direct model expansion and Bayesian evaluation.

Kruschke, J. K. (in press). Posterior predictive check can and should be Bayesian: Comment on Gelman and Shalizi (2012a). British Journal of Mathematical and Statistical Psychology.
Get the full manuscript here.

One-group version of BEST (Bayesian estimation supersedes the t test)

The in-press article, Bayesian estimation supersedes the t test, focuses on the two-group case. Various readers have wanted a one-group version, which is now available. It is in the zip file with the two-group programs.

Here is an example of the code for using the program:

# Specify the data
y = c(101,100,102,104,102,97,105,105,98,101,100,123,105,103,100,95,102,106,
       109,102,82,102,100,102,102,101,102,102,103,103,97,97,103,101,97,104,
       96,103,124,101,101,100,101,101,104,100,101)

# Run the Bayesian analysis:
source("BEST1G.R")
mcmcChain = BEST1Gmcmc( y )

# Display the results:
BEST1Gplot( y , mcmcChain , compValm=100 , ROPEeff=c(-0.1,0.1) , pairsPlot=TRUE )

The function BEST1Gplot returns detailed numerical summaries of the posterior distribution (not shown here), and it also produces graphical output like this:

Workshop at Michigan State University, Sept. 14-15.

Workshop on doing Bayesian data analysis at Michigan State University, East Lansing, September 14-15. See details and registration info here.

A list of some previous workshops is here.

Gamma Likelihood Parameterized by MODE and SD

In a previous post I showed that it's more intuitive to think of a gamma distribution in terms of its mode and standard deviation (sd) than its mean and sd because the gamma distribution is typically skewed. But I did not show how to estimate the mode and sd in the context of a working JAGS/BUGS program, which I do here.

Suppose we have some data, y, that we wish to describe with a gamma distribution. (This requires that all the y values are non-negative, of course.) The dgamma function in JAGS/BUGS and R is parameterized by shape and rate parameters, not by mean, mode, or sd. Therefore we must reparameterize the shape and rate into equivalent mean, mode or sd. If we want to reparameterize by the mean of the gamma distribution, a JAGS/BUGS model statement could look like this:

  model {
    for ( i in 1:N ) {
      y[i] ~ dgamma( sh , ra )
    }
    # parameterized by mean (m) and standard deviation (sd)
    sh <- pow(m,2) / pow(sd,2)
    ra <-     m    / pow(sd,2)
    m ~ dunif(0,100)
    sd ~ dunif(0,100)
  }

If we want to reparameterize by the mode of the gamma distribution, a JAGS/BUGS model statement could look like this:

  model {
    for ( i in 1:N ) {
      y[i] ~ dgamma( sh , ra )
    }
    # parameterized by mode (m) and standard deviation (sd):
    sh <- 1 + m * ra
    ra <- ( m + sqrt( m^2 + 4*sd^2 ) ) / ( 2 * sd^2 )
    m ~ dunif(0,100)
    sd ~ dunif(0,100)
  }
In both model statements, the m parameter is given a dunif(0,100) prior, which implies that the priors in the two models are not equivalent because the mean is not the same as the mode. With vague priors the difference may be inconsequential for the posterior, but it is an important conceptual distinction. In fact, the whole point of reparameterizing is to make the prior more intuitive to specify! It's more intuitive, for most of us, to put a meaningful prior on the mode and sd than on the shape and rate.

Here are the results for the two models:

The left panels show the data (which are the same for both models) with a smattering of posterior predicted gamma distributions superimposed. The upper row shows the estimates for the model parameterized by the mean and sd, while the lower row shows the estimates for the model parameterized by mode and sd. You can see that the estimated mean (about 8.08) is larger than the estimated mode (about 7.01), but the estimated sd is essentially the same in the two models.


Appendix: Full program listed here, merely FYI.

graphics.off()
rm(list=ls(all=TRUE))
fileNameRoot="GammaModeJags" # for constructing output filenames
if ( .Platform$OS.type != "windows" ) {
  windows <- function( ... ) X11( ... )
}
require(rjags)         # Kruschke, J. K. (2011). Doing Bayesian Data Analysis:
                       # A Tutorial with R and BUGS. Academic Press / Elsevier.
#------------------------------------------------------------------------------
# THE MODEL.

mType = c("Mean","Mode")[1]

if ( mType == "Mean" ) {
  modelstring = "
  model {
    for ( i in 1:N ) {
      y[i] ~ dgamma( sh , ra )
    }
    # parameterized by mean (m) and standard deviation (sd)
    sh <- pow(m,2) / pow(sd,2)
    ra <-     m    / pow(sd,2)
    m ~ dunif(0,100)
    sd ~ dunif(0,100)
  }
  " # close quote for modelstring
} 

if ( mType == "Mode" ) {
  modelstring = "
  model {
    for ( i in 1:N ) {
      y[i] ~ dgamma( sh , ra )
    }
    # parameterized by mode (m) and standard deviation (sd):
    sh <- 1 + m * ra
    ra <- ( m + sqrt( m^2 + 4*sd^2 ) ) / ( 2 * sd^2 )
    m ~ dunif(0,100)
    sd ~ dunif(0,100)
  }
  " # close quote for modelstring
}
 
# Write model to a file, and send to BUGS:
writeLines(modelstring,con="model.txt")

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

# Load the data:

N = 87 ; m=8 ; sd=3
set.seed(47401)
y = rgamma( N , shape=m^2/sd^2 , rate=m/sd^2 )

# Specify the data as a list:
dataList = list( y = y , N = N )

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

initsList = list( m = mean(y) , sd = sd(y) )

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

parameters = c( "m" , "sd" ) 
adaptSteps = 1000              # 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).
nPerChain = 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=nPerChain , thin=thinSteps )
# resulting codaSamples object has these indices:
#   codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]

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

checkConvergence = F
if ( checkConvergence ) {
  windows()
  autocorr.plot( codaSamples , ask=F )
}

# 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 )
chainLength = NROW(mcmcChain)

source("plotPost.R")

windows(width=10,height=3)
layout( matrix(1:3,nrow=1) )
# Data with superimposed posterior predictive gamma's:
hist( y , xlab="y" , ylim=c(0,0.15) , main="Data" , col="grey" , prob=TRUE )
xComb=seq(min(y),max(y),length=501)
nPPC = 15
for ( idx in round(seq(1,chainLength,length=nPPC)) ) {
  if ( mType=="Mean" ) {
    m = mcmcChain[idx,"m"]
    sd = mcmcChain[idx,"sd"]
    sh = m^2 / sd^2
    ra = m   / sd^2
  }
  if ( mType=="Mode" ) {
    m = mcmcChain[idx,"m"]
    sd = mcmcChain[idx,"sd"]
    ra = ( m + sqrt( m^2 + 4*sd^2 ) ) / ( 2 * sd^2 )
    sh = 1 + m * ra
  }
  lines( xComb , dgamma( xComb , sh , ra ) , col="skyblue" , lwd=2 )
}
# Posterior estimates of parameters:
par( mar=c(3,1,2.5,0) , mgp=c(2,0.7,0) )
plotPost( mcmcChain[,"m"] , xlab=mType , main="Posterior Est." , showMode=T )
plotPost( mcmcChain[,"sd"] , xlab="sd" , main="Posterior Est." , showMode=T )
savePlot(file=paste(fileNameRoot,mType,sep=""),type="jpg")

Metropolis Algorithm: Discrete Position Probabilities

I was asked by a reader how I created Figure 7.2 of the book, reproduced at right, which shows the progression of discrete position probabilities at each step in a simple Metropolis algorithm. I couldn't find the original program, so I made a new one, reported here, with an additional example that illustrates a bimodal target distribution.

The main point is mentioned in the middle of p. 128: "At every time step, we multiply the current position probability vector w by the transition probability matrix T to get the position probabilities for the next time step. We keep multiplying by T over and over again to derive the long-run position probabilities. This process is exactly what generated the graphs in Figure 7.2."

The code consists of two main parts. First, build the transition matrix. Second, specify an initial position vector and repeatedly multiply by the transition matrix. Here is the code:


# For Figure 7.2, p. 123 of Doing Bayesian Data Analysis.
fileNameRoot = "DBDAfigure7pt2"

# Specify the target probability distribution.
pTargetName = "Linear" # for filename of saved graph
nSlots = 7
pTarget = 1:nSlots
pTarget = pTarget/sum(pTarget)

# Uncomment following lines for different example
#pTargetName = "Bimodal" # for filename of saved graph
#pTarget = c(1,2,2,1,3,3,1)
#nSlots = length(pTarget)
#pTarget = pTarget/sum(pTarget)

# Construct the matrix of proposal probabilities.
# Row is from position, column is to position.
proposalMatrix = matrix( 0 , nrow=nSlots , ncol=nSlots )
for ( fromIdx in 1:nSlots ) {
  for ( toIdx in 1:nSlots ) {
    if ( toIdx == fromIdx-1 ) { proposalMatrix[fromIdx,toIdx] = 0.5 }
    if ( toIdx == fromIdx+1 ) { proposalMatrix[fromIdx,toIdx] = 0.5 }
  }
}

# Construct the matrix of acceptance probabilities.
# Row is from position, column is to position.
acceptMatrix = matrix( 0 , nrow=nSlots , ncol=nSlots )
for ( fromIdx in 1:nSlots ) {
  for ( toIdx in 1:nSlots ) {
    acceptMatrix[fromIdx,toIdx] = min( pTarget[toIdx]/pTarget[fromIdx] , 1 )
  }
}

# Compute the matrix of overall move probabilities:
moveMatrix = proposalMatrix * acceptMatrix

# Compute the transition matrix, including the probability of staying in place:
transitionMatrix = moveMatrix
for ( diagIdx in 1:nSlots ) {
  transitionMatrix[diagIdx,diagIdx] = 1.0 - sum(moveMatrix[diagIdx,])
}
show( transitionMatrix )

# Specify starting position vector:
positionVec = rep(0,nSlots)
positionVec[round(nSlots/2)] = 1.0

# Loop through time steps and display position probabilities:
windows(height=7,width=7) # change to x11 for non-WindowsOS
layout( matrix(1:16,nrow=4) )
par( mar = c( 2.8 , 2.8 , 1.0 , 0.5 ) , mgp = c( 1.4 , 0.5 , 0 ) )
for ( timeIdx in 1:15 ) {
  # Display position probability:
  plot( 1:nSlots , positionVec , xlab="Position" , ylab="Probability" ,
        type="h" , lwd=5 , col="skyblue" )
  text( 1 , max(positionVec) , bquote(t==.(timeIdx)) , adj=c(0,1) )
  # Update next position probability:
  positionVec = positionVec %*% transitionMatrix
}
# Plot target distribution
plot( 1:nSlots , pTarget , xlab="Position" , ylab="Probability" ,
      type="h" , lwd=5 , col="skyblue" )
text( 1 , max(positionVec) , bquote(target) , adj=c(0,1) )

savePlot( filename=paste(fileNameRoot,pTargetName,sep="") , type="jpg" ) 

For the linearly increasing target distribution, the transition matrix looks like this:

     [,1]      [,2]      [,3]  [,4]      [,5]       [,6]      [,7]
[1,] 0.50 0.5000000 0.0000000 0.000 0.0000000 0.00000000 0.0000000
[2,] 0.25 0.2500000 0.5000000 0.000 0.0000000 0.00000000 0.0000000
[3,] 0.00 0.3333333 0.1666667 0.500 0.0000000 0.00000000 0.0000000
[4,] 0.00 0.0000000 0.3750000 0.125 0.5000000 0.00000000 0.0000000
[5,] 0.00 0.0000000 0.0000000 0.400 0.1000000 0.50000000 0.0000000
[6,] 0.00 0.0000000 0.0000000 0.000 0.4166667 0.08333333 0.5000000
[7,] 0.00 0.0000000 0.0000000 0.000 0.0000000 0.42857143 0.5714286

The row is the "from" position and the column is the "to" position. For example, the probability of moving from position 4 to position 3 is 0.375. 

The resulting graph looks like this:

If you uncomment the lines that specify a bimodal target, the resulting graph looks like this:

Program updates and changes

The plotPost() function, which plots histograms of MCMC chains annotated with HDI and other info, was changed a few weeks ago. plotPost() now returns summary information about the chain instead of summary information about the histogram.  A reader caught a conflict this caused in the program BernMetropolisTemplate.R, which has now been changed to accommodate the new plotPost(). The new programs can be found at the usual site: http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/Programs/?M=D.The zip file with all the programs has also been updated. Please let me know right away if you find other programs that have conflicts with the new plotPost(). Greatly appreciated!