Introduction

compareMCMCs is a package for running, managing, and comparing results from different MCMC packages. You can configure a set of MCMCs to run, have them automatically timed, and then generate html output with comparisons of efficiency and posterior distributions. Throughout this document, we use “MCMC” as a short-hand for an MCMC kernel or the samples it generates.

This system started life as part of the nimble package. Although it remains somewhat nimble-centric, it is entirely possible to use it without nimble involved. The system is highly configurable, allowing use of other MCMC systems and generation of other output metrics and graphical summaries fairly easily.

Use of other MCMCs is supported by a plugin system. Plugins are provided for JAGS and Stan. Draft plugins for WinBUGS and OpenBUGS are provided but need testing. Since nimble, JAGS, WinBUGS and OpenBUGS use different dialects of the same model language, it is sometimes possible to compare them using the same model code. It is possible to write new plugins to call other MCMC packages fairly easily.

compareMCMCs provides:

the compareMCMCs function to run one or more MCMCs and manage the results;
the MCMCresult class to manage results by storing samples, timing information, and metrics or summaries of performance;
a plugin system to include new MCMC engines;
a plugin system for new metrics for comparison among MCMCs;
a system for applying parameter conversions, in case difference MCMCs use different parameterizations and/or parameter names;
a system for generating html pages with figures from comparison metrics, including a plugin system to provide new page components;

Beware that you might need to call the compareMCMCs function as compareMCMCs::compareMCMCs. That’s because there is a deprecated compareMCMCs function in nimble. Eventually that will be removed from nimble. compareMCMCs is not backward-compatible with nimble’s original features for comparing MCMC methods.

An example using `nimble`, `jags`, and `stan`.

Let’s say you want to compare the performance of adaptive random-walk Metropolis-Hastings sampling with nimble (its default sampler in this case), slice sampling with nimble, JAGS (which will also use slice sampling in this case), and Stan. We’ll use a toy model with only one observation following a gamma distribution, known rate parameter, and unknown shape parameter with a uniform prior distribution. The following code accomplishes this:

# This model code will be used for both nimble and JAGS
modelCode <- nimbleCode({
  a ~ dunif(0, 100)
  y ~ dgamma(a, 2)
})
modelInfo <- list(
  code = modelCode,
  constants = list(y = 2), # data can be included in constants or data list
  inits = list(a = 1)
)
# Here is a custom MCMC configuration function for nimble
configure_nimble_slice <- function(model) {
  configureMCMC(model, onlySlice = TRUE)
}
# Here is information to set the model up for Stan
stan_code <- c("data {real y;}",
               "parameters {real a;}",
               "model {target += uniform_lpdf(a | 0, 100);",
               "       target += gamma_lpdf(y | a, 2);}")
# By default, work_dir is R's temporary session directory, tempdir().
# If you want to change it, modify the first code chunk of the R markdown source.
writeLines(stan_code, file.path(work_dir, "toy.stan"))
# An alternative way to provide the Stan model is to make a
# list entered below as stan_model_args (similar to how sampling_args
# appears in the externalMCMCinfo list).  The elements of stan_model_args
# (e.g. 'file' or 'model_code') will be passed to the 'stan_model` function.
stan_sampling_args <- list(data = list(y=2),
                           init = list(list(a=1)), # Only one chain will be run
                           warmup = 200, # Stan was not happy with warmup = 100
                           iter = 2000)
