MCMC diagnostics problem and proposed fix

[adsmithca]

Hello, I found a problem on line 160 of of your fit_jabba() function here: https://github.com/jabbamodel/JABBA/blob/master/R/fit_jabba.R

I found it while trying to see if I could modify your code to produce bilingual (English and French) figures and also use the plotting functions from bayesplot and shinystan on JABBA output by saving the mod object created by R2jags::jags() on line 139 in the output of fit_jabba().

The line in question is:
posteriors = mod$BUGSoutput$sims.list

This code concatenates the chains and randomizes the order of the draws. This is because of a bug in R2jags and is discussed here by, among others, Martyn Plummer, the author of rjags etc. https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/cc61b820/. Essentially .$sims.list and .$sims.matrix do not maintain draw order. The author of the R2jags package subsequently modified the R2jags:::as.mcmc.jags() function to maintain the correct order of the draws.

As a result, the model convergence diagnostics from coda, coda::geweke.diag() and coda::heidel.diag(), as well as visual inspection of the mcmc chains from jbplot_mcmc() cannot be used as criteria for model convergence. However, this bug does not change parameter estimates so no worries there.

This can easily be fixed.
Line 160 can be replaced with :

mod_mcmc <- R2jags:::as.mcmc.rjags(mod)
modmcmc <- as.data.frame(as.matrix(mod_mcmc))

These objects can then be saved along with all the other fit_jabba output. For example, with these lines of code after line 432:

jabba$jags <- mod
jabba$mcmc <- mod_mcmc

From this, the model convergence diagnostics tools you provide can be used internally or externally. In addition, the saved mcmc and fitted jags model could subsequently be used by other packages such as tidybayes, bayesplot, shinystan, etc.

I also have another proposal for your code. It could eventually be included as an argument in the function that be turned on or off. The code essentially updates the jags model until Gelman and Rubin's Rhat diagnostic for all estimated parameters is below 1.1. Note that I use the pipe operator %>% from magrittr.

 # Andrew added this next bit ----
  # Effective Sample (ESS) should be as large as possible, although for most applications,
  # an effective sample size greater than 1000 is sufficient for stable estimates (Bürkner, 2017).
  # The ESS corresponds to the number of independent samples with the same estimation power as the N auto-correlated samples.
  # It is is a measure of “how much independent information there is in auto-correlated chains” (Kruschke 2015, p182-3).

  # get model summary and examine Rhat, DIC, effective sample size, and Gelman Rubin diagnostic
  cat("Initial JAGS run summary: \n")

  cat("\n DIC = ", mod$BUGSoutput$DIC, "\n")

  mod_summary <- mod$BUGSoutput$summary %>% as.data.frame()
summary(mod_summary)

  cat("\n Checking convergence via the Gelman-Rubin statistic (potential scale reduction) comparing within and between chain variance ... \n")
  if (any(mod_summary$Rhat > 1.1)) {
    message("\n Rhat values for one or more of the estimated values are above 1.1 indicating the model has not converged ... updating model with more iterations and treating initial as burn in \n")
    message("\n Mean effective sample size among parameters = ", round(mean(mod_summary$n.eff)),"\n")
    cat("\n DIC = ", mod$BUGSoutput$DIC, "\n")
    # update until model shows no signs of non convergence
    mod <- R2jags::autojags(mod, n.iter = 500000, n.thin = 50, n.update = 3)
  }else{
    message("\n Rhat values < 1.1 indicate lack of evidence for non convergence \n")
    message(" \n Mean effective sample size among parameters = ", round(mean(mod_summary$n.eff)),"\n")
    cat("\n DIC = ", mod$BUGSoutput$DIC, "\n")
    mod <- mod
  }
rm(mod_summary)

# Do a second round of n = 3 if still no convergence
  mod_summary <- mod$BUGSoutput$summary %>% as.data.frame()

  cat("\n Double checking convergence via the Gelman-Rubin statistic (potential scale reduction) comparing within and between chain variance ... \n")
  if (any(mod_summary$Rhat > 1.1)) {
    message("\n Rhat values for one or more of the estimated values are above 1.1 indicating the model has not converged ... updating model with more iterations and treating initial as burn in \n")
    message("\n Mean effective sample size among parameters = ", round(mean(mod_summary$n.eff)),"\n")
    cat("\n DIC = ", mod$BUGSoutput$DIC, "\n")
    # update until model shows no signs of non convergence
    mod <- R2jags::autojags(mod, n.iter = 500000, n.thin = 50, n.update = 3)
  }else{
    message("\n Rhat values < 1.1 indicate lack of evidence for non convergence \n")
    message(" \n Mean effective sample size among parameters = ", round(mean(mod_summary$n.eff)),"\n")
    cat("\n DIC = ", mod$BUGSoutput$DIC, "\n")
    mod <- mod
  }

# Extract mcmc chains (consider tidybayes in the future)
  mod_mcmc <- R2jags:::as.mcmc.rjags(mod)
  modmcmc <- as.data.frame(as.matrix(mod_mcmc))

End of tweak ----

Happy to discuss and hope this is of use to you.

Cheers,

Andrew