---
title: "Overview of compareMCMCs"
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Overview of compareMCMCs}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
## Note: special steps are needed to build this vignette.
## See README.txt.
##
## If rstan generates C++ output with an error compiling foo.c,
## try re-installing rstan. Clues suggest that resulting updates to StanHeaders
## seems to fix the problem.
knitr::opts_chunk$set(
  echo = TRUE,
  collapse = TRUE,
  comment = "#>"
)
library(compareMCMCs)
this_system_has_rjags <- requireNamespace("rjags")
if(!this_system_has_rjags) message("Portions of this vignette use package rjags.  That is not installed, so those portions will be skipped.")
this_system_has_rstan <- requireNamespace("rstan")
if(!this_system_has_rstan) message("Portions of this vignette use package rstan.  That is not installed, so those portions will be skipped.")
work_dir <- tempdir() # This line is the default to respect CRAN policy that a vignette should not generate local files.
work_dir <- getwd() # This is what you should use to actually build the vignette and have the internal links work. By uncommenting this line, you are confirming that you accept having this script write jpg and html files locally (besides the knitr output of this file itself).
```

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

```{r, message=FALSE}
# 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.
)
file.remove(work_dir, "toy.stan")
make_MCMC_comparison_pages(res, dir = work_dir,  modelName = 'example1')
```

Results are [here](example1.html).  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:

1. 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.)
2. 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:
    i. 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).
    ii. The character *name* of a function as described in (i).
    iii. 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:

```{r, eval=!this_system_has_rjags, echo=FALSE}
message("Skipping this code because rjags is not installed.")
```

```{r, eval=this_system_has_rjags}
res_jags <- compareMCMCs(modelInfo,
                MCMCs = c('jags'),
                MCMCcontrol = list(inits = list(a = 1),
                                   niter = 2000,
                                   burnin = 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))
res <- c(res_jags, res_nimble)
make_MCMC_comparison_pages(res, dir = work_dir, modelName = 'example2')
```

Results are [here](example2.html). 

# 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:
```{r}
res$nimble$metrics
```

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:
```{r, eval = FALSE}
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:

1. By calling `addMetrics` with a list of metrics, which can be functions or registered metric names.
2. 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:

```{r}
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:

```{r}
addMetrics(res, list(MCMCmetric_quartiles))
res$nimble$metrics
```

To register a metric so that it can be called by name, use
```{r}
registerMetrics(
  list(quartiles = MCMCmetric_quartiles)
)
```

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.

```{r}
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))

# We will look at the result using combineMetrics
# rather than generating new html pages.
combineMetrics(res)
```

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.*

```{r}
reparam <- list(a  = "exp(`log_a`)", log_a = NULL)
conversions <- list(nimble = reparam,
                    nimble_slice = reparam,
                    jags = reparam)
applyConversions(res, conversions)
clearMetrics(res)
addMetrics(res) # use default metrics
combineMetrics(res) # An easy way to see that it worked
```

## 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:
```{r}
combineMetrics(res, include_times = TRUE)
```

### Using a new page component

Before using a new page component, you must register it using:
```{r}
registerPageComponents(
  list(myNewComponent = 
         list(make = "myMakeFunction",
              fileSuffix = "_myPageComponent",
              linkText = "My new page component.")
       )
  )
```

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`).
