README.md
In mvuorre/brmstools: Tools and Helpers for brms Package

brmstools

brmstools was an R package that provided one-liners for drawing figures from regression models fitted with the brms package.

The package is no longer maintained: Its functionality can be replicated by using native functions from the brms package and functions from the tidybayes package. Example code for reproducing brmstools’ figures is shown below.

Required packages

library(brms)
#> Loading required package: Rcpp
#> Registered S3 method overwritten by 'xts':
#>   method     from
#>   as.zoo.xts zoo
#> Loading 'brms' package (version 2.9.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
library(tidyverse)
#> ── Attaching packages ──────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
#> ✔ ggplot2 3.2.0          ✔ purrr   0.3.2     
#> ✔ tibble  2.1.3          ✔ dplyr   0.8.3     
#> ✔ tidyr   0.8.3.9000     ✔ stringr 1.4.0     
#> ✔ readr   1.3.1          ✔ forcats 0.4.0
#> ── Conflicts ─────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag()    masks stats::lag()
library(tidybayes)  # https://github.com/mjskay/tidybayes
library(ggridges)  # https://github.com/clauswilke/ggridges
#> 
#> Attaching package: 'ggridges'
#> The following object is masked from 'package:ggplot2':
#> 
#>     scale_discrete_manual

Forest plots

After fitting a meta-analytic model, creating a forest plot is easy with tidybayes and ggridges.

First, fit a meta-analytic model with some example data from the metafor package.

data("dat.bangertdrowns2004", package = "metafor")
dat <- dat.bangertdrowns2004 %>%
  mutate(study = paste0(author, " (", year, ")"), sei = sqrt(vi)) %>%
  select(study, yi, sei) %>%
  slice(1:15)
dat$study <- str_replace(dat$study, ",", "")  # Remove commas in study names
fit_rem <- brm(
  yi | se(sei) ~ 1 + (1|study),
  data = dat,
  cores = 4, 
  control=list(adapt_delta = .99),
  file = "fit-rem"
)

For an explanation of tidybayes::spread_draws(), refer to http://mjskay.github.io/tidybayes/articles/tidy-brms.html.

# Study-specific effects are deviations + average
out_r <- spread_draws(fit_rem, r_study[study,term], b_Intercept) %>% 
  mutate(b_Intercept = r_study + b_Intercept) 
# Average effect
out_f <- spread_draws(fit_rem, b_Intercept) %>% 
  mutate(study = "Average")
# Combine average and study-specific effects' data frames
out_all <- bind_rows(out_r, out_f) %>% 
  ungroup() %>%
  # Ensure that Average effect is on the bottom of the forest plot
  mutate(study = fct_relevel(study, "Average"))
# Data frame of summary numbers
out_all_sum <- group_by(out_all, study) %>% 
      mean_qi(b_Intercept)
#> Warning: unnest() has a new interface. See ?unnest for details.
#> Try `cols = c(.lower, .upper)`, with `mutate()` needed
# Draw plot
out_all %>%   
  ggplot(aes(b_Intercept, study)) +
  geom_density_ridges(
    rel_min_height = 0.01, 
    col = NA,
    scale = 1
  ) +
  geom_pointintervalh(
    data = out_all_sum, size = 1
  ) +
  geom_text(
    data = mutate_if(out_all_sum, is.numeric, round, 2),
    # Use glue package to combine strings
    aes(label = glue::glue("{b_Intercept} [{.lower}, {.upper}]"), x = Inf),
    hjust = "inward"
  )
#> Picking joint bandwidth of 0.0226

Panel plots

First, fit a multilevel model to multiple subjects’ data

data(sleepstudy, package = "lme4")
fit_ml <- brm(
  Reaction ~ Days + (Days|Subject),
  data = sleepstudy,
  cores = 4,
  file = "fit-ml"
)

Draw a panel plot using brms::marginal_effects()

marginal_effects(
  fit_ml,
  "Days", 
  conditions = distinct(sleepstudy, Subject),
  re_formula = NULL
  )

Spaghetti plots

Spaghetti plot combines marginal_effects() of the random and fixed effects’ regression lines.

out_f <- marginal_effects(
  fit_ml,
  "Days"
  )[[1]]
out_r <- marginal_effects(
  fit_ml,
  "Days", 
  conditions = distinct(sleepstudy, Subject),
  re_formula = NULL
  )[[1]]
out_r %>% 
  ggplot(aes(Days, estimate__)) +
  geom_ribbon(
    data = out_f,
    aes(ymin = lower__, ymax = upper__),
    alpha = .33
  ) +
  geom_line(data = out_f, size = 2) +
  geom_line(aes(group = Subject))