```
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.0023 | -0.0058 |

[-0.3948, 0.3891] | [-0.0079, -0.0036] | |

PClass mean(2nd) - mean(1st) | -0.0032 | -0.19 |

[-0.3933, 0.3954] | [-0.27, -0.12] | |

PClass mean(3rd) - mean(1st) | 0.0029 | -0.38 |

[-0.3904, 0.4015] | [-0.45, -0.31] | |

SexCode mean(1) - mean(0) | 0.00093 | 0.49 |

[-0.40089, 0.41583] | [0.44, 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. :-4.008e-03 Min. :-0.821792 Min. :0.0000152
1st Qu.: 0.000e+00 1st Qu.:-0.352447 1st Qu.:0.0161285
Median : 0.000e+00 Median :-0.096142 Median :0.0940996
Mean : 2.443e-06 Mean :-0.211660 Mean :0.2115576
3rd Qu.: 0.000e+00 3rd Qu.:-0.016274 3rd Qu.:0.3669070
Max. : 3.929e-03 Max. :-0.000024 Max. :0.8087957
predicted_lo predicted_hi predicted tmp_idx
Min. :0.0000000 Min. :0.000e+00 Min. :0.000e+00 Min. : 1.0
1st Qu.:0.0000000 1st Qu.:0.000e+00 1st Qu.:0.000e+00 1st Qu.:189.8
Median :0.0000159 Median :2.253e-05 Median :1.709e-05 Median :378.5
Mean :0.0163892 Mean :1.591e-02 Mean :1.611e-02 Mean :378.5
3rd Qu.:0.0026629 3rd Qu.:3.071e-03 3rd Qu.:2.663e-03 3rd Qu.:567.2
Max. :0.2052452 Max. :2.052e-01 Max. :2.052e-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. :-1.409e-01 Min. :-0.218324 Min. :-0.065843 Min. :0.02483
1st Qu.:-1.227e-02 1st Qu.:-0.037970 1st Qu.:-0.008329 1st Qu.:0.21583
Median : 5.224e-04 Median :-0.008359 Median : 0.005345 Median :0.29603
Mean :-9.452e-05 Mean :-0.034518 Mean : 0.034909 Mean :0.33292
3rd Qu.: 2.625e-02 3rd Qu.: 0.010550 3rd Qu.: 0.091578 3rd Qu.:0.45950
Max. : 1.382e-01 Max. : 0.050567 Max. : 0.223201 Max. :0.89534
predicted_hi predicted tmp_idx PClass
Min. :0.02671 Min. :0.02515 Min. : 1.0 Length:4536
1st Qu.:0.22420 1st Qu.:0.21583 1st Qu.:189.8 Class :character
Median :0.31426 Median :0.29827 Median :378.5 Mode :character
Mean :0.33286 Mean :0.33289 Mean :378.5
3rd Qu.:0.42074 3rd Qu.:0.44007 3rd Qu.:567.2
Max. :0.90777 Max. :0.89534 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
```