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:
compareMCMCs function to run one or more MCMCs and
manage the results;MCMCresult class to manage results by storing
samples, timing information, and metrics or summaries of
performance;If you want to run the code below yourself, first set the
work_dir variable using one of the following options. This
will be the directory where html and jpg files generated in the exampes
below will be placed.
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 stan_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
#>
#>
|
| | 0%
|
|+ | 2%
|
|++ | 4%
|
|+++ | 6%
|
|++++ | 8%
|
|+++++ | 10%
|
|++++++ | 12%
|
|+++++++ | 14%
|
|++++++++ | 16%
|
|+++++++++ | 18%
|
|++++++++++ | 20%
|
|+++++++++++ | 22%
|
|++++++++++++ | 24%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 28%
|
|+++++++++++++++ | 30%
|
|++++++++++++++++ | 32%
|
|+++++++++++++++++ | 34%
|
|++++++++++++++++++ | 36%
|
|+++++++++++++++++++ | 38%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 42%
|
|++++++++++++++++++++++ | 44%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 48%
|
|+++++++++++++++++++++++++ | 50%
|
|++++++++++++++++++++++++++ | 52%
|
|+++++++++++++++++++++++++++ | 54%
|
|++++++++++++++++++++++++++++ | 56%
|
|+++++++++++++++++++++++++++++ | 58%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 62%
|
|++++++++++++++++++++++++++++++++ | 64%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 68%
|
|+++++++++++++++++++++++++++++++++++ | 70%
|
|++++++++++++++++++++++++++++++++++++ | 72%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|++++++++++++++++++++++++++++++++++++++ | 76%
|
|+++++++++++++++++++++++++++++++++++++++ | 78%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 82%
|
|++++++++++++++++++++++++++++++++++++++++++ | 84%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 88%
|
|+++++++++++++++++++++++++++++++++++++++++++++ | 90%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 92%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
#>
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
#>
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
#>
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 8e-06 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.08 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.03 seconds (Warm-up)
#> Chain 1: 0.141 seconds (Sampling)
#> Chain 1: 0.171 seconds (Total)
#> Chain 1:
#> Warning: There were 94 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(file.path(work_dir, "toy.stan"))
#> [1] TRUE
make_MCMC_comparison_pages(res, dir = file.path(work_dir, "example1"), modelName = 'example1',
control = list(res = 75))
# Open "example1.html" from a browser, or do
# browseURL(file.path(work_dir, "example1", example1.html")) Results (run on one of the author’s machines) 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.
Default res (passed to function jpeg) is
300 but is reduced here to reduce the size of this vignette.
If one wants to get hold of the figures used in the html output, they
are returned as a list of ggplot objects (using package
ggplot2) from make_MCMC_comparison_pages.
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 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.
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.
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:
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.)nimbleMCMCdefs list with the
same name as entered in MCMCs (e.g., “nimble_slice”). That
element can take one of three formats:
configureMCMC, possibly modified using its
addSampler and removeSamplers methods).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.So far, we have assumed the inputs to compareMCMCs to be
used by nimble are the ingredients from which
compareMCMCs will build a nimble model and any MCMCs,
compile them all, and run the compiled MCMCs. It is also possible to
build the model yourself (using nimbleModel) and provide
that (named “model”) in the modelInfo list.
Furthermore, it is possible to compile the model yourself and provide
the compiled version as “model” in the modelInfo list. In
that case, you must also compile any MCMCs yourself and provide them in
the nimbleMCMCdefs list. The steps to do this (modified
from above) are:
nimbleMCMCdefs list with the
same name as entered in MCMCs. That element should be the
compiled MCMC object, such as returned by
compileNimble(MCMC, project = model) where
MCMC is the result of buildMCMC,
model is the result of nimbleModel, and a
compiled version of model has already been created.The package comes with the following configurations, which are
similar to options in nimble’s configureMCMC function.
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
#>
#>
|
| | 0%
|
|+ | 2%
|
|++ | 4%
|
|+++ | 6%
|
|++++ | 8%
|
|+++++ | 10%
|
|++++++ | 12%
|
|+++++++ | 14%
|
|++++++++ | 16%
|
|+++++++++ | 18%
|
|++++++++++ | 20%
|
|+++++++++++ | 22%
|
|++++++++++++ | 24%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 28%
|
|+++++++++++++++ | 30%
|
|++++++++++++++++ | 32%
|
|+++++++++++++++++ | 34%
|
|++++++++++++++++++ | 36%
|
|+++++++++++++++++++ | 38%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 42%
|
|++++++++++++++++++++++ | 44%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 48%
|
|+++++++++++++++++++++++++ | 50%
|
|++++++++++++++++++++++++++ | 52%
|
|+++++++++++++++++++++++++++ | 54%
|
|++++++++++++++++++++++++++++ | 56%
|
|+++++++++++++++++++++++++++++ | 58%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 62%
|
|++++++++++++++++++++++++++++++++ | 64%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 68%
|
|+++++++++++++++++++++++++++++++++++ | 70%
|
|++++++++++++++++++++++++++++++++++++ | 72%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|++++++++++++++++++++++++++++++++++++++ | 76%
|
|+++++++++++++++++++++++++++++++++++++++ | 78%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 82%
|
|++++++++++++++++++++++++++++++++++++++++++ | 84%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 88%
|
|+++++++++++++++++++++++++++++++++++++++++++++ | 90%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 92%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
#>
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
#>
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
# 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 = file.path(work_dir, "example2"), modelName = 'example2',
control = list(res = 75))Results are here.
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 419122.7 419122.7
#>
#> $byParameter
#> MCMC Parameter mean median sd CI95_low CI95_upp ESS
#> 1 nimble a 4.890002 4.696869 1.902419 1.66887 9.0571 838.2455
#> efficiency
#> 1 419122.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:
addMetrics with a list of metrics, which can
be functions or registered metric names.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 419122.7 419122.7 2.626136
#>
#> $byParameter
#> MCMC Parameter mean median sd CI95_low CI95_upp ESS
#> 1 nimble a 4.890002 4.696869 1.902419 1.66887 9.0571 838.2455
#> efficiency p25 p75
#> 1 419122.7 3.512159 6.138295
#>
#> $other
#> list()To register a metric so that it can be called by name, use
Now the name “quartiles” can be included in the metrics
argument to compareMCMCs or addMetrics, and
the registered MCMCmetric_quartiles will be called.
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.
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
#>
#>
|
| | 0%
|
|+ | 2%
|
|++ | 4%
|
|+++ | 6%
|
|++++ | 8%
|
|+++++ | 10%
|
|++++++ | 12%
|
|+++++++ | 14%
|
|++++++++ | 16%
|
|+++++++++ | 18%
|
|++++++++++ | 20%
|
|+++++++++++ | 22%
|
|++++++++++++ | 24%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 28%
|
|+++++++++++++++ | 30%
|
|++++++++++++++++ | 32%
|
|+++++++++++++++++ | 34%
|
|++++++++++++++++++ | 36%
|
|+++++++++++++++++++ | 38%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 42%
|
|++++++++++++++++++++++ | 44%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 48%
|
|+++++++++++++++++++++++++ | 50%
|
|++++++++++++++++++++++++++ | 52%
|
|+++++++++++++++++++++++++++ | 54%
|
|++++++++++++++++++++++++++++ | 56%
|
|+++++++++++++++++++++++++++++ | 58%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 62%
|
|++++++++++++++++++++++++++++++++ | 64%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 68%
|
|+++++++++++++++++++++++++++++++++++ | 70%
|
|++++++++++++++++++++++++++++++++++++ | 72%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|++++++++++++++++++++++++++++++++++++++ | 76%
|
|+++++++++++++++++++++++++++++++++++++++ | 78%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 82%
|
|++++++++++++++++++++++++++++++++++++++++++ | 84%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 88%
|
|+++++++++++++++++++++++++++++++++++++++++++++ | 90%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 92%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
#>
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
#>
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
#> ===== 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.538417 1.593251 0.4710176 0.4777214 2.297639
#> 2 nimble log_a 1.517911 1.590057 0.4662310 0.3734571 2.245108
#> 3 nimble_slice log_a 1.503304 1.561714 0.4721760 0.4205572 2.237085
#> ESS efficiency
#> 1 1354.0812 135408.1
#> 2 411.9158 205957.9
#> 3 2293.5906 573397.6
#>
#> $byMCMC
#> MCMC min_efficiency mean_efficiency
#> 1 jags 135408.1 135408.1
#> 2 nimble 205957.9 205957.9
#> 3 nimble_slice 573397.6 573397.6We 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.
compareMCMCsSay 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: 2.67952055004207 2.68620692599888 6.22924980096044 7.864 ...
#> 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: 4.42247227738948 5.15937649146681 5.94280698562883 5.460 ...
#> 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: 5.56368011701955 2.2550413636501 2.46604884066722 8.0955 ...
#> 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 5.131928 4.919719 2.110617 1.612406 9.950711 1088.4257
#> 2 nimble a 5.015437 4.904027 2.010647 1.452849 9.441434 390.7049
#> 3 nimble_slice a 4.955322 4.766984 2.043282 1.522848 9.366010 1913.9187
#> efficiency
#> 1 108842.6
#> 2 195352.4
#> 3 478479.7
#>
#> $byMCMC
#> MCMC min_efficiency mean_efficiency
#> 1 jags 108842.6 108842.6
#> 2 nimble 195352.4 195352.4
#> 3 nimble_slice 478479.7 478479.7The 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 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.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'.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 5.131928 4.919719 2.110617 1.612406 9.950711 1088.4257
#> 2 nimble a 5.015437 4.904027 2.010647 1.452849 9.441434 390.7049
#> 3 nimble_slice a 4.955322 4.766984 2.043282 1.522848 9.366010 1913.9187
#> efficiency
#> 1 108842.6
#> 2 195352.4
#> 3 478479.7
#>
#> $byMCMC
#> MCMC min_efficiency mean_efficiency
#> 1 jags 108842.6 108842.6
#> 2 nimble 195352.4 195352.4
#> 3 nimble_slice 478479.7 478479.7
#>
#> $times
#> burnin post-burnin total sampling
#> jags 3e-03 0.0070 0.010
#> nimble 1e-04 0.0019 0.002
#> nimble_slice 2e-04 0.0038 0.004Before using a new page component, you must register it using:
registerPageComponents(
list(myNewComponent =
list(make = "myMakeFunction",
fileSuffix = "_myPageComponent",
linkText = "My new page component.")
)
)
#> <environment: 0x7feb1fae9288>Then the name “myNewComponent” can be included in the character
vector argument pageComponents to
make_MCMC_comparison_pages.
The pre-defined page components include:
If the pageComponents argument to
make_MCMC_comparison_pages is NULL, all of
these except efficiencySummary (which is included in
efficiencySummaryAllParams) will be included.
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).
nimble, jags, and stan.compareMCMCs.