Visualize Bayesian model

Hello there,

I’m trying to visualize a Bayesian model in R using likelihood, prior, and posterior (see picture below):


My R code is the following:

prior <- set_prior("student_t(3,0,2.5)", class = "b")
fit <- brm(F1 ~  0  + Intercept + vowel +  group + vowel:group + (1 | participant) , data = test, prior = prior, sample_prior = "yes")

I’m struggling with the code, so any ideas are welcome.



The prior and posterior are fairly straightforward to extract from your fitted model using the prior_draws and as_draws_df function. I’ve made an attempt at the “likelihood,” but I’m less confident that I’ve done this correctly. Briefly, my approach was to re-fit the same model using lme4::lmer, then use profile.merMod to calculate the zeta function, and finally fit a Gaussian kernel density estimator on zeta to estimate the probability distribution. I’m hoping someone who knows this better can weigh in and/or identify a more elegant solution.


# get data
data(sleepstudy, package = "lme4")

# fit example model
prior1 <- prior(normal(0, 3), class = b)

fit_mcmc <- brm(Reaction ~ Days + (1 | Subject),
                data = sleepstudy, prior = prior1,
                cores = 2, backend = "cmdstanr",
                sample_prior = "yes")

fit_reml <- lmer(Reaction ~ Days + (1 | Subject),
                 data = sleepstudy)

# extract prior and posterior samples for Days
prior_draws <- tibble(b_Days = prior_draws(fit_mcmc, variable = "b")$b)
post_draws <- as_tibble(as_draws_df(fit_mcmc))

draws <- bind_rows(prior = prior_draws, posterior = post_draws, .id = "dist")

# extract parameter density
prof <- densityplot(profile(fit_reml, "Days"), bw = "SJ", kernel = "gaussian")
lik <- tibble(b_Days = seq(prof$x.limits[[1]][1], prof$x.limits[[1]][2], length.out = 201),
              density = prof$call$data$density,
              dist = "likelihood")

# plot
ggplot(draws, aes(x = b_Days, fill = dist))+
  stat_slab(alpha = 0.5)+
  geom_ribbon(data = lik, aes(ymax = density), ymin = 0, color = "transparent", alpha = 0.5)+
  scale_fill_manual(name = "", values = c("#4477AA", "#228833", "#CCBB44"))+

1 Like

There’s a normalization conflict between the distributions that come from samples and that from the likelihood kernel density. stat_slab offers some normalization options. normalize = "none" looks the most correct to me despite the warning in the documentation not to use it except for values in [0,1].

1 Like

For examples using much simpler models (e.g., teaching), see also examples from @Solomon’s recoding of Kruschke’s Doing Bayesian Data Analysis:
5 Bayes’ Rule | Doing Bayesian Data Analysis in brms and the tidyverse.

1 Like

Adding in the likelihood will be difficult for anything other than the simplest of models. I recommend focusing on plotting the prior and/or the posterior.