Posthoc contrasts with emmeans, tidybayes, and brms

bayesian_posthoc.R

library(tidyverse)
library(emmeans)
library(brms)
library(tidybayes)

warp.brms <- brm(breaks ~ wool * tension, data = warpbreaks)

#get the adjusted means
warp_em <- emmeans (warp.brms,  ~ tension | wool)
warp_em

#get all possible contrasts
cont <- contrast(warm_em, "tukey")
cont

#get the posterior draws from the contrasts
cont_posterior <- gather_emmeans_draws(cont)

#plot
ggplot(cont_posterior,
       aes(y = contrast, x = .value)) +
  geom_halfeyeh() +
  facet_wrap(~wool) +
  geom_vline(xintercept = 0, color = "red", lty = 2)


ggplot(cont_posterior,
       aes(y = contrast, x = .value, fill = wool, group = wool)) +
  geom_halfeyeh(alpha = 0.5)+
  geom_vline(xintercept = 0, color = "red", lty = 2)

#need to figure out how to make this look all posterior-y
emmip(warp.brms,  ~ tension | wool, CIs = TRUE, cov.reduce = range)

This is great!

Here's some quick brainstorming inspired by your the comment on your call to emmip...

warp.brms %>%
  emmeans( ~ tension | wool) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = tension, y = .value)) +
  stat_lineribbon() +
  scale_fill_brewer(palette = "Greys") +
  facet_grid(~ wool) +
  theme_light()

warp.brms %>%
  emmeans( ~ tension | wool) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = tension, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  facet_grid(~ wool) +
  theme_light()

(reminds me I need to implement the vertical geom_halfeye sometime... or just a generalization of all of the interval+density variants)

warp.brms %>%
  emmeans( ~ tension | wool) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = tension, y = .value, group = .draw)) +
  geom_line(alpha = .01) +
  scale_fill_brewer() +
  facet_grid(~ wool) +
  theme_light()

warp.brms %>%
  emmeans( ~ tension | wool) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = tension, y = .value, fill = wool, color = wool)) +
  stat_lineribbon(alpha = 1/4) +
  theme_light()

Too many intervals here probably, it's a little mind-bending...

warp.brms %>%
  emmeans( ~ tension | wool) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = tension, y = .value, fill = wool, color = wool)) +
  stat_lineribbon(alpha = 1/3, .width = .95) +
  theme_light()

Or you could go the other way and approximate a gradient plot (mind you this brings all the perceptual issues with gradients...)

warp.brms %>%
  emmeans( ~ tension | wool) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = tension, y = .value, fill = wool, color = wool)) +
  stat_lineribbon(alpha = 1/41, .width = ppoints(40)) +
  theme_light()

Heh, the way stat_lineribbon interleaves layers makes that last one kind of painterly. Neat, but I wouldn't use it in practice :)

warp.brms %>%
  emmeans( ~ tension | wool) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = tension, y = .value)) +
  stat_lineribbon(fill = "gray25", color = "gray25", size = .5, alpha = 1/41, .width = ppoints(40)) +
  stat_pointinterval() +
  facet_grid(~ wool) +
  theme_light()

warp.brms %>%
  emmeans( ~ tension | wool) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = tension, y = .value)) +
  geom_line(aes(group = .draw), alpha = .005) +
  stat_pointinterval() +
  facet_grid(~ wool) +
  theme_light()

Hey thank you so much for this code! This definitely helped me to get going myself. I am wondering though is there any reason why you used contrast(..., "tukey") instead of pairwise? I am getting the same results for both for my model but I want to be sure that I am not missing anything about using tukey for Bayesian mdoels.