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

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

predC = jabba$est.catch
 plot(years, (jabba$catch), type = "n", ylim = c(0, max(predC, 
            na.rm = T)), lty = 1, lwd = 1.3, xlab = "Year", ylab = paste0("Catch ", 
            jabba$settings$catch.metric), main = "")

ylim = c(0, max(predC, na.rm = T)

ylim = c(0, max(predC[,-1], na.rm = T)

> checking for missing documentation entries ... WARNING
  Undocumented code objects:
    'iccat'
  Undocumented data sets:
    'iccat'
  All user-level objects in a package should have documentation entries.
  See chapter 'Writing R documentation files' in the 'Writing R
  Extensions' manual.

❯ checking for missing documentation entries ... WARNING
  Undocumented code objects:
    'catch' 'catch.n' 'iccat' 'index' 'om'
  Undocumented data sets:
    'iccat' 'catch' 'catch.n' 'index' 'om'
  All user-level objects in a package should have documentation entries.
  See chapter 'Writing R documentation files' in the 'Writing R
  Extensions' manual.

 Running examples in 'JABBA-Ex.R' failed
  The error most likely occurred in:
  
  > base::assign(".ptime", proc.time(), pos = "CheckExEnv")
  > ### Name: jbplot_ensemble
  > ### Title: jbplot_ensemble()
  > ### Aliases: jbplot_ensemble
  > 
  > ### ** Examples
  > 
  > data(iccat)
  > bet <- iccat$bet 
  > 
  > # Fit Fox and Schaefer
  > jb1 <- build_jabba(catch=bet$catch, 
  +                    cpue=bet$cpue,
  +                    se=bet$se,
  +                    scenario = "Fox",
  +                    model.type="Fox")
  
  ><> Prepare JABBA input data <><
  
  
  
  ><> Assume Catch with error CV = 0.1 <><
  
  
  
  ><> Model type:Fox <><
  
  
  ><> Shape m =1.001
  
  
  ><> K prior mean =1079485.2368and CV =1(log.sd = 0.832554611157698)
  
  
  ><> r prior mean =0.2and CV =0.532940350027788(log.sd = 0.5)
  
  
  ><> Psi (B1/K) prior mean =0.9and CV =0.25withlnormdestribution
  
  
  
  
  ><> ALWAYS ENSURE to adjust default settings to your specific stock <><
  
  
  > jb2 <- build_jabba(catch=bet$catch,
  +                    cpue=bet$cpue,
  +                    se=bet$se,
  +                    scenario = "Schaefer",
  +                    model.type="Schaefer")
  
  ><> Prepare JABBA input data <><
  
  
  
  ><> Assume Catch with error CV = 0.1 <><
  
  
  
  ><> Model type:Schaefer <><
  
  
  ><> Shape m =2
  
  
  ><> K prior mean =1079485.2368and CV =1(log.sd = 0.832554611157698)
  
  
  ><> r prior mean =0.2and CV =0.532940350027788(log.sd = 0.5)
  
  
  ><> Psi (B1/K) prior mean =0.9and CV =0.25withlnormdestribution
  
  
  ><> ALWAYS ENSURE to adjust default settings to your specific stock <><
  
  
  > fit1 <- fit_jabba(jb1,
  +                   quickmcmc=TRUE,
  +                   verbose=TRUE)
  Compiling model graph
     Resolving undeclared variables
     Allocating nodes
  Graph information:
     Observed stochastic nodes: 199
     Unobserved stochastic nodes: 358
     Total graph size: 4305
  
  Initializing model
  
  
  ><> Produce results output of Fox model for jabba Fox <><
  
  
  
  ><> Scenario Fox_Fox completed in 0 min and 33 sec <><
  
  > fit2 <- fit_jabba(jb2,
  +                   quickmcmc=TRUE,
  +                   verbose=TRUE)
  Compiling model graph
     Resolving undeclared variables
     Allocating nodes
  Graph information:
     Observed stochastic nodes: 199
     Unobserved stochastic nodes: 358
     Total graph size: 4305
  
  Initializing model
  
  
  ><> Produce results output of Schaefer model for jabba Schaefer <><
  
  
  
  ><> Scenario Schaefer_Schaefer completed in 0 min and 33 sec <><
  
  > # Compare
  > jbplot_ensemble(list(fit1,fit2))
  Error in kb[kb$type == "prj", ] : incorrect number of dimensions
  Calls: jbplot_ensemble -> pmin
  Execution halted

✔  checking examples (15m 33.7s)
   Examples with CPU (user + system) or elapsed time > 5s
                     user system elapsed
   hindcast_jabba  264.82   6.22  271.19
   jbhcxval        263.20   6.10  269.39
   jbretro          88.49   2.36   90.94
   jbplot_retro     88.10   2.41   90.55
   jbplot_ensemble  59.08   2.69   61.79
   jbplot_Crisk     41.68   2.03   43.75
   mp_jabba         34.64   0.67   35.33
   fit_jabba        22.39   0.86   23.27
   jbrunstest       22.54   0.58   23.14
   jbplot_runstest  22.43   0.55   22.97