Background

The epidist package uses Bayesian inference to estimate delay distributions and other quantities. Doing Bayesian inference amounts to approximating the posterior distribution of each parameter in the statistical model¹. A range of methods exist to perform this approximation.

By default, epidist uses the No-U-Turn Sampler [NUTS; Hoffman et al. (2014)] Hamiltonian Monte Carlo (HMC) algorithm. NUTS is an example of a broader class of Markov chain Monte Carlo (MCMC) methods. These methods work by simulating from a Markov chain with the target posterior distribution as its stationary distribution. When MCMC algorithms are run for sufficiently many iterations, and have reached convergence, the samples can be treated as being drawn from the posterior distribution. Relevant posterior quantities such as expectations may then be computed using these samples. A drawback of MCMC methods, like NUTS, is that simulations can be quite computational intensive, especially for complex models or large data.

The epidist package is built using brms (Bürkner 2017), which stands for “Bayesian Regression Models using Stan”, where Stan (Carpenter et al. 2017) is a probabilistic programming language. Although NUTS is the primary inference algorithm used in Stan, additional options are available. These additional inference algorithms are also available in epidist due to its dependence on brms.

In this vignette, we first briefly describe the alternative algorithms available (Section @ref(other)) as well as directing you to more detailed resources. Then (Section @ref(demo)) we demonstrate their application to fitting simulated data, before extracting and comparing posterior distributions. By comparing the resulting inferences to those from NUTS, we hope to help you make informed decisions about which algorithm to use in your applied problem.

Alternative approximate inference algorithms

Here we describe three alternative approximate Bayesian inference algorithms that are available to use in brms, and therefore also available in epidist. It’s worth noting that further inference algorithms may have become available since this vignette was last updated. Please check brms package updates if interested!

Demonstration

In this demonstration, we use the following packages:

library(epidist)
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
library(tibble)
library(tidybayes)
library(cmdstanr) # nolint

To access the approximate inference methods used in this vignette we will need to use the cmdstanr backend for brms (we generally recommend using this backend for fitting models). To do this, we first need to install CmdStan (see the README for more details). We can check we have everything we need as follows:

cmdstanr::cmdstan_version()
#> [1] "2.40.0"

We can simulate data to use for fitting models. The example data simulation process follows that used in the Getting started with epidist vignette, so we will not detail exactly what is happening here, but please consult that vignette if interested:

meanlog <- 1.8
sdlog <- 0.5
obs_time <- 25
sample_size <- 200

obs_cens_trunc_samp <- simulate_gillespie(seed = 101) |>
  simulate_secondary(
    meanlog = meanlog,
    sdlog = sdlog
  ) |>
  mutate(
    ptime_lwr = floor(.data$ptime),
    ptime_upr = .data$ptime_lwr + 1,
    stime_lwr = floor(.data$stime),
    stime_upr = .data$stime_lwr + 1,
    obs_time = obs_time
  ) |>
  filter(.data$stime_upr <= .data$obs_time) |>
  slice_sample(n = sample_size, replace = FALSE)

We now prepare the data for fitting with the marginal model. We first pre-compile the Stan model so that compilation time is excluded from the timing comparisons below.

linelist_data <- as_epidist_linelist_data(
  obs_cens_trunc_samp$ptime_lwr,
  obs_cens_trunc_samp$ptime_upr,
  obs_cens_trunc_samp$stime_lwr,
  obs_cens_trunc_samp$stime_upr,
  obs_time = obs_cens_trunc_samp$obs_time
)

data <- as_epidist_marginal_model(linelist_data)

# Pre-compile the model so compilation time is not included
# in the timing comparisons
fit_compile <- epidist(
  data = data, backend = "cmdstanr",
  chains = 1, iter = 5
)

We now perform inference with HMC:

t <- proc.time()
fit_hmc <- update(fit_compile, chains = 4, iter = 2000)
time_hmc <- proc.time() - t

Note that for clarity above we specify algorithm = "sampling", but if you were to call epidist(data = data) the result would be the same since "sampling" (i.e. HMC) is the default value for the algorithm argument.

Now, we fit² the same marginal model using each method in Section @ref(other). To match the four Markov chains of length 1000 in HMC above, we then draw 4000 samples from each approximate posterior.

t <- proc.time()
fit_laplace <- update(
  fit_compile,
  algorithm = "laplace", draws = 4000
)
time_laplace <- proc.time() - t

