runMCMCbtadjust Presentation

Frédéric Gosselin

2023-11-12

Introduction

This file is meant to present the function runMCMC_btadjust() in the runMCMCbtadjust package. The aim of this function is to run a Markov Chain Monte Carlo (MCMC) for a specified Bayesian model while adapting automatically the burn-in and thinning parameters to meet pre-specified targets in terms of MCMC convergence and number of effective values of MCMC outputs - where the term “number of effective values” for the MCMC outputs refers to sample size adjusted for autocorrelation. This is done in only one call to the function that repeatedly calls the MCMC until criteria for convergence and number of effective values are met. This allows to obtain a MCMC output that is out of the transient phase of the MCMC (convergence) and that contains a pre-specified number of nearly independent draws from the posterior distribution (number of effective values).

This function has four main advantages: (i) it saves the analyst’s programming time since he/she does not have to repeatedly diagnose and re-run MCMCs until desired levels of convergence and number of effective values are reached; (ii) it allows a minimal, normalized quality control of MCMC outputs by allowing to meet pre-specified levels in terms of convergence and number of quasi-independent values; (iii) it may save computer’s time when compared to cases where we have to restart the MCMC from the beginning if it has not converged or reached the specified number of effective values (as e.g. with runMCMC function in NIMBLE); and (iv) it can be applied with different MCMC R languages.

Indeed, runMCMC_btadjust() uses other Bayesian packages to fit the MCMC. At present, only the JAGS, NIMBLE and greta languages can be used as these are the main Bayesian languages in R known by the package author and that permit to continue an already fitted MCMC - which is required for numerical efficiency. We will here show how to fit and compare a very simple model under these three languages, using the possibilities allowed by runMCMC_btadjust(). Our model is one of the simplest statistical model we could think of: inspired from @Kery_2010, we model data of weights of 1,000 Pilgrim falcons (Falco peregrinus) simulated from a Gaussian distribution with mean 600 grams and standard error 30 grams:

set.seed(1)
y1000<-rnorm(n=1000,mean=600,sd=30)

This document is made to provide simple examples with the three languages - JAGS, NIMBLE and greta. Yet only the languages that are available/installed on the computer compiling this document will be developed. And if NIMBLE is not available, no example will be developed given that Nimble is the reference example herein.

NIMBLE

