Vincent Arel-Bundock - Prior Predictive Checks with marginaleffects and brms

Bayesians often advocate for the use of prior predictive checks (Gelman et al. 2020). The idea is to simulate from the model, without using the data, in order to refine the model before fitting. For example, we could draw parameter values from the priors, and use the model to simulate values of the outcome. Then, could inspect those to determine if the simulated outcomes (and thus the priors) make sense substantively. Prior predictive checks allow us to iterate on the model without looking at the data multiple times.

One major challenge lies in interpretation: When the parameters of a model are hard to interpret, the analyst will often need to transform before they can assess if the generated quantities make sense, and if the priors are an appropriate representation of available information.

In this post I show how to use the marginaleffects and brms packages for R to facilitate this process. The benefit of the approach described below is that it allows us to conduct prior predictive checks on the actual quantities of interest. For example, if the ultimate quantity that we want to estimate is a contrast or an Average Treatment Effect, then we can use marginaleffects to simulate the specific quantity of interest using just the priors and the model.

In this example, we create two model objects with brms. In one of them, we set sample_prior="only" to indicate that we do not want to use the dataset at all, and that we only want to use the priors and model for simulation:

library(brms)
library(ggplot2)
library(marginaleffects)
library(modelsummary)
options(brms.backend = "cmdstanr")
theme_set(theme_minimal())

titanic <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/Stat2Data/Titanic.csv")
titanic <- subset(titanic, PClass != "*")

f <- Survived ~ SexCode + Age + PClass

mod_prior <- brm(f,
    data = titanic,
    prior = c(prior(normal(0, .2), class = b)),
    cores = 4,
    sample_prior = "only")

mod_posterior <- brm(f,
    data = titanic,
    cores = 4,
    prior = c(prior(normal(0, .2), class = b)))

Now, we use the avg_comparisons() function from the marginaleffects package to compute contrasts of interest:

cmp <- list(
    "Prior" = avg_comparisons(mod_prior),
    "Posterior" = avg_comparisons(mod_posterior))

Finally, we compare the results with and without the data in tables and plots:

modelsummary(
    cmp,
    output = "markdown",
    statistic = "conf.int",
    fmt = fmt_significant(2),
    gof_map = NA,
    shape = term : contrast ~ model)

	Prior	Posterior
Age +1	0.0017	-0.0058
	[-0.3884, 0.4063]	[-0.0080, -0.0037]
PClass 2nd - 1st	0.00048	-0.19
	[-0.38974, 0.38988]	[-0.26, -0.12]
PClass 3rd - 1st	-0.0023	-0.38
	[-0.3957, 0.4042]	[-0.45, -0.30]
SexCode 1 - 0	0.0023	0.49
	[-0.3967, 0.3964]	[0.43, 0.55]

draws <- lapply(names(cmp), \(x) transform(posteriordraws(cmp[[x]]), Label = x))
draws <- do.call("rbind", draws)

ggplot(draws, aes(x = draw, color = Label)) +
    xlim(c(-1, 1)) +
    geom_density() +
    facet_wrap(~term + contrast, scales = "free")

This kind of approach is particularly useful with more complicated models, such as this one with categorical outcomes. In such models, it would be hard to know if a normal prior is appropriate for the different parameters:

modcat_posterior <- brm(
    PClass ~ SexCode + Age,
    prior = c(
        prior(normal(0, 3), class = b, dpar = "mu2nd"),
        prior(normal(0, 3), class = b, dpar = "mu3rd")),
    family = categorical(link = logit),
    cores = 4,
    data = titanic)

modcat_prior <- brm(
    PClass ~ SexCode + Age,
    prior = c(
        prior(normal(0, 3), class = b, dpar = "mu2nd"),
        prior(normal(0, 3), class = b, dpar = "mu3rd")),
    family = categorical(link = logit),
    sample_prior = "only",
    cores = 4,
    data = titanic)

pd <- posteriordraws(comparisons(modcat_prior))

comparisons(modcat_prior) |> summary()


 Group    Term          Contrast  Estimate   2.5 % 97.5 %
   1st SexCode mean(1) - mean(0) -1.79e-05 -0.1793 0.1963
   1st Age     mean(+1)           1.38e-03 -0.0338 0.0302
   2nd SexCode mean(1) - mean(0) -3.34e-04 -0.2155 0.2054
   2nd Age     mean(+1)           1.58e-03 -0.0349 0.0343
   3rd SexCode mean(1) - mean(0)  1.73e-06 -0.2297 0.2352
   3rd Age     mean(+1)           1.82e-03 -0.0352 0.0333

Columns: group, term, contrast, estimate, conf.low, conf.high

comparisons(modcat_posterior) |> summary()


 Group    Term          Contrast Estimate    2.5 %    97.5 %
   1st SexCode mean(1) - mean(0)   0.0995  0.03926  1.60e-01
   1st Age     mean(+1)            0.0128  0.01110  1.43e-02
   2nd SexCode mean(1) - mean(0)   0.0231 -0.04077  8.70e-02
   2nd Age     mean(+1)           -0.0020 -0.00417  7.98e-05
   3rd SexCode mean(1) - mean(0)  -0.1226 -0.18895 -5.45e-02
   3rd Age     mean(+1)           -0.0107 -0.01283 -8.51e-03

Columns: group, term, contrast, estimate, conf.low, conf.high

References

Gelman, Andrew, Aki Vehtari, Daniel Simpson, Charles C. Margossian, Bob Carpenter, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner, and Martin Modrák. 2020. “Bayesian Workflow.” https://arxiv.org/abs/2011.01808.