library(brms)
library(ggplot2)
library(marginaleffects)
library(modelsummary)
options(brms.backend = "cmdstanr")
theme_set(theme_minimal())
<- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/Stat2Data/Titanic.csv")
titanic <- subset(titanic, PClass != "*")
titanic
<- Survived ~ SexCode + Age + PClass
f
<- brm(f,
mod_prior data = titanic,
prior = c(prior(normal(0, .2), class = b)),
cores = 4,
sample_prior = "only")
<- brm(f,
mod_posterior data = titanic,
cores = 4,
prior = c(prior(normal(0, .2), class = b)))
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:
Now, we use the avg_comparisons()
function from the marginaleffects
package to compute contrasts of interest:
<- list(
cmp "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 mean(+1) | 0.004 | -0.0059 |
[-0.378, 0.387] | [-0.0080, -0.0037] | |
PClass mean(2nd) - mean(1st) | 0.0016 | -0.19 |
[-0.3941, 0.3960] | [-0.27, -0.12] | |
PClass mean(3rd) - mean(1st) | 0.005 | -0.38 |
[-0.383, 0.402] | [-0.45, -0.31] | |
SexCode mean(1) - mean(0) | -0.0013 | 0.49 |
[-0.3910, 0.3884] | [0.43, 0.55] |
<- lapply(names(cmp), \(x) transform(posteriordraws(cmp[[x]]), Label = x))
draws <- do.call("rbind", draws)
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:
<- brm(
modcat_posterior ~ SexCode + Age,
PClass 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)
<- brm(
modcat_prior ~ SexCode + Age,
PClass 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)
<- posteriordraws(comparisons(modcat_prior))
pd
comparisons(modcat_prior) |> summary()
rowid term group contrast
Min. : 1.0 Length:4536 Length:4536 Length:4536
1st Qu.:189.8 Class :character Class :character Class :character
Median :378.5 Mode :character Mode :character Mode :character
Mean :378.5
3rd Qu.:567.2
Max. :756.0
estimate conf.low conf.high
Min. :-3.397e-03 Min. :-0.8089043 Min. :0.0000085
1st Qu.: 0.000e+00 1st Qu.:-0.3127726 1st Qu.:0.0174662
Median : 0.000e+00 Median :-0.1015558 Median :0.0994095
Mean : 6.530e-06 Mean :-0.2002571 Mean :0.2011653
3rd Qu.: 0.000e+00 3rd Qu.:-0.0168535 3rd Qu.:0.3291096
Max. : 4.388e-03 Max. :-0.0000266 Max. :0.8221026
predicted_lo predicted_hi predicted tmp_idx
Min. :0.000e+00 Min. :0.000e+00 Min. :0.000e+00 Min. : 1.0
1st Qu.:0.000e+00 1st Qu.:0.000e+00 1st Qu.:0.000e+00 1st Qu.:189.8
Median :1.865e-05 Median :2.029e-05 Median :1.865e-05 Median :378.5
Mean :1.582e-02 Mean :1.570e-02 Mean :1.569e-02 Mean :378.5
3rd Qu.:2.232e-03 3rd Qu.:2.604e-03 3rd Qu.:2.324e-03 3rd Qu.:567.2
Max. :1.983e-01 Max. :1.983e-01 Max. :1.983e-01 Max. :756.0
PClass SexCode Age
Length:4536 Min. :0.000 Min. : 0.17
Class :character 1st Qu.:0.000 1st Qu.:21.00
Mode :character Median :0.000 Median :28.00
Mean :0.381 Mean :30.40
3rd Qu.:1.000 3rd Qu.:39.00
Max. :1.000 Max. :71.00
comparisons(modcat_posterior) |> summary()
rowid term group contrast
Min. : 1.0 Length:4536 Length:4536 Length:4536
1st Qu.:189.8 Class :character Class :character Class :character
Median :378.5 Mode :character Mode :character Mode :character
Mean :378.5
3rd Qu.:567.2
Max. :756.0
estimate conf.low conf.high predicted_lo
Min. :-0.1400135 Min. :-0.217763 Min. :-0.067315 Min. :0.02483
1st Qu.:-0.0122534 1st Qu.:-0.043127 1st Qu.:-0.008258 1st Qu.:0.21587
Median : 0.0005224 Median :-0.008488 Median : 0.005366 Median :0.29646
Mean :-0.0001023 Mean :-0.035289 Mean : 0.035449 Mean :0.33291
3rd Qu.: 0.0260982 3rd Qu.: 0.010511 3rd Qu.: 0.093869 3rd Qu.:0.45905
Max. : 0.1379084 Max. : 0.049241 Max. : 0.223777 Max. :0.89539
predicted_hi predicted tmp_idx PClass
Min. :0.02703 Min. :0.02515 Min. : 1.0 Length:4536
1st Qu.:0.22474 1st Qu.:0.21587 1st Qu.:189.8 Class :character
Median :0.31442 Median :0.29885 Median :378.5 Mode :character
Mean :0.33293 Mean :0.33291 Mean :378.5
3rd Qu.:0.41974 3rd Qu.:0.44152 3rd Qu.:567.2
Max. :0.90783 Max. :0.89539 Max. :756.0
SexCode Age
Min. :0.000 Min. : 0.17
1st Qu.:0.000 1st Qu.:21.00
Median :0.000 Median :28.00
Mean :0.381 Mean :30.40
3rd Qu.:1.000 3rd Qu.:39.00
Max. :1.000 Max. :71.00