We start with fitting the example with NIMBLE (cf. https://r-nimble.org/).


library(runMCMCbtadjust)
library(nimble)
#> Warning: le package 'nimble' a été compilé avec la version R 4.3.1

As NIMBLE distinguishes data that have random distributions from other data, we specify two distinct lists to contain these:


ModelData <-list(mass = y1000)
ModelConsts <- list(nobs = length(y1000))

We then write our Bayesian code within R with the nimbleCode function in the nimble package:

 ModelCode<-nimbleCode(
  {
    # Priors
    population.mean ~ dunif(0,5000)
    population.sd ~ dunif(0,100)
    
    # Normal distribution parameterized by precision = 1/variance in Nimble
    population.variance <- population.sd * population.sd
    precision <- 1 / population.variance
  
    # Likelihood
    for(i in 1:nobs){
      mass[i] ~ dnorm(population.mean, precision)
    }
  })

Our -optional- next step is to specify starting values for model’s parameters. This is done by first writing a function that is repetitively called for each chain. We - also optionally - indicate the names of parameters to be saved and diagnosed in a vector called params:

ModelInits <- function()
{list (population.mean = rnorm(1,600,90), population.sd = runif(1, 1, 30))}
  
Nchains <- 3

set.seed(1)
Inits<-lapply(1:Nchains,function(x){ModelInits()})

#specifying the names of parameters to analyse and save:
params <- c("population.mean", "population.sd") 

We are now ready to launch runMCMC_btadjust(): since we use NIMBLE, we must specify arguments code, data, constants (see below) as well as MCMC_language="Nimble". We first do it on one chain (argument Nchains=1) using in the control list argument neff.method="Coda" to use the Coda method to calculate the number of effective parameters and convtype="Geweke" to use the Geweke method to diagnose convergence, with the pre-specified maximum - over analyzed parameters - convergence of 1.05 (conv.max=1.05) and the minimum - over analyzed parameters - number of effective values of 1,000 (neff.min=1000). Other arguments that are the same for all MCMC languages include params (parameter names to diagnose and save), inits (initial values), niter.min (minimum number of iterations), niter.max (maximum number of iterations), nburnin.min, nburnin.max thin.min, thin.max (minimum and maximum number of iterations for respectively the burn-in and thinning parameters):

out.mcmc.Coda.Geweke<-runMCMC_btadjust(code=ModelCode, constants = ModelConsts, data = ModelData, MCMC_language="Nimble",
    Nchains=1, params=params, inits=Inits[1],
    niter.min=1000, niter.max=300000,
    nburnin.min=100, nburnin.max=200000, 
    thin.min=1, thin.max=1000,
    conv.max=1.05, neff.min=1000,
    control=list(neff.method="Coda", convtype="Geweke"),
    control.MCMC=list(showCompilerOutput=FALSE))
#> ===== Monitors =====
#> thin = 1: population.mean, population.sd
#> ===== Samplers =====
#> RW sampler (2)
#>   - population.mean
#>   - population.sd
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> [1] "Case of niter update: Non convergence"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> [1] "raw multiplier of thin:  25.63"
#> [1] "###################################################################################"
#> [1] "Case of niter update: Convergence and trying to reach end of MCMC at the end of next cycle"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> [1] "Case of niter update: Non convergence"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|

We then run the MCMC with Nchains MCMC chains, the -default- Gelman-Rubin diagnostic of convergence and the -default- rstan method for calculating the number of effective values:

out.mcmc<-runMCMC_btadjust(code=ModelCode, constants = ModelConsts, data = ModelData, MCMC_language="Nimble",
    Nchains=Nchains, params=params, inits=Inits,
    niter.min=1000, niter.max=300000,
    nburnin.min=100, nburnin.max=200000, 
    thin.min=1, thin.max=1000,
    conv.max=1.05, neff.min=1000,
    control.MCMC=list(showCompilerOutput=FALSE))
#> ===== Monitors =====
#> thin = 1: population.mean, population.sd
#> ===== Samplers =====
#> RW sampler (2)
#>   - population.mean
#>   - population.sd
#> ===== Monitors =====
#> thin = 1: population.mean, population.sd
#> ===== Samplers =====
#> RW sampler (2)
#>   - population.mean
#>   - population.sd
#> ===== Monitors =====
#> thin = 1: population.mean, population.sd
#> ===== Samplers =====
#> RW sampler (2)
#>   - population.mean
#>   - population.sd
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> [1] "Case of niter update: Non convergence"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> [1] "raw multiplier of thin:  4.791"
#> [1] "###################################################################################"
#> [1] "Case of niter update: Convergence and trying to reach end of MCMC at the end of next cycle"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> [1] "raw multiplier of thin:  1.338"
#> [1] "###################################################################################"
#> [1] "Case of niter update: Convergence and trying to reach end of MCMC at the end of next cycle"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|

We compare the characteristics of the two MCMCs, both in terms of burn-in, thinning parameter, number of iterations and in terms of time (both total time and CPU time).

Table: Comparison of the efficiency of first two NIMBLE models:

Nimble.Coda.Geweke Nimble.default
converged 1.000000 1.000000
neffs.reached 1.000000 1.000000
burnin 102.000000 575.000000
thin 26.000000 5.000000
niter.tot 41701.000000 3065.000000
duration 31.190509 60.680501
duration.MCMC.preparation 24.184504 54.622998
duration.MCMC.transient 0.009065 0.158048
duration.MCMC.asymptotic 3.694719 0.684417
duration.btadjust 3.302221 5.215039
CPUduration 10.440000 15.410000
CPUduration.MCMC.preparation 6.190000 13.360000
CPUduration.MCMC.transient 0.009105 0.165090
CPUduration.MCMC.asymptotic 3.710895 0.714910
CPUduration.btadjust 0.530000 1.170000

We acknowledge that the Coda.Geweke algorithm takes much less time (rows named duration and CPUduration in the previous table) than the classical setting to prepare data (rows named duration.MCMC.preparation and CPUduration.MCMC.preparation)- as NIMBLE takes quite a lot of time to prepare each MCMC chain - and we have 3 chains to prepare in the default setting compared to 1 with Geweke.

We also notice that the Coda.Geweke algorithm uses more time (duration.MCMC.transient and CPUduration.MCMC.transient) and iterations (burnin) to converge while the default setting takes more time for the asymptotic phase (duration.MCMC.asymptotic and CPUduration.MCMC.asymptotic), linked to a greater thinning parameter (thin). This transient part should be linked to the different behaviors between Geweke and Gelman-Rubin convergence diagnostics, while the differences in thinning parameters might be linked to the different methods in calculating the number of effective parameters. We therefore run a third MCMC on one chain with Geweke diagnostic but the default, rstan method for number of effective values.


out.mcmc.Geweke<-runMCMC_btadjust(code=ModelCode, constants = ModelConsts, data = ModelData, MCMC_language="Nimble",
    Nchains=1, params=params, inits=Inits[1],
    niter.min=1000, niter.max=300000,
    nburnin.min=100, nburnin.max=200000, 
    thin.min=1, thin.max=1000,
    conv.max=1.05, neff.min=1000,
    control=list(convtype="Geweke"),
    control.MCMC=list(showCompilerOutput=FALSE))
#> ===== Monitors =====
#> thin = 1: population.mean, population.sd
#> ===== Samplers =====
#> RW sampler (2)
#>   - population.mean
#>   - population.sd
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> [1] "Case of niter update: Non convergence"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> [1] "raw multiplier of thin:  4.533"
#> [1] "###################################################################################"
#> [1] "Case of niter update: Convergence and trying to reach end of MCMC at the end of next cycle"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> [1] "raw multiplier of thin:  1.251"
#> [1] "###################################################################################"
#> [1] "Case of niter update: Convergence and trying to reach end of MCMC at the end of next cycle"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|

We compare the characteristics of the three NIMBLE MCMCs,

Table: Comparison of the efficiency of the three NIMBLE models:

Nimble.Coda.Geweke Nimble.Geweke Nimble.default
converged 1.000000 1.000000 1.000000
neffs.reached 1.000000 1.000000 1.000000
burnin 102.000000 1525.000000 575.000000
thin 26.000000 5.000000 5.000000
niter.tot 41701.000000 8720.000000 3065.000000
duration 31.190509 24.313382 60.680501
duration.MCMC.preparation 24.184504 18.539179 54.622998
duration.MCMC.transient 0.009065 0.143196 0.158048
duration.MCMC.asymptotic 3.694719 0.675605 0.684417
duration.btadjust 3.302221 4.955402 5.215039
CPUduration 10.440000 5.670000 15.410000
CPUduration.MCMC.preparation 6.190000 4.380000 13.360000
CPUduration.MCMC.transient 0.009105 0.143406 0.165090
CPUduration.MCMC.asymptotic 3.710895 0.676594 0.714910
CPUduration.btadjust 0.530000 0.470000 1.170000

Results do not corroborate our above expectations: indeed, first the Geweke model does not converge within the same number of iterations as the Coda.Geweke one (row burnin), which is strange since they use the same method for diagnosing convergence and the same seed. Second, the thinning parameter was not increased when changing from Coda.Geweke to Geweke as expected above, but actually decreased (row thin).

We now turn to the comparison of the statistical parameter outputs. We use two sample Kolmogorov-Smirnov tests to compare each parameter by pairs of MCMC methods:

Table: P-values of paired Kolmogorov-Smirnov tests of output parameters (columns) between the first three NIMBLE models (rows):

mean sd
default vs. Geweke 0.01245 0.19455
Coda.Geweke vs. Geweke 0.11718 0.52490
Default vs. Coda.Geweke 0.29689 0.68044

The p-values associated to the KS tests are not very small - only one out of six is near 0.05. This indicates that the MCMC outputs can be considered as being drawn from the same distributions.

These parameters are summarized in the next tables.

Table: Summary of the statistical parameters of the Nimble Coda.Geweke model:

Mean SD Naive SE Time-series SE
population.mean 599.602 1.111 0.028 0.029
population.sd 31.090 0.726 0.018 0.018

Table: Summary of the statistical parameters of the Nimble Geweke model:

Mean SD Naive SE Time-series SE
population.mean 599.688 0.979 0.026 0.028
population.sd 31.081 0.697 0.018 0.021

Table: Summary of the statistical parameters of the Nimble default model:

Mean SD Naive SE Time-series SE
population.mean 599.589 0.939 0.024 0.025
population.sd 31.112 0.691 0.018 0.020

We notice that parameter values are very close, that naive standard errors (SEs) are very close to Time-series SEs - which is linked to the automatic tuning of the thinning parameter which produces output samples which are nearly independent - and that differences between mean estimators are within several units of Time series-SEs - which we interpret is mostly due to the control of convergence.

JAGS

We now turn to analyzing the same data with the same statistical model using JAGS with runMCMC_btadjust(). We rely on the data simulated above. In JAGS, we now put all the data in the same list:


ModelData <-list(mass = y1000, nobs = length(y1000))

We then propose the use of JAGS with a specification of the model from within R - which we find more convenient. We therefore write the model within R as a character chain:


modeltotransfer<-"model {

		# Priors
			population.mean ~ dunif(0,5000)
			population.sd ~ dunif(0,100)

			# Normal distribution parameterized by precision = 1/variance in Jags
    	population.variance <- population.sd * population.sd
      precision <- 1 / population.variance

			# Likelihood
			for(i in 1:nobs){
			  mass[i] ~ dnorm(population.mean, precision)
			}
		}"