t <- proc.time()
fit_advi <- update(
  fit_compile,
  algorithm = "meanfield", draws = 4000
)
time_advi <- proc.time() - t

For the Pathfinder algorithm we will set num_paths = 1.

t <- proc.time()
fit_pathfinder <- update(
  fit_compile,
  algorithm = "pathfinder", draws = 4000,
  chains = 1
)
time_pathfinder <- proc.time() - t

We now extract the posterior distribution for the delay parameters from the fitted model for each inference method. Thankfully, each algorithm is implemented to sample draws from the posterior distribution, making post-processing straightforward.

fits <- list(
  HMC = fit_hmc,
  Laplace = fit_laplace,
  ADVI = fit_advi,
  Pathfinder = fit_pathfinder
)

draws <- imap(fits, function(fit, name) {
  draws <- fit |>
    epidist_strata() |>
    add_delay_parameter_draws(fit) |>
    add_summaries() |>
    as.data.frame() |>
    pivot_longer(
      cols = c("mu", "sigma", "mean", "sd"),
      names_to = "parameter",
      values_to = "value"
    ) |>
    filter(parameter %in% c("mu", "sigma")) |>
    mutate(method = as.factor(name))
  return(draws)
})

draws <- bind_rows(draws)

Comparison of parameter posterior distributions

The mean estimated value of each parameter, from each method, is as follows.

pars <- draws |>
  group_by(method, parameter) |>
  summarise(value = mean(value)) |>
  pivot_wider(names_from = parameter, values_from = value) |>
  ungroup()

pars
#> # A tibble: 4 × 3
#>   method        mu sigma
#>   <fct>      <dbl> <dbl>
#> 1 HMC         1.75 0.460
#> 2 Laplace     1.74 0.455
#> 3 ADVI        1.77 0.420
#> 4 Pathfinder  1.74 0.456

More comprehensively, the estimated posterior distributions are shown in Figure @ref(fig:posterior).

Click to expand for code to create posterior distribution plot

p_posterior <- draws |>
  ggplot(aes(x = value, col = method)) +
  stat_slabinterval(density = "histogram", breaks = 30, alpha = 0.8) +
  scale_colour_manual(values = c("#56B4E9", "#009E73", "#E69F00", "#CC79A7")) +
  facet_grid(method ~ parameter, scales = "free_x") +
  theme_minimal() +
  guides(fill = "none") +
  labs(x = "", y = "", col = "Method") +
  theme(legend.position = "bottom")

(ref:posterior) Estimated posterior distributions for the mu and sigma parameters using each inference method, shown using tidybayes::stat_slabinterval().

p_posterior

(#fig:posterior)(ref:posterior)

Comparison of resulting delay distributions

Figure @ref(fig:delay-pdf) shows how the different mu and sigma posterior mean estimates from each inference method alter an estimated delay distribution.

Click to expand for code to create delay PDF plot

p_delay_pdf <- pmap_df(
  filter(pars), ~ tibble(
    x = seq(0, 25, by = 0.1),
    method = ..1, density = dlnorm(x, ..2, ..3)
  )
) |>
  ggplot(aes(x = x, y = density, col = method)) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#56B4E9", "#009E73", "#E69F00", "#CC79A7")) +
  theme_minimal() +
  labs(x = "", y = "", col = "Method") +
  theme(legend.position = "bottom")

(ref:delay-pdf) Delay probability density functions obtained based on the posterior mean estimated mu and sigma parameters.

p_delay_pdf

(#fig:delay-pdf)(ref:delay-pdf)

times <- list(
  HMC = time_hmc,
  Laplace = time_laplace,
  ADVI = time_advi,
  Pathfinder = time_pathfinder
)

times
#> $HMC
#>    user  system elapsed
#>  12.551   0.206  12.771
#>
#> $Laplace
#>    user  system elapsed
#>   1.296   0.110   1.502
#>
#> $ADVI
#>    user  system elapsed
#>   1.029   0.082   1.122
#>
#> $Pathfinder
#>    user  system elapsed
#>   0.396   0.083   0.528

Conclusion

The range of alternative approximation algorithms available, and their ease of use, is an attractive feature of brms. We found that these algorithms do produce reasonable approximations in far less time than HMC. Of course, this vignette only includes one example, so results may differ for other models and datasets.