knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(here)
library(brms)
library(brmstools)
library(dplyr)

Introduction

Forest plots display estimated parameters from multiple sources (studies, participants, etc.) in one figure. They are most commonly used in meta-analysis, where individual studies are used to inform an average, or meta-analytic, overall estimate. However, they can be seamlessly applied to other types of multilevel models--models in which parameters are assumed to vary among units. brmstools' forest() function draws forest plots from brmsfit objects. They should be most useful for meta-analytic models, but can be produced from any brmsfit with one or more varying parameters.

The forest() function uses the fantastic ggridges R package in the backend, and assumes you've installed it. If you haven't, forest() will return an error.

Random effects meta-analysis

We illustrate using a data set 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)

brms allows flexible specification of meta-analytic models.

fit_rem <- brm(
  yi | se(sei) ~ 1 + (1|study),
  data = dat,
  cores = 4, 
  control=list(adapt_delta = .99)
)
# Save time by using locally saved models
save(fit_rem, file = here("vignettes/forest-plots/fit_rem.rda"))
load(here("vignettes/forest-plots/fit_rem.rda"))

Use forest() to draw the forest plot:

forest(fit_rem)

There are various options (see ?forest)

forest(fit_rem, 
       level = .80, 
       av_name = "Meta-Analytic\nEstimate", 
       col_ridge = "purple", 
       fill_ridge = "grey90")

Data points can also be shown (note this probably only makes sense with a meta-analytic model):

forest(fit_rem, show_data = T)

Multilevel model

The forest() function can be seamlessly applied to any multilevel model.

We use example data from the lme4 package.

data(sleepstudy, package = "lme4")
head(sleepstudy)

A multilevel model with varying intercepts and slopes (effect of Days):

fit_ml <- brm(
  Reaction ~ Days + (Days|Subject),
  data = sleepstudy,
  cores = 4
)
save(fit_ml, file = here("vignettes/forest-plots/fit_ml.rda"))
load(here("vignettes/forest-plots/fit_ml.rda"))

If there are multiple varying parameters, users can input a variable name:

forest(fit_ml, pars = "Days")

Or let the function automatically draw a plot with all the variables:

forest(fit_ml, digits=0)

You can also turn off the ridgeline plots (densities)

forest(fit_ml, density = F, digits=0)


mvuorre/brmstools documentation built on Oct. 22, 2018, 5:23 p.m.