The other objects useful or required for running runMCMC_btadjust with JAGS are similar to those required with NIMBLE (Inits, Nchains, params) and are not repeated here.

We then launch runMCMC_btadjust() with MCMC_language="Jags", specifying arguments code and data which are required in this case:


set.seed(1)
out.mcmc.Jags<-runMCMC_btadjust(code=modeltotransfer,  data = ModelData, MCMC_language="Jags", 
    Nchains=Nchains, params=params, inits=Inits,
    niter.min=1000,niter.max=300000,
    nburnin.min=100,nburnin.max=200000,
    thin.min=1,thin.max=1000,
		conv.max=1.05,neff.min=1000)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 1000
#>    Unobserved stochastic nodes: 2
#>    Total graph size: 1009
#> 
#> Initializing model
#> 
#> [1] "###################################################################################"

Note that if we had written the JAGS code in a text file named "ModelJags.txt", we would just have replaced in the command above code=modeltotransfer by code="ModelJags.txt".

Table: Summary of the statistical parameters of the Jags model:

Mean SD Naive SE Time-series SE
population.mean 599.645 1.015 0.020 0.023
population.sd 31.078 0.700 0.013 0.017

Results seem in line with those of NIMBLE. We check this using a paired Kolmogorov-Smirnov tests with NIMBLE models:

