---
title: "Computation framework"
output: rmarkdown::html_vignette
bibliography: ../inst/REFERENCES.bib
vignette: >
  %\VignetteIndexEntry{Computation framework}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## JAGS vs. NIMBLE vs. Stan

We use both JAGS [@plummer_2003] and NIMBLE [@valpine_et_al_2017]. NIMBLE uses a derivative of the [JAGS/BUGS language](https://r-nimble.org/manual/cha-writing-models.html#sec:supp-feat-bugs), so that the models can be used interchangeably with only minor modifications. Still, JAGS and NIMBLE are conceptually quite different. While it is theoretically possible to extend the functionality of JAGS by custom functions, this process is not straightforward and requires advanced programming skills. Therefore, the available functions and samplers have to be taken as they are. In contrast, NIMBLE offers a great suite of customization possibilities. This is especially true for the samplers used even for individual nodes. We tried to tune the model parameters, but with only limited success. We found that the default samplers of NIMBLE exhibit poorer mixing than those of JAGS and thus need longer chains to reach the same level in terms of quality measures. Overall, we observed that the same model with the same data takes longer to run in NIMBLE than in JAGS. Additionally, NIMBLE is less stable than JAGS, especially in the multicore setup.

For the reasonse given above, we propose to use JAGS for the univariate cases and cases with multiple normal likelihood, while for multivariate applications with multinormal likelihood function, NIMBLE is currently the only viable option.

Despite the fact that the models of JAGS and NIMBLE are almost identical, the differences are subtle enough so that the same data and parameters will produce slightly different results. This is because JAGS and NIMBLE use different samplers and the way the ordinal thresholds are ordered works differently.

We interface JAGS with a wrapper function which in turn relies on `runjags` [@derwood_2016] to run JAGS. NIMBLE is a package in its own right [@valpine_et_al_2017]. The processing of the chains is done with the `coda` package [@plummer_et_al_2006]. For further analysis, `baytaAAR` provides custom functions but in the vignettes we also show how to analyze the data with the packages `bayesplot` [@gabry_et_al_2019] and `tidybayes` [@kay_2024].

We do not use [`Stan`](https://mc-stan.org) which is currently probably the most popular choice when it comes to working with Bayesian models in R. There are two main reasons for this: Firstly, `Stan` is much more restrictive with respect to unobserved nodes. While JAGS and NIMBLE simply treat those as parameters to be estimated, they have to be explicitly declared in `Stan`. Secondly, and more importantly, `Stan` does not provide the interval-distribution (`dinterval`) which is instrumental to our application because only with the `dinterval`-distribution is it possible to implement the multivariate model.


## Parallel processing and number of chains

On modern hardware with usually more than one core, parallel processing is very useful as it cuts down time significantly. In most cases, you would use as many cores as chains. The default for JAGS and NIMBLE is three chains. More chains can help in reaching quality criteria like `PSRF` faster. However, it also increases the size of the MCMC samples so that post-processing can become limited by available memory. Therefore, and following J. Kruschke [-@ref_77051], we control the number of saved iterations by dividing `numSavedSteps` by the number of chains. So, increasing the number of chains will automatically shorten the number of saved iterations which will usually lead to poorer mixing. In this respect, we observed that a higher number of chains does not outweigh the disadvantage of shorter chains. We therefore do not recommend increasing the number of chains.

Parallel processing is straightforward in `bay.ta()`. You simply specify the mode of operation with the parameter `multicore` which can be `FALSE` (single-core-processing) or `TRUE` (parallel processing). This runs out-of-the-box with JAGS but parallel processing in NIMBLE is much more complicated. The [NIMBLE manual](https://r-nimble.org/examples/parallelizing_NIMBLE.html) gives a general account of the necessary steps. For a package like ours, there is the additional difficulty that we provide custom functions which have to be supplied in the environment of each cluster. However, simply copying the functions to the global environment would be unacceptable behaviour of a package. We use the R-package `parallel`, which is part of base R, for multi-core processing in NIMBLE.

In the end, we hope that we succeeded to make parallel processing with NIMBLE as unobtrusive as possible by providing a wrapper function which takes care of the different modes of operation. Please note, though, that you will get no feedback while your model is running in parallel. Therefore, it is good practice to start with a low number of steps and see how long these take, 'guesstimate' the required number of steps from the quality measures (especially `ESS`) and change the settings accordingly.

Generally, JAGS is very stable in parallel processing so apart from your fans possibly spinning up, you will notice very little difference between single- and multicore-processing. For NIMBLE, this is unfortunately not always the case. Apart from the lack of feedback, there is also the possibility of a memory leak. This makes the cores consume more and more memory until the machine becomes unusable and might eventually crash. Therefore, we would strongly advise to monitor the cores closely to see if their memory consumption remains stable or steadily increases. If the latter is the case, the model should be stopped and the R session should be restarted. Furthermore, it might be necessary to terminate the processes manually using your system monitor.

---

## References