Skip to content

Instantly share code, notes, and snippets.

@andrewheiss
Created November 17, 2022 02:16
Show Gist options
  • Save andrewheiss/a4e0c0ab2d735625ac17ec8a081f0f32 to your computer and use it in GitHub Desktop.
Save andrewheiss/a4e0c0ab2d735625ac17ec8a081f0f32 to your computer and use it in GitHub Desktop.
library(tidyverse)
library(tidybayes)
library(brms)
library(ggtext)

priors <- c(prior(normal(20, 5), class = Intercept),
            prior(normal(0, 2), class = b),
            prior(exponential(1), class = sigma),
            prior(exponential(1), class = sd),
            prior(lkj(2), class = cor),
            prior(beta(1, 1), class = hu))

priors %>% 
  parse_dist() %>% 
  # K = dimension of correlation matrix; 
  # ours is 2x2 here because we have one random slope
  marginalize_lkjcorr(K = 2) %>%
  mutate(nice_title = glue::glue("**{class}**<br>{prior}")) %>% 
  ggplot(aes(y = 0, dist = .dist, args = .args, fill = prior)) +
  stat_slab(normalize = "panels") +
  scale_fill_viridis_d(option = "plasma", end = 0.9) +
  facet_wrap(vars(nice_title), scales = "free_x") +
  guides(fill = "none") +
  labs(x = NULL, y = NULL) +
  theme_bw() +
  theme(strip.text = element_markdown(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())

# Run actual model
model <- brm(
  bf(outcome ~ year_centered + (1 + year_centered | country),
     decomp = "QR"),
  data = data,
  family = hurdle_lognormal(),
  prior = priors
)

Created on 2022-11-16 with reprex v2.0.2

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment