Working with posterior samples

The output of a call to epidist is compatible with typical Stan workflows. We recommend use of the posterior package for working with samples.

library(posterior)
draws <- as_draws_df(fit, variable = c("Intercept", "Intercept_sigma"))
head(draws)

## # A draws_df: 6 iterations, 1 chains, and 2 variables
##   Intercept Intercept_sigma
## 1       1.9           -0.43
## 2       1.7           -0.64
## 3       1.7           -0.64
## 4       1.7           -0.59
## 5       1.8           -0.60
## 6       1.8           -0.43
## # ... hidden reserved variables {'.chain', '.iteration', '.draw'}

library(posterior)
# Convert model parameters to rvars
rv <- as_draws_rvars(fit_location)

# Look at the intercept in rvar form
rv$b_Intercept

## rvar<500,2>[1] mean ± sd:
## [1] 1.6 ± 0.066

# Calculate the mean delay (on natural scale) for each location
mean_delay_loc0 <- exp(rv$b_Intercept + 0.5 * exp(rv$b_sigma_Intercept)^2)
mean_delay_loc1 <- exp(rv$b_Intercept +
  rv$b_location + 0.5 * (exp(rv$b_sigma_Intercept + rv$b_sigma_location))^2)

# Summarize the distributions
summary(mean_delay_loc0)

## # A tibble: 1 × 10
##   variable         mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>           <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mean_delay_loc0  5.53   5.47 0.453 0.400  4.91  6.36  1.00     841.     727.

summary(mean_delay_loc1)

## # A tibble: 1 × 10
##   variable         mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>           <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mean_delay_loc1  14.1   11.5  11.6  3.14  8.20  25.4  1.02     179.     137.

# Calculate the difference between locations
delay_diff <- mean_delay_loc1 - mean_delay_loc0
delay_diff

## rvar<500,2>[1] mean ± sd:
## [1] 8.6 ± 12

How can I assess whether sampling has converged?

The output of a call to epidist is compatible with typical Stan workflows. We recommend use of the bayesplot package for sampling diagnostic plots. For example, the function bayesplot::mcmc_trace() can be used to produce traceplots for specified parameters.

library(bayesplot)
mcmc_trace(fit, pars = c("Intercept", "Intercept_sigma"))

(#fig:unnamed-chunk-4)plot of chunk unnamed-chunk-4

We also provide a function epidist_diagnostics() which can be used to obtain common diagnostics used to assess the quality of a fitted model.

epidist_diagnostics(fit)

## # A tibble: 1 × 8
##    time samples max_rhat divergent_transitions per_divergent_transitions
##   <dbl>   <dbl>    <dbl>                 <dbl>                     <dbl>
## 1  3.97    1000     1.01                     0                         0
## # ℹ 3 more variables: max_treedepth <dbl>, no_at_max_treedepth <int>,
## #   per_at_max_treedepth <dbl>

How can I perform posterior predictive checks?

Posterior predictive checks are a useful tool for assessing model fit by comparing observed data to simulated data from the posterior predictive distribution. brms::pp_check() is a convenient function for this, but it does not automatically handle the aggregated/weighted data structure used by epidist models. Because brms::pp_check() ignores the weights argument when plotting observed data, passing aggregated data directly will result in an incorrect comparison between observed and predicted distributions.

To use pp_check() correctly with aggregated data, you must first expand the data to an individual-level format (one row per observation) using tidyr::uncount(). Then, pass this expanded data to the newdata argument of pp_check().

# Expand the aggregated data to individual-level data
# Note: We must ensure the weight column 'n' is present but set to 1
data_expanded <- data |>
  tidyr::uncount(weights = n) |>
  mutate(n = 1)

## Warning: ! Dropping the <epidist_marginal_model> class because the object no longer
##   meets its requirements:
## ✖ Assertion on 'names(data)' failed: Names must include the elements
##   {'delay_lwr','delay_upr','delay_min','relative_obs_time','pwindow','swindow','n'},
##   but is missing elements {'n'}.
## ℹ Use the matching `as_epidist_*()` function to recreate the object.

# Run pp_check with the expanded data
pp_check(fit, newdata = data_expanded, ndraws = 100)

(#fig:unnamed-chunk-6)plot of chunk unnamed-chunk-6

For more advanced custom visualizations, you can also use tidybayes::add_predicted_draws() as demonstrated in the “Advanced features with Ebola data” vignette.

What is an observation in the meta model?

as_epidist_meta_model() fits the summaries a study computed from the same delays jointly, so one row of the model data is a group of summaries rather than a single reported value. This is the same situation as aggregated data, where one row stands for many individual delays. brms::log_lik() and loo::loo() work at the row level, so their values are per group, and pp_check() needs the same care as described above. Predictions for summary rows are the reported summary itself rather than a delay, so restrict pp_check() to the individual level rows of a mixed fit.

I’d like to run a simulation study

We recommend use of the purrr package for running many epidist models, for example as a part of a simulation study. We particularly highlight two functions which might be useful:

1. purrr::map() (and other similar functions) for iterating over a list of inputs.
2. purrr::safely() which ensures that the function called “always succeeds”. In other words, if there is an error it will be captured and output, rather than ending computation (and potentially disrupting a call to purrr::map()).

For an example use of these functions, have a look at the epidist-paper repository containing the code for Park et al. (2024). (Note that in that codebase, we use map as a part of a targets pipeline.)

How are the default priors for epidist chosen?

brms provides default priors for all parameters. However, some of those priors do not make sense in the context of our application. Instead, we used prior predictive checking to set epidist-specific default priors which produce epidemiological delay distribution mean and standard deviation parameters in a reasonable range.

For example, for the brms::lognormal() latent individual model, we suggest the following prior distributions for the brms mu and sigma intercept parameters:

# Note that we export lognormal() as part of epidist hence no need for brms::
family <- lognormal()

epidist_family <- epidist_family(data, family)
epidist_formula <- epidist_formula(
  data,
  family = epidist_family,
  formula = mu ~ 1
)

# NULL here means no replacing priors from the user!
epidist_prior <- epidist_prior(
  data = data,
  family = family,
  formula = epidist_formula,
  prior = NULL
)

epidist_prior

##                 prior     class coef group resp  dpar nlpar lb ub tag  source
##    student_t(3, 5, 3) Intercept                                       default
##  student_t(3, 0, 2.5) Intercept                 sigma                 default

(Note that the functions epidist_family() and epidist_prior() are mostly for internal use!)

Here are the distributions on the delay distribution mean and standard deviation parameters that these prior distributions imply:

set.seed(1)
fit_ppc <- epidist(
  data = data,
  formula = mu ~ 1,
  family = lognormal(),
  sample_prior = "only",
  seed = 1,
  backend = "cmdstanr"
)

pred <- fit_ppc |>
  epidist_strata() |>
  add_delay_parameter_draws(fit_ppc) |>
  add_summaries()

plot(pred, pars = c("mean", "sd")) +
  scale_x_log10(labels = comma)

(#fig:unnamed-chunk-9)plot of chunk unnamed-chunk-9

quantile(pred$mean, c(0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99))

##         1%        10%        25%        50%        75%        90%        99%
##  0.3172760  0.8667011  1.6163277  3.2461963  6.6291746 12.0852973 34.0415723

quantile(pred$sd, c(0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99))

##         1%        10%        25%        50%        75%        90%        99%
##  0.1335500  0.3950604  0.7687361  1.7107506  3.9206460  8.1670947 33.4423697

How can I assess how sensitive the fitted posterior distribution is to the prior distribution used?

We recommend use of the priorsense package (Kallioinen et al. 2024) to check how sensitive the posterior distribution is to perturbations of the prior distribution and likelihood using power-scaling analysis:

library(priorsense)
powerscale_plot_dens(fit, variable = c("Intercept", "Intercept_sigma")) +
  theme_minimal()

(#fig:unnamed-chunk-10)plot of chunk unnamed-chunk-10

What do the parameters in my model output correspond to?

The epidist package uses brms to fit models. This means that the model output will include brms-style names for parameters. Here, we provide a table giving the correspondence between the distributional parameter names used in brms and those used in standard R functions for some common likelihood families.

Note that all families in brms are parameterised with some measure of centrality mu as their first parameter. This parameter does not necessarily correspond to the mean, hence the provision of a function add_summaries() within epidist to add columns containing the natural scale mean, standard deviation and quantiles to a data.frame of draws.

How can I generate predictions with my fitted epidist model?

It is possible to generate predictions manually by working with samples from the model output. However this is tricky to do, and so where possible we recommend using the tidybayes package. In particular, following functions may be useful:

1. tidybayes::add_epred_draws() for predictions of the expected value of a delay.
2. tidybayes::add_linpred_draws() for predictions of the delay distributional parameter linear predictors.
3. tidybayes::add_predicted_draws() for predictions of the observed delay.

To see these functions demonstrated in a vignette, see “Advanced features with Ebola data”. As a short example, to generate 4000 predictions (equal to the number of draws) of the delay that would be observed with a double censored observation process (in which the primary and secondary censoring windows are both one) then:

library(tidybayes)
draws_pmf <- tibble::tibble(
  relative_obs_time = Inf, pwindow = 1, swindow = 1, delay_upr = NA,
  delay_min = 0
) |>
  add_predicted_draws(fit)

## Warning: Found infinite values in the data, which may cause issues for Stan.

ggplot(draws_pmf, aes(x = .prediction)) +
  geom_bar(aes(y = after_stat(count / sum(count)))) +
  labs(x = "Delay", y = "PMF") +
  scale_x_continuous(limits = c(0, 30)) +
  theme_minimal()

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_bar()`).

(#fig:unnamed-chunk-11)plot of chunk unnamed-chunk-11

Importantly, this functionality is only available for epidist models using brms families that have a log_lik_censor method implemented internally in brms. If you are using another family, consider submitting a pull request to implement these methods!

library(marginaleffects)

# Compare outcomes between location categories (0 and 1)
avg_comparisons(
  fit_location,
  variables = list(location = c(0, 1))
)

## Warning: The training data could not be extracted reliably from the model object. Column types, factor levels, or values may have been coerced or altered. It is safer to use `set_modeldata(model, data)` to attach the training data explicitly. Original warning from `insight::get_data()`: Could not recover model data from environment. Please make sure your
##   data is available in your workspace.
##   Trying to retrieve data from the model frame now. Original warning from `insight::get_data()`: Following potential variables could not be found in the data: n) +
##   vreal(relative_obs_time, pwindow, swindow, delay_upr, delay_min This warning appears once per session.

##
##  Estimate 2.5 % 97.5 %
##      6.02  2.12   29.8
##
## Term: location
## Type: response
## Comparison: 1 - 0

cmdstanr::cmdstan_version()

## [1] "2.40.0"