Batched Resilient and Rapid Estimation Workflow through Approximate Bayesian Computation • BRREWABC

BRREWABC (Batched Resilient and Rapid Estimation Workflow through Approximate Bayesian Computation, pronounced “brew abc”: /bruː ˌeɪ.biːˈsiː/) : an R package designed to facilitate inference through parallelized Approximate Bayesian Computation algorithms (Sequential Monte Carlo (ABC SMC), rejection). This package streamlines the process of conducting Bayesian inference for complex models by implementing efficient parallelization techniques.

Overview

ABC-rejection

The ABC rejection algorithm approximates the posterior distribution by sampling parameters from a prior, simulating data, and then calculating a distance between the summary statistics of the simulated and observed data. If this distance falls within a predefined tolerance, the parameters are accepted. This process is repeated until a sufficient number of parameter sets are accepted, approximating the posterior distribution of the parameters given the observed data.

ABC-SMC

The algorithm used corresponds to an Approximate Bayesian Computation approach using a Sequential Monte Carlo sampler. This iterative algorithm enhances the basic ABC algorithm by incorporating two main steps: weighted resampling of simulated particles and a gradual reduction in tolerance. Similar to the ABC rejection approach, a prior distribution is defined, aiming to estimate a posterior distribution. In ABC-SMC, this estimation is achieved sequentially by constructing intermediate distributions in each iteration, converging towards the posterior distribution.

The specific implementation used in the package improves upon Del Moral et al. (2006) original algorithm ¹ in three ways:

an adaptive threshold schedule selection based on quantiles of distances between simulated and observed data ² ³
an adaptive perturbation kernel width during the sampling step, dependent on the previous intermediate posterior distribution ⁴ ⁵
and the capability to use multiple criteria simultaneously.

Features

Parallelized ABC SMC: Conduct inference using an Approximate Bayesian Computation Sequential Monte Carlo algorithm, parallelized for enhanced computational efficiency.
Flexible Model Specification: Easily specify complex models.
Simulation Diagnostics: Optionally retain summary statistics and detailed model outputs in Parquet files, including long tabular trajectories.
Customizable Settings: Fine-tune algorithm parameters to suit specific modeling needs and computational resources.
Scalable: Utilize parallel computing capabilities to handle large datasets and complex models with ease.
Comprehensive Documentation: Detailed documentation and examples to guide users through package functionality and usage.

Installation

You can install the development version of BRREWABC from GitHub with:

# install.packages("devtools")
devtools::install_github("GaelBn/BRREWABC")

Usage

For a basic example, see the Get Started section.

Retaining summary statistics and model outputs

Model functions may return only a numeric distance vector, as in previous versions, or a structured result containing distances, summary statistics, and detailed outputs:

compute_dist <- function(x, ss_obs) {
  trajectory <- data.frame(
    timestep = 1:3,
    pop_id = "A",
    S = c(900, 880, 850),
    E = c(50, 55, 60),
    I = c(40, 50, 65),
    R = c(10, 15, 25)
  )

  list(
    distances = c(dist1 = sum((trajectory$I - ss_obs)^2)),
    summaries = list(
      epidemic_trajectory = trajectory,
      peak_infectious = max(trajectory$I)
    ),
    outputs = list(
      final_state = trajectory[nrow(trajectory), ]
    )
  )
}

Each named summary or output can be a vector, matrix, or data frame. For a given name, tabular column names and types must remain fixed across simulations, while the number of rows may vary.

Choose independently which summary statistics and outputs to retain:

result <- abcsmc(
  model_list = list(m1 = compute_dist),
  prior_dist = prior_dist,
  ss_obs = ss_obs,
  store_summaries = "all",
  store_outputs = "retained"
)

The available policies are "none", "retained", "accepted", and "all". Stored objects are buffered by each worker and written as bounded, atomic Parquet fragments under res/parquet. Buffer limits are independent of the simulation batch size and can be adjusted with storage_chunk_rows and storage_chunk_mb.

Storage and performance warning

Using store_summaries = "all" or store_outputs = "all" stores data for every tested particle, including rejected particles. This can require substantial disk space, especially when model outputs contain long trajectories or when the acceptance rate is low.

For large analyses, prefer "retained" or "accepted" unless rejected simulations are required for diagnostics. Ensure that sufficient local storage is available. Writing many fragments to synchronized or network storage may substantially reduce performance.

Parallel simulations are executed in bounded batches. Each local batch starts an isolated R process, so package loading and process startup are paid once per batch. When individual simulations are inexpensive, use a sufficiently large batch_size (for example 100 or more) to amortize this overhead. Smaller batches provide more responsive scheduling for expensive or highly variable simulations.

list_abc_stored_data(result)

accepted_trajectories <- read_summary_statistics(
  result,
  name = "epidemic_trajectory",
  generation = 5,
  status = "accepted"
)

rejected_trajectories <- read_summary_statistics(
  result,
  name = "epidemic_trajectory",
  status = "rejected"
)

retained_outputs <- read_model_outputs(
  result,
  name = "final_state",
  status = "retained"
)

For transfer or archival, fragments can be consolidated on demand into one Parquet file per named object and generation. The default export mode leaves the operational fragmented storage unchanged:

consolidate_abc_storage(
  result,
  kind = "summaries",
  generation = 5
)

Use dry_run = TRUE to inspect the planned files and volume. With mode = "replace", the consolidated files atomically replace the fragments in the active manifest; source fragments are retained unless keep_fragments = FALSE is requested.

The stored tables include attempt_id, generation, job_id, accepted, and retained, allowing each trajectory to be linked to the corresponding ABC attempt. See the storing-summaries-and-outputs vignette for further details.

Getting help

# Access package documentation
help(package = "BRREWABC")

Contributing

Contributions to BRREWABC are welcome! If you encounter any issues, have feature requests, or would like to contribute enhancements, please feel free to contact us.

Project status

This project is constantly evolving, according to needs and suggestions, at a pace that depends on the time that can be devoted to it.

Authors and acknowledgment

BRREWABC was developed by Gaël Beaunée thanks to the advice of a number of people, many thanks to them.

Troubleshooting

If you encounter any issues, please feel free to contact us.