Cautionary note on the use of Rhat alone to decide on sufficient chain length

[Marc Kery]

Dear all,

here is something that I just learned while translating code in the BPA book into NIMBLE and comparing the NIMBLE solution with that from JAGS (this is from BPA Section 13.3.2). When using 'small' MCMC settings (ni <- 10000  ;  nt <- 8  ;  nb <- 2000  ;  nc <- 3) I got quite discrepant solutions. Note though that the values of Rhat are quite respectable.

### Solution with NIMBLE

            mean    sd median    l95     u95  Rhat MCEpc

alpha.occ -0.711 0.321 -0.718 -1.326  -0.095 1.001 2.563

alpha.p    0.396 0.291  0.402 -0.194   0.957 0.988 4.093

beta.occ   4.794 1.152  4.681  2.772   7.079 1.006 2.956

beta.p    -2.398 0.480 -2.393 -3.319  -1.434 1.003 3.653

occ.fs    88.322 6.715 89.000 75.000 101.000 1.000 3.045

### Solution with JAGS

             mean     sd    2.5%     50%   97.5% overlap0     f  Rhat n.eff

alpha.occ  -0.698  0.289  -1.239  -0.704  -0.137    FALSE 0.986 1.002   975

beta.occ    3.141  0.737   1.881   3.067   4.763    FALSE 1.000 1.001  2306

alpha.p     0.185  0.278  -0.343   0.178   0.728     TRUE 0.748 1.001  2127

beta.p     -2.365  0.453  -3.309  -2.347  -1.526    FALSE 1.000 1.000  3000

occ.fs     80.520  7.935  65.000  81.000  95.000    FALSE 1.000 1.001  1435

When using 'larger' MCMC settings (ni <- 50000  ;  nt <- 40  ;  nb <- 10000  ;  nc <- 4), I get quite different solutions, which show a much better agreement between the two engines:

### Solution with NIMBLE

            mean    sd median    l95    u95  Rhat MCEpc

alpha.occ -1.080 0.256 -1.077 -1.572 -0.577 1.016 1.788

alpha.p    0.423 0.296  0.427 -0.162  1.010 0.996 1.448

beta.occ   3.277 0.569  3.242  2.236  4.431 1.003 1.585

beta.p    -1.761 0.435 -1.756 -2.640 -0.950 0.998 1.633

occ.fs    71.930 5.120 72.000 62.000 81.000 1.000 1.848

### Solution with JAGS

             mean     sd    2.5%     50%   97.5% overlap0    f  Rhat n.eff

alpha.occ  -1.071  0.259  -1.591  -1.068  -0.563    FALSE 1.00 1.000  4000

beta.occ    3.256  0.569   2.245   3.214   4.463    FALSE 1.00 1.001  2499

alpha.p     0.422  0.290  -0.128   0.421   1.002     TRUE 0.93 1.001  2471

beta.p     -1.759  0.431  -2.618  -1.757  -0.922    FALSE 1.00 1.000  4000

occ.fs     71.935  5.011  63.000  72.000  82.000    FALSE 1.00 1.000  3813

My take-home message is this: Rhat < 1.1 is not all, but ESS should not be ignored: it really pays to have large effective sample size.

Viele Gruesse  --- Marc

[falaf...@gmail.com]

Hi Marc,

This is a very good point. A low Rhat alone only tells you so much!

A few months ago I adapted some code (below) to run an "auto" and parallel version of NIMBLE that checks both for Rhat and ESS, for a case where I needed 5,000 independent samples of the posterior. 

Jeff

Jeffrey A. Hostetler, Ph.D.
Research Biologist
U.S. Geological Survey
Eastern Ecological Science Center
12100 Beech Forest Road
Laurel, MD 20708
USGS Staff Profile | ORCID

library(parallel)
library(MCMCvis)
library(nimble)

...

cl<-makeCluster(nc)

clusterExport(cl, c("nimblecode", "inits", "nimbledata", "nimbleconstants",
                    "params", "nt", "nb", "ni.it", "data.array"))
for (j in seq_along(cl)) {
  set.seed(j)
  init <- inits()
  clusterExport(cl[j], "init")
}

#run short chains in parallel on clusters
out1 <- clusterEvalQ(cl, {
  library(nimble)
  library(coda)
  model <- nimbleModel(code = nimblecode, name = "model",
                       constants = nimbleconstants, data = nimbledata,
                       inits = init)
  Cmodel <- compileNimble(model)
  modelConf <- configureMCMC(model)
  modelConf$resetMonitors()
  modelConf$addMonitors(params)
  modelMCMC <- buildMCMC(modelConf)
  CmodelMCMC <- compileNimble(modelMCMC, project = model)
  CmodelMCMC$run(nb, reset = TRUE, thin=nt)
  return(as.mcmc(as.matrix(CmodelMCMC$mvSamples)))
})

out.mcmc <- mcmc.list(out1)
out.mcmc <- window(out.mcmc, start=2) #remove first row of NA's
out.summary <- MCMCsummary(out.mcmc)

# extend chains if not converged
start.time = Sys.time()  # Set timers
place <- nb
while ((max(out.summary$Rhat, na.rm = T)>1.1 ||
        min(out.summary$n.eff[out.summary$n.eff>0], na.rm = T)<5000) && place<niter.max) {
  cat(place, round(max(out.summary$Rhat, na.rm = T), 2),
      min(out.summary$n.eff[out.summary$n.eff>0], na.rm = T), "\n")
  out2 <- clusterEvalQ(cl, {
    CmodelMCMC$run(ni.it, reset = FALSE, thin=nt)
    return(as.mcmc(as.matrix(CmodelMCMC$mvSamples)))
  })
 
  out.mcmc.update1 <- mcmc.list(out2)
  #Remove burn-in
  out.mcmc.update1 <- window(out.mcmc.update1, start=place + 1)
 
  # Save summary
  options(max.print = 5000)
  out.summary <- MCMCsummary(out.mcmc.update1)
  sink(file=paste0(name2, ".summary.txt"))
  print(out.summary)
  sink()
  place <- place + ni.it
}

cat(place, round(max(out.summary$Rhat, na.rm = T), 2),

    min(out.summary$n.eff[out.summary$n.eff>0], na.rm = T), "\n")

stopCluster(cl)

[Marc Kery]

Dear Jeff,

thanks a lot ! Will try out sometimes (I am still mostly a JAGS user ;)...). This is surely very useful. We're about to distribute a NIMBLE translation of the code for the BPA book. For this, we were privileged to have Ken Kellner write us code that produces output for NIMBLE almost like what jagsUI does (i.e., including n.eff and Rhat).

Best regards  ---- Marc

---

From: hmec...@googlegroups.com <hmec...@googlegroups.com> on behalf of falaf...@gmail.com <falaf...@gmail.com>
Sent: Tuesday, June 24, 2025 21:07
To: hmecology: Hierarchical Modeling in Ecology <hmec...@googlegroups.com>
Subject: Re: Cautionary note on the use of Rhat alone to decide on sufficient chain length

 

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/hmecology/63b159c7-0452-40a5-8e7e-e6af2f814ca5n%40googlegroups.com.