Bayesian fitting packages

[adamkucharski]

Do people have a preferred package when running MCMC (e.g. Metropolis Hastings) on a generic R function that takes a parameter vector and data.frame fun(param, data), then returns a likelihood (e.g. like the likelihood functions in epichains, or a short script based on something else). So rather than a package like rstan, rbi or pomp that has custom syntax for the model, something that takes an existing function as input.

Have had a couple of situations where it's come up for collaborative projects where Epiverse functionality could be helpful, but thought it might be worth some alignment if these could potentially become vignettes at a later date.

A few options I've come across:

• MCMCpack. Lots of functions and contributors and works with coda functions. Some example code below, although can't see capacity for adaptive MCMC.

• adaptMCMC. One core adaptive MCMC function.

• LaplacesDemon. Widely used, but current maintenance structure unclear.

There's also some related discussion about posteriors in this quickfit issue.

# Example using MCMCmetrop1R to fit a distribution to simulated data
library(MCMCpack)

data_in <- rgamma(100, shape = 2, rate = 0.5) # Simulate data

f <- function(param,data) {
  if (any(param <= 0)) return(-Inf) # Ensure positive parameters
# Note: to implement properly, would need to include transformation (e.g. lognormal and effect on transition probabilities – this example is just to get some code working)
  
  # Extract parameters
  alpha <- param[1]
  beta <- param[2]
  
  # Assuming non-informative priors for alpha and beta
  log_prior <- 0 # Could add informative priors here
  
  # Calculate log-likelihood
  log_likelihood <- sum(dgamma(data, shape = alpha, rate = beta, log = TRUE))
  
  # Return log-posterior (log-likelihood + log-prior)
  return(log_likelihood + log_prior)
}

initial_param <- c(1, 1) # Initial guess for shape and rate

output <- MCMCmetrop1R(fun = f, 
                       theta.init = initial_param,
                       mcmc = 10000, # Number of MCMC iterations
                       burnin = 1000, # Burn-in period
                       verbose = TRUE,
                       data = data_in) # Add data argument

plot(output)

[sbfnk]

I have found LaplacesDemon useful in the past but agree there doesn't seem to be a huge amount of development.

[Bisaloo]

With kind of the same interface as adaptMCMC (same original author), I'm maintaining:

• https://cran.r-project.org/package=simpleMH; https://hugogruson.fr/simpleMH/
• https://cloud.r-project.org/package=mcmcensemble; https://hugogruson.fr/mcmcensemble/

Both work with coda as well.

[adamkucharski]

Thanks, will have a look. @dchodge and @thimotei also mentioned the following:

• LazyMCMC
• DrJacoby
• Bayesian Tools
• greta

Plus ptmc by @dchodge, which can also handle discrete spaces.

Might be worth a head-to-head at some point, to identify key benefits/limitations - although I suppose if we're using for simple-ish use cases, shouldn't be too hard to switch between them if they've all got a reasonably plug-and-play interface.

[thimotei]

Thanks @adamkucharski. Although I think this might be the wrong Github user...mine is thimotei!

[thimotei]

Out of those I've used, I'd say greta is probably the best bet as it lets you use any R function for the likelihood, its simple to read and write but its fast under the hood. Only downside is that it can be tricky to install, given that it depends on tensorflow under the hood (from a Python virtual environment within R)

[Bisaloo]

although I suppose if we're using for simple-ish use cases, shouldn't be too hard to switch between them if they've all got a reasonably plug-and-play interface.

I think this is key and where we are trying to get with Epiverse (slides on interoperability). By providing standard interfaces, with standard R objects (e.g., R functions), we allow flexibility for users and let use quickly switch between alternatives depending on their exact requirements, without having to rewrite their entire code.

The MCMC package ecosystem is a good model to try and replicate where possible

[adamkucharski]

Thanks all. Going to explore some comparisons in a draft .Rmd, partly just to get my head around what's available, with potential to turn into a vignette as a future advanced use case in Epiverse. Maybe using {epichains} likelihoods as a starting point as have quite a bit of cluster size distribution data floating around from previous projects/talks (unless you've already got something in progress @sbfnk @jamesmbaazam?)

One thing I've noticed from an initial look at packages is that most put quite a simple adjustment for parameter positivity (i.e. set likelihood = -Inf if a sampled parameter is below zero, which can create some mixing inefficiencies compared to a transform-based approach (e.g. like the below). But perhaps a simple step-function prior is OK for a simple walkthrough example.

params <- exp(log_params)  # Transform parameters to original scale
jacobian_correction <- sum(log_params)  # Diagonal Jacobian with log transform
original_lik <- lik_function(params, data)
adjusted_lik <- original_lik + jacobian_correction # Adjusted likelihood step
return(adjusted_lik)

[sbfnk]

Yes, good point - alternative for MCMC is to use a truncated proposal distribution

[adamkucharski]

I've started a draft vignette on {epichains} to explore and compare options. {MCMCpack} and {BayesianTools} were pretty plug-and-play, but didn't manage to get the formatting of inputs right for {mcmcensemble} yet @Bisaloo. Also haven't got all the dependencies lined up for {greta}, so this also disabled in current .Rmd.

Would be interested to hear thoughts/comments!

[Bisaloo]

I'll have a look but on a side note, I don't think this type of documentation that doesn't directly talk about the package should be a vignette. It's quite complex to get started with our packages and it's good that the key information is very clearly signposted: "here is what epichains does, here are the functions, and how to use them".

The nice to have extras, related topics, advanced asides, should probably be somewhere else to keep the package documentation on topic. For this specific task, I could see it well as a blog post, or as a training material somewhere.

[adamkucharski]

Yep, epichains not necessarily the optimal final home for this. Although given the package includes likelihoods, it does seem useful to have a vignette to show user how to fit with these likelihoods (e.g. syntax required to link with data and input parameters etc). Perhaps just using one fitting package, so most of vignette is focused on epichains functionality?

Can definitely fitting comparisons consider for future training materials - but I think there's also value in having shorter quick start type scripts for users who are more familiar with broad methods, but less so the package specific inputs and outputs. Happy to move over the fuller comparison to a blog post once get mcmcensemble and Greta working - can then reference in relevant places in other packages.

[adamkucharski]

Having discussed with @avallecam, we're going to make the howto repo a home for most of these examples, with some more applied examples being converted into case studies.

[Bisaloo]

@adamkucharski, thanks for the report on mcmcensemble. By design, the method only makes sense if you have more than 1 walker. I will update the package to return a clearer error message.

Here the full corrected chunk:

# Number of ensembles
n_walkers <- 2

# Define random initial values for ensemble estimation
init_ensemble <- data.frame(
                    "R" = rnorm(n_walkers, init_param[["R"]], 0.1),
                    "k" = rnorm(n_walkers, init_param[["k"]], 0.1)
                  )

start_time <- Sys.time()

# Run MCMC to estimate parameters
result_mcmcensemble <- MCMCEnsemble(lik_function, 
                                    inits = init_ensemble,
                                    max.iter = n_iter, 
                                    n.walkers = n_walkers, # ensemble size
                                    method = "stretch", # proposal method
                                    coda = TRUE)

# Store time and effective sample size
time_mcmcensemble <- Sys.time() - start_time 
ess_mcmcensemble <- effectiveSize(result_mcmcensemble$samples)

# Plot posterior estimates
plot(result_mcmcensemble$samples)