Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?

Complete documentation here: https://cran.r-project.org/web/packages/brms/vignettes/brms_missings.html

library(tidyverse)
library(brms)
library(mice)
library(Amelia)
library(broom)
library(tictoc)

data("nhanes", package = "mice")
head(nhanes)
#>   age  bmi hyp chl
#> 1   1   NA  NA  NA
#> 2   2 22.7   1 187
#> 3   1   NA   1 187
#> 4   3   NA  NA  NA
#> 5   1 20.4   1 113
#> 6   3   NA  NA 184

# Impute beforehand with mice
tic()
imp_mice <- mice(nhanes, m = 5, print = FALSE)
fit_mice <- brm_multiple(bmi ~ age*chl, data = imp_mice, chains = 4)
toc()
#> 71.602 sec elapsed

# Impute beforehand with Amelia
tic()
imp_amelia <- amelia(x = nhanes, m = 5, p2s = 0)
fit_amelia <- brm_multiple(bmi ~ age*chl, data = unclass(imp_amelia$imputations), chains = 4)
toc()
#> 67.132 sec elapsed

# Work with missing data on-the-fly with brms
# Main model = bmi ~ age * chl
# There's missing data in bmi and chl, so we tell brms to model the missingness with `| mi()`
# and we use a second formula to explain how to model missing chl (with age)
tic()
imp_formula <- bf(bmi | mi() ~ age * mi(chl)) +
  bf(chl | mi() ~ age) + set_rescor(FALSE)
fit_on_the_fly <- brm(imp_formula, data = nhanes)
toc()
#> 79.215 sec elapsed
# 70.686 seconds

estimates <- bind_rows(
  tidy(fit_mice) %>% mutate(model = "mice") %>% filter(str_detect(term, "^b_")),
  tidy(fit_amelia) %>% mutate(model = "amelia") %>% filter(str_detect(term, "^b_")),
  tidy(fit_on_the_fly) %>% mutate(model = "on-the-fly") %>% 
    filter(str_detect(term, "^b_bmi") | term == "bsp_bmi_michl:age")
) %>% 
  mutate(term = str_replace_all(term, "bsp_bmi_michl:age", "b_age:chl"),
         term = str_remove_all(term, "bmi_"))

ggplot(estimates, aes(x = estimate, y = term, color = model, shape = model)) +
  geom_pointrange(aes(xmin = lower, xmax = upper), position = position_dodge(width = 0.3))

Created on 2020-06-29 by the reprex package (v0.3.0)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.