Table: P-values of paired Kolmogorov-Smirnov tests of output parameters (columns) of the Jags model with the three NIMBLE models (rows):

mean sd
Nimble.Geweke vs. Jags 0.1685 0.4334
Nimble.Coda.Geweke vs. Jags 0.3362 0.8447
Nimble.Default vs. Jags 0.1240 0.5962

Our results do confirm that the JAGS result cannot be considered as stemming from a different probability distribution than NIMBLE results.

We finally compare the efficiency of the JAGS and default NIMBLE MCMCs:

Table: Comparison of the efficiency of the default NIMBLE model and the Jags model:

Nimble.default Jags
converged 1.0000 1.0000
neffs.reached 1.0000 1.0000
burnin 575.0000 101.0000
thin 5.0000 1.0000
niter.tot 3065.0000 1000.0000
duration 60.6805 6.0822
duration.MCMC.preparation 54.6230 2.8543
duration.MCMC.transient 0.1580 0.2110
duration.MCMC.asymptotic 0.6844 1.8783
duration.btadjust 5.2150 1.1386
CPUduration 15.4100 4.3000
CPUduration.MCMC.preparation 13.3600 2.1600
CPUduration.MCMC.transient 0.1651 0.2121
CPUduration.MCMC.asymptotic 0.7149 1.8879
CPUduration.btadjust 1.1700 0.0400

The conclusion is that JAGS is much faster than NIMBLE on this example (row named duration in the previous table), due to much less time devoted to MCMC preparation - as well as to burn-in/thinning adjustment (rows named duration.MCMC.preparation and duration.btadjust in the previous table). Actually there is no adjustment with JAGS (niter.tot is equal to the initial number of iterations). Yet, NIMBLE is quicker regarding MCMC updating by iteration since it took NIMBLE less time than JAGS for the transient phase (respectively less than twice time for the asymptotic phase) although using more than seven (resp. seventeen) times more iterations than JAGS.

At first sight, we would also conclude that MCMC efficiency per effective value is also better with NIMBLE since both languages had the same target for the minimum number of effective value - 1,000 - and the total MCMC time was lower with NIMBLE. Yet, a more rigorous comparison finds the reverse result:

Table: Comparison of the number of effective values between the default NIMBLE model and the JAGS model:

Nimble.default Jags
Min. Number Eff. values 1074.0000000 1728.0000000
MCMC CPU time per Effective Value 0.0008194 0.0012153

Indeed, JAGS with just the first iterations produced a high number of effective values - actually much bigger than the targeted neff.min - which renders the MCMC time per effective value now lower with JAGS than NIMBLE with this model (cf. table above).