# Here is the call to compareMCMCs
res <- compareMCMCs(modelInfo,
                    MCMCs = c(if(this_system_has_rjags) 'jags' else NULL,
                              'nimble',       # nimble with default samplers
                              'nimble_slice', # nimble with slice samplers
                              if(this_system_has_rstan) 'stan' else NULL),
                    nimbleMCMCdefs = 
                      list(nimble_slice = 'configure_nimble_slice'),
                    MCMCcontrol = list(inits = list(a = 1),
                                       niter = 2000,
                                       burnin = 100),
                    externalMCMCinfo = list(stan = list(
                      file = file.path(work_dir, "toy.stan"),
                      sampling_args = stan_sampling_args))
                    # Stan was not happy with warmup = 100, which would match the others.
)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 1
#>    Unobserved stochastic nodes: 1
#>    Total graph size: 5
#> 
#> Initializing model
#> 
#> 
#> SAMPLING FOR MODEL 'toy' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 1.1e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  201 / 2000 [ 10%]  (Sampling)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Sampling)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Sampling)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Sampling)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.005292 seconds (Warm-up)
#> Chain 1:                0.018267 seconds (Sampling)
#> Chain 1:                0.023559 seconds (Total)
#> Chain 1:
#> Warning: There were 29 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> ===== Monitors =====
#> thin = 1: a
#> ===== Samplers =====
#> RW sampler (1)
#>   - a
#> ===== Monitors =====
#> thin = 1: a
#> ===== Samplers =====
#> slice sampler (1)
#>   - a
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
file.remove(work_dir, "toy.stan")
#> Warning in file.remove(work_dir, "toy.stan"): cannot remove file '/Users/perry/
#> Projects/GitHub/nimble-dev/compareMCMCs/compareMCMCs/inst/vignette_source',
#> reason 'Directory not empty'
#> [1] FALSE  TRUE
make_MCMC_comparison_pages(res, dir = work_dir,  modelName = 'example1')

Results are here. In this case, all samplers are so efficient that the comparison between them is not very stable and not meaningful for such a toy model. On the machine this is being written on, nimble_slice is the most efficient MCMC.

Data vs. constants in nimble vs. JAGS and WinBUGS/OpenBUGS

When using the same model code for nimble and any of JAGS, WinBUGS or OpenBUGS, it is important to note that nimble makes a distinction between “data” and “constants”, which the other packages group together as “data”. In nimble, constants are values necessary to define a model, such as “N” in “for(i in 1:N)”, while data are observed scientific data. Since nimble models are objects that can be manipulated, it is possible to change data values, although that won’t happen just from MCMC sampling.

