Add `nchains` option to `mcmc`

[danielinteractive]

So far in https://github.com/Roche/crmPack/blob/main/R/mcmc.R#L75 we are not using the nchains argument of rjags::jags.model, therefore we can only run a single MCMC chain at the moment.

We will need some downstream updates when adding the nchains functionality, where maybe we can borrow inspiration from the way rstan or cmdstanr handles this.

To do:

• Design doc
• Separate small issues and PRs as much as possible

Ideas:

• Put chains as the last dimension of Samples and derived objects
  • e.g. a matrix with nsamples rows and nchains columns
• Have a collapse() method on the Samples objects to collapse all chains into a single one upon request
• For plots, by default if there is more than one chain present, use faceting across chains to visualize the different results and compare them
  • e.g. multiple dose-response curves with one set of median and credible interval lines for each chain
• adapt the get() method to transform correctly into ggmcmc object representing multiple chains
• adapt the tidy() methods to convert multiple chains to tidy tibbles

[Puzzled-Face]

I fully support this idea. It will allow us to add more diagnostics to confirm that the chains have converged correctly.

[Puzzled-Face]

My first impression is that a comprehensive solution will be more difficult than it might at first seem. Adding the option to run more than one chain - and adapting downstream operations to deal with multiple chains - should be relatively straightforward. The greater difficulty may well be in in initialising the chains appropriately.

Because of the way JAGS works, it may not be necessary to introduce a new parameter. JAGS will run one chain for each list of initialisation parameters it is passed. At the moment, crmPack passes a list with one element for each parameter that JAGS needs. So if we pass a list of such lists as the init parameter of jags.model, then we will get a different chain of samples for each element of the outer list.

To test for convergence of the MCMC process robustly, different chains need to have different initial values. Ideally, the starting values should be in some way "extreme" to give credence to the belief that the chain will always converge (and to a constant value), regardless of how it is initialised. Thus, the initial values of a chain's parameters are model dependent, and should be different even for chains that are fitting the same model.

I cannot yet think of a simple way to do this, let alone to do this in a manner that does not break existing code. I believe this problem has two main parts:

• The fact that crmPack uses S4 objects
• The way in which crmPack handles the interface between R and the JAGS executable

Currently, JAGS obtains initial values for model parameters via the init slot of a Model object. For example, here is the init slot for LogisticLogNormal:

init = function() {
  list(theta = c(0, 1))
}

So the same initial values are used for all LogisticLogNormal fits, regardless of the model's prior. The same is true for all other Model classes.

So, first we need to be able to generate starting values that are "extreme" relative the prior defined in this specific Model object. Initially, I thought we could do this using the same object's prior model slot. But this is not possible for two reasons:

1. There appears to be no way for one slot of an S4 object to access another slot of the same object. S4 has no concept of self:: or $this.
2. The priormodel slot contains pseudocode. It's JAGS code, not R code. A call similar to model@priormodel() with throw an exception, even if the model is correctly instantiated. For example:

my_model <- LogisticLogNormal(
  mean = c(-0.85, 1),
  cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2),
  ref_dose = 56
)

my_model@priormodel()
Error in my_model@priormodel() : object 'theta' not found

In this case, the priormodel slot is

priormodel = function() {
   theta ~ dmnorm(mean, prec)
   alpha0 <- theta[1]
   alpha1 <- exp(theta[2])
 }

theta ~ dmnorm(mean, prec) is not syntactically valid R code and, even if it were, mean and prec are not defined in scope.

One workaround might be to introduce a new slot chain_initialiser or similar, which would be a function that could return appropriate initial values for the particular instance of the model. It could pick up appropriate parameter values from the parameters passed to the Model's constructor.

The issue then would be to ensure consistency following subsequent changes to the model's prior.

I will think some more on this and start drafting a design document once #829 is complete.