# 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.

Best,

George

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.

``````library(tidyverse)
library(brms)
library(tidybayes)
library(lme4)
library(lattice)

# 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"))+
theme_linedraw()
``````

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.