When both data and constants are provided as “constants” to nimble, it attempts to determine which values should be treated as data. Correspondingly, the constants will be passed to the other packages as data. Therefore, it is advisable to provide both data and constants as “constants” when using both nimble and JAGS, WinBUGS or OpenBUGS. See the NIMBLE User Manual (https://r-nimble.org) for more information about data and constants in nimble.

MCMC efficiency and pace

MCMC efficiency for a parameter is defined as the effective sample size divided by computation time (in seconds). It is the number of effectively independent samples generated per second.

MCMC pace is defined as the inverse of MCMC efficiency. It is the time (in seconds) required to generate one effectively independent sample.

What is timed

For both of these metrics, one may take interest in different components of computation time. compareMCMCs records time in three components: setup time, burnin time (“warmup” in Stan), and post-burnin time. The samples that are saved for analysis are the post-burnin samples. A fourth time, sampling time, is defined as the sum of burnin and post-burnin times. Setup time includes any steps prior to sampling. It is possible to control which time is used for calculating MCMC efficiency and pace. The default is sampling time. This focuses on comparisons once algorithms are up and running.

New MCMC configurations in nimble

NIMBLE provides the capability to customize an MCMC configuration by choosing samplers, writing new samplers, and arranging sampler order. To have compareMCMCs use a custom configuration, the steps are:

Choose a name for the configuration, such as “nimble_slice” above. Include that name in the MCMCs argument. Some names are pre-defined, such as “nimble” and “jags”. (Actually, “nimble_slice” is provided as a pre-defined configuration. We have just used it as an example.)
Provide an element in the nimbleMCMCdefs list with the same name as entered in MCMCs (e.g., “nimble_slice”). That element can take one of three formats:
1. A function, which should take a nimble model as its first and only argument and return an MCMC configuration (e.g., what is returned by configureMCMC, possibly modified using its addSampler and removeSamplers methods).
2. The character name of a function as described in (i).
3. An R code object such as that created by quote. This can assume there is an object called model that is a nimble model. It should return an MCMC configuration. An example that would work above is nimbleMCMCdefs = list(nimble_slice = quote({configureMCMC(model, onlySlice = TRUE)})). The curly braces can contain as many lines of code as needed, separated by semi-colons.

Pre-defined configurations

The package comes with the following configurations, which are similar to options in nimble’s configureMCMC function.

“nimble_slice” : Use only slice sampling where possible.
“nimble_noConj”: Do not use conjugate (Gibbs) samplers even where they are valid. Instead use adaptive random-walk Metropolis-Hastings samplers.
“nimble_RW”: Use only adaptive random-walk Metropolis-Hastings samplers.

Combining results from different calls to `compareMCMCs`.

Results are stored as a list of MCMCresult objects, so it is possible to simply concatenate lists to combine results. This makes it possible to run and save results separately and combine them later for comparison purposes.

For example, the above example (omitting Stan for simplicity) could be created as follows:

res_jags <- compareMCMCs(modelInfo,
                MCMCs = c('jags'),
                MCMCcontrol = list(inits = list(a = 1),
                                   niter = 2000,
                                   burnin = 100))
#> building nimble model...
#> Defining model
#>   [Note] Using 'y' (given within 'constants') as data.
#> Building model
#> Setting data and initial values
#> Running calculate on model
#>   [Note] Any error reports that follow may simply reflect missing values in model variables.
#> Checking model sizes and dimensions
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 1
#>    Unobserved stochastic nodes: 1
#>    Total graph size: 5
#> 
#> Initializing model
# To illustrate that each MCMC can have its own settings, we run nimble and jags
# different numbers of iterations.
res_nimble <- compareMCMCs(modelInfo,
                MCMCs = c('nimble', 'nimble_slice'),
                nimbleMCMCdefs = list(nimble_slice = 'configure_nimble_slice'),
                MCMCcontrol = list(inits = list(a = 1),
                                   niter = 4000,
                                   burnin = 200))
#> building nimble model...
#> Defining model
#>   [Note] Using 'y' (given within 'constants') as data.
#> Building model
#> Setting data and initial values
#> Running calculate on model
#>   [Note] Any error reports that follow may simply reflect missing values in model variables.
#> Checking model sizes and dimensions
#> ===== Monitors =====
#> thin = 1: a
#> ===== Samplers =====
#> RW sampler (1)
#>   - a
#> ===== Monitors =====
#> thin = 1: a
#> ===== Samplers =====
#> slice sampler (1)
#>   - a
#> Compiling
#>   [Note] This may take a minute.
#>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
#> Compiling
#>   [Note] This may take a minute.
#>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
res <- c(res_jags, res_nimble)
make_MCMC_comparison_pages(res, dir = work_dir, modelName = 'example2')

Results are here.

Writing new comparison metrics.

The result of compareMCMCs is a list of MCMCresult objects. In addition to the samples (MCMC output), each MCMCresult has elements for times and metrics. The metrics is a list of metrics organized by their relationship to model parameters. Metrics that have a single value for each MCMC sample are stored in the byMCMC element of metrics. Metrics that have a value for each parameter in each MCMC sample are stored in the byParameter element of metrics. Finally, there is an other element of metrics that is a list available for storing any arbitrary kinds of metrics (or any kind of information really) for each MCMC.

For example, here are the metrics of the result for nimble with default samplers above:

res$nimble$metrics
#> $byMCMC
#>     MCMC min_efficiency mean_efficiency
#> 1 nimble       427454.7        427454.7
#> 
#> $byParameter
#>     MCMC Parameter     mean   median      sd CI95_low CI95_upp      ESS
#> 1 nimble         a 4.933169 4.748909 2.01564 1.524679 9.271893 854.9093
#>   efficiency
#> 1   427454.7
#> 
#> $other
#> list()

A new comparison metric can be provided as a function. It takes as input a MCMCresult object and returns a list with elements named byMCMC, byParameter, and/or other. A result in byParameter should be a list whose names are metric names. Each list element should be a named vector with names corresponding to parameters.

For example, the metric function for median looks like this:

MCMCmetric_median <- function(result, ...) {
  res <- apply(result$samples, 2, median)
  list(byParameter = list(median = res))
}

When this is called by other functions, the argument result will be the object returned from combineMetrics, such as combineMetrics(res) (see below for an example). When called from addMetrics or compareMCMCs, the corresponding element of their options or metricOptions arguments, respectively, will be passed as the options argument to this function.

There are two ways metrics can be used:

By calling addMetrics with a list of metrics, which can be functions or registered metric names.
By registering each metric with a name and including those names in the metrics argument to compareMCMCs.

Let’s look at an example of each. Say we want to make a metric that calculates the 25th and 75th percentiles (i.e. the quartiles) for each parameter and also records the maximum difference between these percentiles across all parameters. We would do so as follows:

MCMCmetric_quartiles <- function(result, options) {
  p25 <- apply(result$samples, 2, quantile, probs = 0.25)
  p75 <- apply(result$samples, 2, quantile, probs = 0.75)
  # q25 and q75 are named vectors with names matching model parameters
  # i.e. column names of result$samples
  maxDiff <- max(p75-p25)
  list(byParameter = list(p25 = p25,
                          p75 = p75),
       byMCMC = list(maxQuartileDiff = maxDiff))
}

Note that familiarity with the contents of a MCMCresults object is needed to provide a new metric. Often it is sufficient to know that its samples field contains the matrix of samples.

We can now add the metric to the set of stored metric results from our previous results:

addMetrics(res, list(MCMCmetric_quartiles))
res$nimble$metrics
#> $byMCMC
#>     MCMC min_efficiency mean_efficiency maxQuartileDiff
#> 1 nimble       427454.7        427454.7        2.796365
#> 
#> $byParameter
#>     MCMC Parameter     mean   median      sd CI95_low CI95_upp      ESS
#> 1 nimble         a 4.933169 4.748909 2.01564 1.524679 9.271893 854.9093
#>   efficiency      p25      p75
#> 1   427454.7 3.469379 6.265743
#> 
#> $other
#> list()

To register a metric so that it can be called by name, use

registerMetrics(
  list(quartiles = MCMCmetric_quartiles)
)
#> <environment: 0x7fb3a0420560>

Now the name “quartiles” can be included in the metrics argument to compareMCMCs or addMetrics, and the registered MCMCmetric_quartiles will be called.

Converting parameters to match each other across different MCMCs

Sometimes different MCMCs for essentially the same model use different parameter names and/or parameterizations. For example, a common difference in parameterization occurs for the variance of a normal distribution, which might be parameterized as variance, standard deviation, log standard deviation, precision, log precision, and so on. compareMCMCs provides a system to support converting parameters to make them comparable across MCMCs. Like metrics, conversions can be applied as part of compareMCMCs or outside of compareMCMCs. When included in compareMCMCs, they are applied before metrics are calculated.

Example

To continue using the toy example above, let us say that one wants to use log(a) for comparisons, and make that conversion for all MCMCs. Typically one would need different conversions for different MCMCs, but in this case we will apply the same conversion to all.

To make this conversion as part of compareMCMCs, we would do the following.

reparam <- list(log_a  = "log(`a`)", a = NULL)
conversions <- list(nimble = reparam,
                    nimble_slice = reparam,
                    jags = reparam)
res <- compareMCMCs(modelInfo,
                    MCMCs = c(if(this_system_has_rjags) 'jags' else NULL,
                              'nimble', 'nimble_slice'),
                    nimbleMCMCdefs = list(nimble_slice = 'configure_nimble_slice'),
                    conversions = conversions,
                    MCMCcontrol = list(inits = list(a = 1),
                                       niter = 2000,
                                       burnin = 100))
#> building nimble model...
#> Defining model
#>   [Note] Using 'y' (given within 'constants') as data.
#> Building model
#> Setting data and initial values
#> Running calculate on model
#>   [Note] Any error reports that follow may simply reflect missing values in model variables.
#> Checking model sizes and dimensions
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 1
#>    Unobserved stochastic nodes: 1
#>    Total graph size: 5
#> 
#> Initializing model
#> 
#> ===== Monitors =====
#> thin = 1: a
#> ===== Samplers =====
#> RW sampler (1)
#>   - a
#> ===== Monitors =====
#> thin = 1: a
#> ===== Samplers =====
#> slice sampler (1)
#>   - a
#> Compiling
#>   [Note] This may take a minute.
#>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
#> Compiling
#>   [Note] This may take a minute.
#>   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|
#> |-------------|-------------|-------------|-------------|
#> |-------------------------------------------------------|

# We will look at the result using combineMetrics
# rather than generating new html pages.
combineMetrics(res)
#> $byParameter
#>           MCMC Parameter     mean   median        sd  CI95_low CI95_upp
#> 1         jags     log_a 1.505575 1.561667 0.4573951 0.4606434 2.212156
#> 2       nimble     log_a 1.530790 1.597528 0.4551483 0.4551026 2.233067
#> 3 nimble_slice     log_a 1.510952 1.564955 0.4810929 0.3721328 2.246992
#>        ESS efficiency
#> 1 1106.377   276594.3
#> 2  388.788   388788.0
#> 3 2288.069   762689.8
#> 
#> $byMCMC
#>           MCMC min_efficiency mean_efficiency
#> 1         jags       276594.3        276594.3
#> 2       nimble       388788.0        388788.0
#> 3 nimble_slice       762689.8        762689.8

We can see that the parameter is now log_a rather than a. The use of backticks in “log(`a`)” ensures that a is evaluated as a variable and is described more below.

Converting parameters outside of `compareMCMCs`

Say we have run compareMCMCs and only later realized we need to convert parameters. For an example, say we want to reverse the conversion to log_a done above. We can do that as follows. Notice that we need to re-compute metrics after applying the conversions, and we must clear the old metrics before re-computing them.

reparam <- list(a  = "exp(`log_a`)", log_a = NULL)
conversions <- list(nimble = reparam,
                    nimble_slice = reparam,
                    jags = reparam)
applyConversions(res, conversions)
#> $jags
#> <MCMCresult>
#>   Public:
#>     addMetricResult: function (metricResult) 
#>     clearMetrics: function (byParameter = TRUE, byMCMC = TRUE) 
#>     clone: function (deep = FALSE) 
#>     initialize: function (...) 
#>     initializeMetrics: function (silent = FALSE) 
#>     MCMC: jags
#>     metrics: list
#>     rename: function (newName, oldName) 
#>     samples: 4.75304848947591 4.43728310468788 2.99346927645235 3.314 ...
#>     sessionInfo: sessionInfo
#>     setSamples: function (samples) 
#>     times: list
#> 
#> $nimble
#> <MCMCresult>
#>   Public:
#>     addMetricResult: function (metricResult) 
#>     clearMetrics: function (byParameter = TRUE, byMCMC = TRUE) 
#>     clone: function (deep = FALSE) 
#>     initialize: function (...) 
#>     initializeMetrics: function (silent = FALSE) 
#>     MCMC: nimble
#>     metrics: list
#>     rename: function (newName, oldName) 
#>     samples: 10.1532747508493 8.34293485243425 6.9661668837256 6.1664 ...
#>     sessionInfo: sessionInfo
#>     setSamples: function (samples) 
#>     times: list
#> 
#> $nimble_slice
#> <MCMCresult>
#>   Public:
#>     addMetricResult: function (metricResult) 
#>     clearMetrics: function (byParameter = TRUE, byMCMC = TRUE) 
#>     clone: function (deep = FALSE) 
#>     initialize: function (...) 
#>     initializeMetrics: function (silent = FALSE) 
#>     MCMC: nimble_slice
#>     metrics: list
#>     rename: function (newName, oldName) 
#>     samples: 3.69872694323803 3.04326108960467 5.5938645997181 6.3191 ...
#>     sessionInfo: sessionInfo
#>     setSamples: function (samples) 
#>     times: list
clearMetrics(res)
addMetrics(res) # use default metrics
combineMetrics(res) # An easy way to see that it worked
#> $byParameter
#>           MCMC Parameter     mean   median       sd CI95_low CI95_upp       ESS
#> 1         jags         a 4.935454 4.766760 1.951325 1.585117 9.135405 1076.3716
#> 2       nimble         a 5.056351 4.940805 1.979402 1.576335 9.328612  328.2171
#> 3 nimble_slice         a 5.012450 4.782458 2.090743 1.450882 9.459239 1927.1341
#>   efficiency
#> 1   269092.9
#> 2   328217.1
#> 3   642378.0
#> 
#> $byMCMC
#>           MCMC min_efficiency mean_efficiency
#> 1         jags       269092.9        269092.9
#> 2       nimble       328217.1        328217.1
#> 3 nimble_slice       642378.0        642378.0

How to specify conversions

The conversions argument to applyConversions (or to compareMCMCs, which is passed to applyConversions) provides a reasonably flexible system. It should be a named list with names corresponding to MCMCs. Each element is itself of a list, which we call the conversions specification. Each element of the conversions specification defines one conversion, and the elements are processed in order. The name of each element states the name of a column that will be created or modified in the MCMC samples. The element itself can be a character string or call that will be evaluated to generate contents of the column. It can also be NULL or "", in which case the the named column will be removed.

In a character string or call, it can be important to use backticks around other column names, since they often include indexing notation (using “[]”) as part of the name. For example, to create the log of “beta[2]”, you would include a conversion specification list element “log_beta2 = log(`beta[2]`)”. The backticks around “beta[2]” identify “beta[2]” as a name rather than the second element of a vector called “beta”. Doing it this way is necessary because columns in MCMC samples have names like “beta[2]”. In the example of above, we do not really need backticks around “a”, but we have used them for illustration.

In the first example above, the first conversion specification element says that log_a will be calculated as log(a). The second element says that a will be removed.

It is also possible to modify the samples by directly manipulating the samples element of an MCMCresult object. It is simply a matrix.

Writing new page component plugins.

Writing a plugin to provide a new page component, such as a new figure, to be drawn by make_MCMC_comparison_pages is more complicated than writing a new metric. A page component plugin is a list with up to five elements:

make is the character name of a function to create the plottable output such as a ggplot object. This is described more below.
fileSuffix is a character suffix to be pasted to the model name to form a filename to which the figure will be output. For example, fileSuffix = '_mySuffix' will lead to an output file called “model_mySuffix.jpg”, if the model is named “model”.
linkText is the character text that will be a hyperlink at the top of the comparison page.
plot is the name of a function to call to plot the output to a jpeg file such as “model_mySuffix.jpg”. The function takes as input the plottable element of the list returned from make. If no plot function name is provided, the default function plot will be used.
control is a list that will be passed to the make function. This allows each plugin to offer users to pass arbitrary information to make.

More about the `make` function.

The make element of a plugin names a function. This function will be called with two arguments. The first is a tidy-formatted combination of metrics from the results list returned by compareMCMCs (or concatenations of such lists). This is created by combineMetrics(results, include_times=TRUE). The second is the control element of the plugin.

The function can return information for either a text component or figure component to the final html. For a figure component, it should return a list with three elements:

plottable should be an object in one of two formats. If the plugin provides a plot element, then plottable can be any object. It will simply be passed to the function named by plot, which should generate a figure. If the plugin does not provide a plot element, then plottable should be an object such that calling plot() with that object as the argument results in generating a plot. For example, an object returned from ggplot will work.
height should be the height in inches to be passed as the height argument to jpeg.
width should be the width in inches to be passed as the width argument to jpeg.

The plotting will be done between a call to jpeg and a corresponding call to dev.off().

For a text component, it should return a list with element:

printable should be either a character string of arbitrary html or an object of class xtable from the xtable package. The latter will be printed with type = 'html'.

Format of combined metrics provided to the `make` function.

An example of the combined metrics from the above example is:

combineMetrics(res, include_times = TRUE)
#> $byParameter
#>           MCMC Parameter     mean   median       sd CI95_low CI95_upp       ESS
#> 1         jags         a 4.935454 4.766760 1.951325 1.585117 9.135405 1076.3716
#> 2       nimble         a 5.056351 4.940805 1.979402 1.576335 9.328612  328.2171
#> 3 nimble_slice         a 5.012450 4.782458 2.090743 1.450882 9.459239 1927.1341
#>   efficiency
#> 1   269092.9
#> 2   328217.1
#> 3   642378.0
#> 
#> $byMCMC
#>           MCMC min_efficiency mean_efficiency
#> 1         jags       269092.9        269092.9
#> 2       nimble       328217.1        328217.1
#> 3 nimble_slice       642378.0        642378.0
#> 
#> $times
#>               burnin post-burnin total sampling
#> jags         0.00000     0.00400          0.004
#> nimble       0.00005     0.00095          0.001
#> nimble_slice 0.00015     0.00285          0.003

Using a new page component

Before using a new page component, you must register it using:

registerPageComponents(
  list(myNewComponent = 
         list(make = "myMakeFunction",
              fileSuffix = "_myPageComponent",
              linkText = "My new page component.")
       )
  )
#> <environment: 0x7fb3a0443718>

Then the name “myNewComponent” can be included in the character vector argument pageComponents to make_MCMC_comparison_pages.

Pre-defined page components and defaults

The pre-defined page components include:

timing: table of setup, burnin, post-burnin, and total sampling (burnin + post-burnin) times.
efficiencySummary: barplot with min and mean MCMC efficiencies.
efficiencySummaryAllParams: efficiencySummary with additional figure of parameter-wise MCMC efficiencies.
paceSummaryAllParams: max, mean, and parameter-wise MCMC paces (pace = 1/efficiency).
efficiencyDetails: parameter-by-parameter barplots of efficiencies.
posteriorSummary: parameter-by-parameter plots of posterior summaries.

If the pageComponents argument to make_MCMC_comparison_pages is NULL, all of these except efficiencySummary (which is included in efficiencySummaryAllParams) will be included.

Writing new MCMC package plugins.

An interface to a new MCMC engine is provided as a function. It accepts the following inputs:

MCMCinfo: The element of externalMCMCinfo named to match this MCMC plugin. This can contain whatever information is needed for the plugin.
MCMCcontrol: The MCMCcontrol argument to compareMCMCs. The seed argument to compareMCMCs will be added to this list with the name seed. Prior to calling the MCMC engine plugin, set.seed() will be called with this seed as its argument.
monitorInfo: A list of names of parameters to monitor (record) in MCMC output. These are provided in two formats. monitorInfo$monitorVars is a character vector of the names of variables for which every element should be monitored. monitorInfo$monitors is a character vector with the name of each scalar element of variables in monitorInfo$monitorVars E.g., if there is a variable x with elements x[1]- x[10], then “x” will be in monitorInfo$monitorVars and “x[1]”, “x[2]”, …, “x[10]” will be in monitorInfo$monitors. Both formats may be needed.
modelInfo : The modelInfo argument to compareMCMCs. If the call to compareMCMCs involved creating a nimble model, it will be added to this list with the name model.

The function should use these arguments to run its MCMC engine and return an MCMCresult object with times (a list with elements setup, burnin, postburnin, and sampling, each with one number) and samples (a matrix with named columns for monitored parameters and rows for MCMC iterations) populated. If the MCMC element is populated with a character string, this will be used as the label of the MCMC in metrics and comparison pages. Otherwise the label used in the MCMCs argument to compareMCMCs will be used as the MCMC element.

A new MCMC engine can be registered with registerMCMCengine(name, fun) where fun is the function just described and name is a name for it (which will be recognized in the MCMCs argument to compareMCMCs).

Overview of compareMCMCs