Greta

We finally run the greta version of our model with runMCMC_btadjust(). greta is rather different from JAGS and NIMBLE in that the model defines objects in R and thus does not require a model code to be passed to runMCMC_btadjust(), nor Data or Constants. We rely on the very data simulated above. The coding with greta is as follows:

#in my setting I need to load not only greta but R6 & tensorflow packages
library(greta)
library (R6)
library(tensorflow)

#first requirement of greta: declaring the data that will be analyzed with the function as_data
Y<-as_data(y1000)

#we then proceed by writing the model directly in R, starting with the priors of the parameters using greta functions for probability distributions - here uniform()
population.mean<-uniform(0,5000)
population.sd<-uniform(0,100)
    
#we then define the distribution of the data - here with the normal distribution - by default parametrized with a standard deviation in greta:
try({distribution(Y)<-normal(population.mean,population.sd) })

#we finally declare the greta model, which will be the object passed to runMCMC_btadjust 
m<-model(population.mean, population.sd)

### we finally have to prepare initial values with a specific greta function - initials:
ModelInits <- function()
    {initials(population.mean = rnorm(1,600,90), population.sd = runif(1, 1, 30))}

set.seed(1)
  Inits<-lapply(1:Nchains,function(x){ModelInits()})

We are now ready to fit the model with runMCMC_btadjust(), specifying MCMC_language="Greta" and giving the argument model instead of code and data:

out.mcmc.greta<-runMCMC_btadjust(model=m, MCMC_language="Greta",
    Nchains=Nchains,params=params,inits=Inits,
		niter.min=1000,niter.max=300000,
    nburnin.min=100,nburnin.max=200000,
		thin.min=1,thin.max=1000,
		conv.max=1.05, neff.min=1000)
#> [1] "###################################################################################"
#> 
#> [1] "raw multiplier of thin:  3.807"
#> [1] "###################################################################################"
#> [1] "Case of niter update: Convergence and trying to reach end of MCMC at the end of next cycle"
#> [1] "###################################################################################"
#> [1] "###################################################################################"
#> 
#> [1] "raw multiplier of thin:  1.414"
#> [1] "###################################################################################"
#> [1] "Case of niter update: Convergence and trying to reach end of MCMC at the end of next cycle"
#> [1] "###################################################################################"
#> [1] "###################################################################################"

Table: Summary of the statistical parameters of the greta model:

Mean SD Naive SE Time-series SE
population.mean 599.666 1.002 0.026 0.032
population.sd 31.099 0.658 0.017 0.016

We first check that estimations are similar to those with NIMBLE and JAGS with paired Kolmogorov-Smirnov tests:

Table: P-values of paired Kolmogorov-Smirnov tests of output parameters of the greta model with the default NIMBLE model and the JAGS model:

mean sd
Nimble vs. greta 0.031532 0.024661
Jags vs. greta 0.384392 0.007422

We then report the efficiency of the MCMCs.

Table: Comparison of the efficiency of the default NIMBLE, the JAGS and the greta models:

Nimble.default Jags Greta
converged 1.0000 1.0000 1.0000
neffs.reached 1.0000 1.0000 1.0000
burnin 575.0000 101.0000 1.0000
thin 5.0000 1.0000 4.0000
niter.tot 3065.0000 1000.0000 2056.0000
duration 60.6805 6.0822 35.7248
duration.MCMC.preparation 54.6230 2.8543 0.3699
duration.MCMC.transient 0.1580 0.2110 10.4235
duration.MCMC.asymptotic 0.6844 1.8783 21.3989
duration.btadjust 5.2150 1.1386 3.5326
CPUduration 15.4100 4.3000 91.6900
CPUduration.MCMC.preparation 13.3600 2.1600 0.0200
CPUduration.MCMC.transient 0.1651 0.2121 29.9710
CPUduration.MCMC.asymptotic 0.7149 1.8879 61.5290
CPUduration.btadjust 1.1700 0.0400 0.1700

MCMC time (rows duration.MCMC.transient & duration.MCMC.asymptotic) was far greater with greta than with JAGS and NIMBLE, for a minimum number of effective values with greta of 1024. Total duration is rather close with greta compared with NIMBLE, due to the great time required by NIMBLE for MCMC preparation - while this preparation is done outside runMCMC_btadjust() with greta. Yet, when we compare CPU total durations (CPUduration), greta gets worse than NIMBLE, simply because greta parallelized its process and therefore required more CPU time per time unit.

