6

I typically call JAGS from rjags with several chains for diagnostic purposes (e.g., 4 chains). After that I often want to do some post processing on the posterior parameter estimates (e.g., use predicted values, compute additional statistics). However, at this point it is a nuisance having the chains stored in a list.

What's a good way to combine the chains into a single parameter list?

Jeromy Anglim
  • 33,939
  • 30
  • 115
  • 173
  • you can combine the chains with `rbind` : `mcmcmtot <- rbind(mcmc[[1]], mcmc[[2]], mcmc[[3]])`and after work with this new object. For example : `pairs(mcmctot)`. – Philippe Jul 22 '14 at 08:12

2 Answers2

6

The runjags package has the function combine.mcmc. Its default setting is to combine one or more chains and return a single chain. E.g.,

library(runjags)
fit <- combine.mcmc(multichainfit)

It also has other options for combining chains.

Jeromy Anglim
  • 33,939
  • 30
  • 115
  • 173
3

Here is the rjags and coda solution. Say you generate two chains like this:

library(rjags)
x <- rnorm(100)
model.str <- 'model {for (i in 1:100) {
  x[i] ~ dnorm(mu, sd)}
  mu ~ dnorm(0, .0001)
  sd ~ dgamma(0.1, 0.1)}'
jags <- jags.model(textConnection(model.str), data = list(x = x),n.chains=2)
smpls <- coda.samples(model=jags,n.iter=2,variable.names=c("mu","sd"))
smpls
# [[1]]
# Markov Chain Monte Carlo (MCMC) output:
# Start = 1 
# End = 2 
# Thinning interval = 1 
#               mu        sd
# [1,] -0.09152588 0.9009973
# [2,]  0.05586651 1.0482184
# 
# [[2]]
# Markov Chain Monte Carlo (MCMC) output:
# Start = 1 
# End = 2 
# Thinning interval = 1 
#              mu        sd
# [1,] 0.06893182 0.7317955
# [2,] 0.13599206 0.8517304
# 
# attr(,"class")
# [1] "mcmc.list"

You can combine the two (or arbitrarily many) chains into a matrix by

do.call(rbind, smpls)
#               mu        sd
# [1,] -0.09152588 0.9009973
# [2,]  0.05586651 1.0482184
# [3,]  0.06893182 0.7317955
# [4,]  0.13599206 0.8517304

If you want an object of class mcmc, simply use the function mcmc: 

mcmc(do.call(rbind, smpls))

# [[1]]
# Markov Chain Monte Carlo (MCMC) output:
# Start = 1 
# End = 4 
# Thinning interval = 1 
#               mu        sd
# [1,] -0.09152588 0.9009973
# [2,]  0.05586651 1.0482184
# [3,]  0.06893182 0.7317955
# [4,]  0.13599206 0.8517304

If you want, you can preserve the start, end, or thin attribute of the original chain. For example to preserve thin:

mcmc(do.call(rbind, smpls), thin=thin(smpls))

(Of course you cannot preserve all three attributes)

sieste
  • 8,296
  • 3
  • 33
  • 48