Joyplot of random effect distributions
# 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()
The problem here is that the x-axis should vary between the two plots. Fix that.
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)) +
labs(y = "Subject", x = NULL) +
geom_joy2() +
facet_grid(~ Parameter, scales = "free_x") 
Now the problem is that the height is scaled for both the facets together instead of the calculating the density separately within each facet.
Don’t use ..density..
I am going to create a function to compute the density/height values manually.
# adapted from ggplot2:::compute_density()
quick_density_df <- function(x, interval_width = 1) {
nx <- length(x)
tail_width <- (1 - interval_width) / 2
qs <- quantile(x, probs = c(tail_width, 1 - tail_width))
dens <- density(
x = x,
bw = "nrd0",
adjust = 1,
kernel = "gaussian",
n = 1024,
from = min(qs),
to = max(qs))
data.frame(
x = dens$x,
density = dens$y,
scaled = dens$y / max(dens$y, na.rm = TRUE),
count = dens$y * nx,
n = nx,
interval_width = interval_width
)
}full_ranef_density <- ranef_samples %>%
group_by(Parameter, PlotSubject, Subject) %>%
do(quick_density_df(.$value, interval_width = 1)) %>%
ungroup() %>%
group_by(Parameter) %>%
mutate(scaled_by_facet = density / max(density)) %>%
ungroup()For some reason, geom_joy() doesn’t work. But the package provides geom_ridgeline().
ggplot(full_ranef_density) +
aes(x = x, y = PlotSubject, height = scaled_by_facet * 1.8,
group = interaction(PlotSubject, Parameter)) +
scale_y_discrete(expand=c(0.01, 0)) +
scale_x_continuous(expand=c(0, 0)) +
facet_grid(~ Parameter, scales = "free_x") +
labs(y = "Subject", x = NULL) +
geom_ridgeline()
We could have the filled portion define the boundaries of, say, the 90% interval.
inner_ranef_density <- ranef_samples %>%
group_by(Parameter, PlotSubject, Subject) %>%
do(quick_density_df(.$value, interval_width = .9)) %>%
ungroup() %>%
group_by(Parameter) %>%
mutate(scaled_by_facet = density / max(density)) %>%
ungroup()
ggplot(full_ranef_density) +
aes(x = x, y = PlotSubject, height = scaled_by_facet * 1.8,
group = interaction(PlotSubject, Parameter)) +
scale_y_discrete(expand=c(0.01, 0)) +
scale_x_continuous(expand=c(0, 0)) +
facet_grid(~ Parameter, scales = "free_x") +
labs(y = "Subject", x = NULL) +
geom_ridgeline(fill = NA) +
geom_ridgeline(data = inner_ranef_density)
devtools::session_info()
#> Session info -------------------------------------------------------------
#> setting value
#> version R version 3.4.0 (2017-04-21)
#> system x86_64, mingw32
#> ui RTerm
#> language (EN)
#> collate English_United States.1252
#> tz America/Chicago
#> date 2017-07-12
#> Packages -----------------------------------------------------------------
#> package * version date source
#> abind 1.4-5 2016-07-21 CRAN (R 3.4.0)
#> acepack 1.4.1 2016-10-29 CRAN (R 3.4.0)
#> arm 1.9-3 2016-11-27 CRAN (R 3.4.1)
#> assertthat 0.2.0 2017-04-11 CRAN (R 3.4.0)
#> backports 1.1.0 2017-05-22 CRAN (R 3.4.0)
#> base * 3.4.0 2017-04-21 local
#> base64enc 0.1-3 2015-07-28 CRAN (R 3.4.0)
#> bindr 0.1 2016-11-13 CRAN (R 3.4.0)
#> bindrcpp * 0.2 2017-06-17 CRAN (R 3.4.0)
#> checkmate 1.8.3 2017-07-03 CRAN (R 3.4.1)
#> cluster 2.0.6 2017-03-10 CRAN (R 3.4.0)
#> coda 0.19-1 2016-12-08 CRAN (R 3.4.0)
#> colorspace 1.3-2 2016-12-14 CRAN (R 3.4.0)
#> compiler 3.4.0 2017-04-21 local
#> data.table 1.10.4 2017-02-01 CRAN (R 3.4.0)
#> datasets * 3.4.0 2017-04-21 local
#> devtools 1.13.2 2017-06-02 CRAN (R 3.4.0)
#> digest 0.6.12 2017-01-27 CRAN (R 3.4.0)
#> dplyr * 0.7.1 2017-06-22 CRAN (R 3.4.0)
#> evaluate 0.10.1 2017-06-24 CRAN (R 3.4.0)
#> forcats 0.2.0 2017-01-23 CRAN (R 3.4.0)
#> foreign 0.8-67 2016-09-13 CRAN (R 3.4.0)
#> Formula 1.2-1 2015-04-07 CRAN (R 3.4.0)
#> ggjoy * 0.1 2017-07-12 Github (clauswilke/ggjoy@3402761)
#> ggplot2 * 2.2.1 2016-12-30 CRAN (R 3.4.0)
#> ghdown 0.0.0.9000 2017-07-12 Github (tjmahr/ghdown@d479b10)
#> glue 1.1.1 2017-06-21 CRAN (R 3.4.0)
#> graphics * 3.4.0 2017-04-21 local
#> grDevices * 3.4.0 2017-04-21 local
#> grid 3.4.0 2017-04-21 local
#> gridExtra 2.2.1 2016-02-29 CRAN (R 3.4.0)
#> gtable 0.2.0 2016-02-26 CRAN (R 3.4.0)
#> Hmisc 4.0-3 2017-05-02 CRAN (R 3.4.0)
#> htmlTable 1.9 2017-01-26 CRAN (R 3.4.0)
#> htmltools 0.3.6 2017-04-28 CRAN (R 3.4.0)
#> htmlwidgets 0.9 2017-07-10 CRAN (R 3.4.0)
#> knitr 1.16 2017-05-18 CRAN (R 3.4.0)
#> labeling 0.3 2014-08-23 CRAN (R 3.4.0)
#> lattice 0.20-35 2017-03-25 CRAN (R 3.4.0)
#> latticeExtra 0.6-28 2016-02-09 CRAN (R 3.4.0)
#> lazyeval 0.2.0 2016-06-12 CRAN (R 3.4.0)
#> lme4 * 1.1-13 2017-04-19 CRAN (R 3.4.0)
#> magrittr 1.5 2014-11-22 CRAN (R 3.4.0)
#> MASS 7.3-47 2017-02-26 CRAN (R 3.4.0)
#> Matrix * 1.2-9 2017-03-14 CRAN (R 3.4.0)
#> memoise 1.1.0 2017-04-21 CRAN (R 3.4.0)
#> methods * 3.4.0 2017-04-21 local
#> minqa 1.2.4 2014-10-09 CRAN (R 3.4.0)
#> munsell 0.4.3 2016-02-13 CRAN (R 3.4.0)
#> nlme 3.1-131 2017-02-06 CRAN (R 3.4.0)
#> nloptr 1.0.4 2014-08-04 CRAN (R 3.4.0)
#> nnet 7.3-12 2016-02-02 CRAN (R 3.4.0)
#> pkgconfig 2.0.1 2017-03-21 CRAN (R 3.4.0)
#> plyr 1.8.4 2016-06-08 CRAN (R 3.4.0)
#> R6 2.2.2 2017-06-17 CRAN (R 3.4.0)
#> RColorBrewer 1.1-2 2014-12-07 CRAN (R 3.4.0)
#> Rcpp 0.12.11 2017-05-22 CRAN (R 3.4.0)
#> reshape2 1.4.2 2016-10-22 CRAN (R 3.4.0)
#> rlang 0.1.1 2017-05-18 CRAN (R 3.4.0)
#> rmarkdown 1.6 2017-06-15 CRAN (R 3.4.0)
#> rpart 4.1-11 2017-03-13 CRAN (R 3.4.0)
#> rprojroot 1.2 2017-01-16 CRAN (R 3.4.0)
#> scales 0.4.1 2016-11-09 CRAN (R 3.4.0)
#> splines 3.4.0 2017-04-21 local
#> stats * 3.4.0 2017-04-21 local
#> stringi 1.1.5 2017-04-07 CRAN (R 3.4.0)
#> stringr 1.2.0 2017-02-18 CRAN (R 3.4.0)
#> survival 2.41-3 2017-04-04 CRAN (R 3.4.0)
#> tibble * 1.3.3 2017-05-28 CRAN (R 3.4.0)
#> tools 3.4.0 2017-04-21 local
#> utils * 3.4.0 2017-04-21 local
#> withr 1.0.2 2016-06-20 CRAN (R 3.4.0)
#> yaml 2.1.14 2016-11-12 CRAN (R 3.4.0)