Can I parallelize a Bayesian hierarchical model using MCMC simulation (specifically, w/the R package ConStruct)?

I'm very new to Bayesian models and MCMC, so please bear with me while I try to explain.

I'm trying to use the R package "conStruct", which provides "A method for modeling genetic data as a combination of discrete layers, within each of which relatedness may decay continuously with geographic distance." The analyses of this package are implemented in rstan.

I'm aiming to parallelize the following code. However, my main issue is that, from trying to understand posts like this and this, parallelization is difficult when these chains have to be run sequentially. Is this possible? And if so, how would I do so?

Importantly, here are the definitions for n.chains and n.iter (the other arguments aren't as important to understand in this context):

n.chains: An integer indicating the number of MCMC chains to be run in the analysis. n.iter: An integer giving the number of iterations each MCMC chain is run. If the number of iterations is greater than 500, the MCMC is thinned so that the number of retained iterations is 500 (before burn-in).

conStruct(spatial = FALSE, 
                    K = 1, 
                    freqs = allele_frqs, 
                    geoDist = NULL, 
                    coords = coords,
                    prefix = "nsp_K1_iter5000_chains3",
                    n.chains = 3,
                    n.iter = 5000)

Alternatively, could I run each chain at a time, but on different cores? How would I change my code to do that?

Thanks so much!

Solution

Subdividing a single chain to run across multiple cores is almost certainly not possible at present (I think it would be theoretically possible to implement in Stan, the engine underneath rstan/ConStruct, but very tricky ...)

However, it is easy and effective to run multiple chains in parallel, for at least a modest speedup. Assuming that the ConStruct package isn't doing anything clever under the hood, setting

options(mc.cores = 3)

should make each of your three chains run on a separate core (assuming you have at least 3 cores on your machine, of course), giving you approximately a three-fold speedup (assuming that you have enough memory to run three instances of your chain at once ...)

(@Axeman comments that

conStruct() passes ... to rstan::sampling, so [you] should be able to pass cores = 3 to conStruct()

)

The startup message for rstan recommends options(mc.cores = parallel::detectCores()), but this will only make a difference if you run more than three chains.

There is a fairly complicated tradeoff as to whether it's worth running more, shorter chains on more cores (the burn-in/adaptation time that's done on every core will effectively be wasted).

(In addition to performance benefits of running multiple chains, the most reliable MCMC convergence diagnostics rely on running multiple chains; I personally think the package author underemphasizes this in the vignette. See e.g. the "diagnose" section of this introduction; the functions shown there should probably work fine with your fitted models ...)