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()
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].
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.
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.