We tried to give a second chance to greta, based on the following post: https://forum.greta-stats.org/t/size-and-number-of-leapfrog-steps-in-hmc/332. The idea was to let greta have more information to adapt its hmc parameters during the warm-up phase by just having more chains to run - hereafter, 15.

Nchains<-15
ModelInits <- function()
    {initials(population.mean = rnorm(1,600,90), population.sd = runif(1, 1, 30))}

set.seed(1)
Inits<-lapply(1:Nchains,function(x){ModelInits()})
  
  out.mcmc.greta.morechains<-runMCMC_btadjust(model=m, MCMC_language="Greta",
    Nchains=Nchains,params=params,inits=Inits,
		niter.min=1000,niter.max=300000,
    nburnin.min=100,nburnin.max=200000,
		thin.min=1,thin.max=1000,
		conv.max=1.05, neff.min=1000)
#> [1] "###################################################################################"

Table: Summary of the statistical parameters of the greta model with 15 chains:

Mean SD Naive SE Time-series SE
population.mean 599.6268 0.9913 0.0081 0.0160
population.sd 31.0860 0.6943 0.0057 0.0118

This run was indeed much faster. Parameter estimates were still not significantly different from those with NIMBLE and JAGS based on paired Kolmogorov-Smirnov tests:

Table: P-values of paired Kolmogorov-Smirnov tests of output parameters of the greta model with 15 chains with the default NIMBLE model and the JAGS model:

mean sd
Nimble vs. greta.morechains 0.03657 0.47788
Jags vs. greta.morechains 0.65206 0.70921

We now report the efficiency of the MCMCs:

Table: Comparison of the efficiency of the default NIMBLE, the JAGS and the greta.morechains models:

Nimble.default Jags Greta.morechains
converged 1.0000 1.0000 1.0000
neffs.reached 1.0000 1.0000 1.0000
burnin 575.0000 101.0000 1.0000
thin 5.0000 1.0000 1.0000
niter.tot 3065.0000 1000.0000 1000.0000
duration 60.6805 6.0822 22.2909
duration.MCMC.preparation 54.6230 2.8543 0.3810
duration.MCMC.transient 0.1580 0.2110 10.2547
duration.MCMC.asymptotic 0.6844 1.8783 10.2342
duration.btadjust 5.2150 1.1386 1.4209
CPUduration 15.4100 4.3000 71.3300
CPUduration.MCMC.preparation 13.3600 2.1600 0.0000
CPUduration.MCMC.transient 0.1651 0.2121 35.6056
CPUduration.MCMC.asymptotic 0.7149 1.8879 35.5344
CPUduration.btadjust 1.1700 0.0400 0.1900

We still observed more CPU duration with greta, although the difference had decreased and the associated number of effective values for greta was now 3291, which rendered MCMC CPU efficiency with greta closer to NIMBLE. Yet, overall, JAGS performed better than greta and NIMBLE on this example.

Conclusion

We hope we have convinced the R user of Bayesian models that runMCMC_btadjust() can help having a more efficient, quality oriented use of these types of models while saving analyst’s and potentially computer time. Indeed, to recap, the aim of this function is to run a Markov Chain Monte Carlo (MCMC) for a specified Bayesian model while adapting automatically the burn-in and thinning parameters to meet pre-specified targets in terms of MCMC convergence and number of effective values of MCMC outputs. This is done in only one call to the function that repeatedly calls the MCMC until criteria for convergence and number of effective values are met. The function has four main advantages:

(i) it saves the analyst’s programming time since he/she does not have to repeatedly diagnose and re-run MCMCs until desired levels of convergence and number of effective values are reached;

(ii) it allows a minimal, normalized quality control of MCMC outputs by allowing to meet pre-specified levels in terms of convergence and number of quasi-independent values;

(iii) it may save computer’s time when compared to cases where we have to restart the MCMC from the beginning if it has not converged or reached the specified number of effective values;

(iv) it can be applied with different MCMC R languages - at present greta, NIMBLE and JAGS. This comes with two positive consequences in practice: first, allowing the user a more rigorous comparison between the three Bayesian fitting languages in terms of comparability of inference and of MCMC efficiency - especially in terms of CPU time per effective value; second, making it easier to develop the same Bayesian model with these different languages, which is to our experience welcome in practical cases, since these different languages have advantages over the other ones that vary from one context to the other.

References