# setup
knitr::opts_chunk$set(
  fig.asp = 0.618,
  fig.width = 6,
  dpi = 300,
  fig.align = "center",
  # out.width = "80%",
  comment = "#>",
  collapse = TRUE)

Code that I am using to make plots for this twitter thread.


Fit a model and get some posterior samples

library(lme4)
#> Loading required package: Matrix
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(tibble)
library(ggplot2)

sleepstudy <- sleepstudy %>%
  as_tibble() %>%
  mutate(Subject = as.character(Subject))

# Add two fake participants to show wide parameter densities
df_sleep <- sleepstudy %>%
  add_row(Reaction = 245, Days = 0, Subject = "373") %>%
  add_row(Reaction = c(286, 288), Days = 0:1, Subject = "374")

# Fit mixed effect models
m <- lmer(Reaction ~ 1 + Days + (1 + Days | Subject), df_sleep)

# Take "informal" posterior samples
sims <- arm::sim(m, 1000)

# Convert to long format
ranef_samples <- ranef(sims) %>%
  reshape2::melt(c("draw", "Subject", "Parameter")) %>%
  as_tibble()

# Sort the participants by median intercept value
subject_order <- ranef_samples %>%
  group_by(Subject, Parameter) %>%
  summarise(median = median(value)) %>%
  filter(Parameter == "(Intercept)") %>%
  ungroup() %>%
  mutate(PlotSubject = factor(Subject) %>% forcats::fct_reorder(median)) %>%
  select(Subject, PlotSubject)

ranef_samples <- left_join(ranef_samples, subject_order)
#> Joining, by = "Subject"

Draw the familiar caterpillar plot of samples

stat_interval <- function(..., conf.int = .95) {
  stat_summary(..., fun.data = median_hilow,
               geom = "linerange", fun.args = list(conf.int = conf.int))
}

stat_median <- function(...) {
  stat_summary(..., fun.y = median, geom = "point")
}

ggplot(ranef_samples) +
  aes(x = PlotSubject, y = value) +
  stat_interval(conf.int = .95, size = 1) +
  stat_interval(conf.int = .80, size = 2) +
  stat_median(size = 3.5) +
  coord_flip() +
  facet_grid(~ Parameter, scales = "free_x") +
  labs(x = "Subject", y = NULL) + 
  labs(caption = "Medians with 95% (thin), 80% (thick) intervals")

Try geom_joy() to draw density curves instead of intervals

The default height variable uses the magic ..density variable which I have never fully understood in ggplot2.

# devtools::install_github("clauswilke/ggjoy")
library(ggjoy)

ggplot(ranef_samples) +
  aes(x = value, y = PlotSubject, height = ..density..,
      group = interaction(PlotSubject, Parameter)) +
  scale_y_discrete(expand=c(0.01, 0)) +
  scale_x_continuous(expand=c(0, 0)) +
  facet_grid(~ Parameter) +
  labs(y = "Subject", x = NULL) +
  geom_joy2()