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)