Imputation before model fitting
There are many approaches allowing us to impute missing data before
the actual model fitting takes place. From a statistical perspective,
multiple imputation is one of the best solutions. Each missing value is
not imputed once but m times leading to a total of
m fully imputed data sets. The model can then be fitted to
each of those data sets separately and results are pooled across models,
afterwards. One widely applied package for multiple imputation is
mice (Buuren & Groothuis-Oudshoorn, 2010) and we
will use it in the following in combination with brms.
Here, we apply the default settings of mice, which
means that all variables will be used to impute missing values in all
other variables and imputation functions automatically chosen based on
the variables’ characteristics.
Now, we have m = 5 imputed data sets stored within the
imp object. In practice, we will likely need more than
5 of those to accurately account for the uncertainty
induced by the missingness, perhaps even in the area of 100
imputed data sets (Zhou & Reiter, 2010). Of course, this increases
the computational burden by a lot and so we stick to m = 5
for the purpose of this vignette. Regardless of the value of
m, we can either extract those data sets and then pass them
to the actual model fitting function as a list of data frames, or pass
imp directly. The latter works because
brms offers special support for data imputed by
mice. We will go with the latter approach, since it is
less typing. Fitting our model of interest with brms to
the multiple imputed data sets is straightforward.
The returned fitted model is an ordinary brmsfit object
containing the posterior draws of all m submodels. While
pooling across models is not necessarily straightforward in classical
statistics, it is trivial in a Bayesian framework. Here, pooling results
of multiple imputed data sets is simply achieved by combining the
posterior draws of the submodels. Accordingly, all post-processing
methods can be used out of the box without having to worry about pooling
at all.
Family: gaussian
Links: mu = identity; sigma = identity
Formula: bmi ~ age * chl
Data: imp (Number of observations: 25)
Draws: 10 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 10000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 14.83 7.99 -0.55 30.78 1.10 64 552
age 1.02 4.93 -8.53 10.82 1.14 46 247
chl 0.09 0.04 0.01 0.17 1.11 57 444
age:chl -0.02 0.02 -0.07 0.02 1.12 54 247
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 3.33 0.57 2.42 4.65 1.11 61 317
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).In the summary output, we notice that some Rhat values
are higher than \(1.1\) indicating
possible convergence problems. For models based on multiple imputed data
sets, this is often a false positive: Chains of
different submodels may not overlay each other exactly, since there were
fitted to different data. We can see the chains on the right-hand side
of
Such non-overlaying chains imply high Rhat values
without there actually being any convergence issue. Accordingly, we have
to investigate the convergence of the submodels separately, which we can
do for example via:
library(posterior)
draws <- as_draws_array(fit_imp1)
# every dataset has nc = 2 chains in this example
nc <- nchains(fit_imp1) / m
draws_per_dat <- lapply(1:m,
\(i) subset_draws(draws, chain = ((i-1)*nc+1):(i*nc))
)
lapply(draws_per_dat, summarise_draws, default_convergence_measures())[[1]]
# A tibble: 8 × 4
variable rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl>
1 b_Intercept 1.01 703. 923.
2 b_age 1.00 667. 903.
3 b_chl 1.00 716. 907.
4 b_age:chl 1.00 649. 791.
5 sigma 1.00 1049. 1191.
6 Intercept 1.00 1140. 1047.
7 lprior 1.00 973. 1186.
8 lp__ 1.00 752. 1380.
[[2]]
# A tibble: 8 × 4
variable rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl>
1 b_Intercept 1.00 746. 737.
2 b_age 1.00 736. 940.
3 b_chl 1.00 774. 703.
4 b_age:chl 1.00 731. 821.
5 sigma 1.00 1172. 1234.
6 Intercept 1.00 1345. 837.
7 lprior 1.00 890. 987.
8 lp__ 1.00 469. 671.
[[3]]
# A tibble: 8 × 4
variable rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl>
1 b_Intercept 1.01 510. 495.
2 b_age 1.01 497. 754.
3 b_chl 1.01 563. 552.
4 b_age:chl 1.01 479. 681.
5 sigma 1.00 1073. 1038.
6 Intercept 1.00 1243. 841.
7 lprior 1.01 1010. 834.
8 lp__ 1.00 728. 930.
[[4]]
# A tibble: 8 × 4
variable rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl>
1 b_Intercept 1.00 732. 871.
2 b_age 1.00 743. 936.
3 b_chl 1.00 765. 916.
4 b_age:chl 1.00 703. 914.
5 sigma 1.00 937. 1080.
6 Intercept 1.01 1435. 1177.
7 lprior 1.00 1018. 1140.
8 lp__ 1.01 622. 761.
[[5]]
# A tibble: 8 × 4
variable rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl>
1 b_Intercept 1.00 574. 851.
2 b_age 1.00 546. 828.
3 b_chl 1.00 599. 806.
4 b_age:chl 1.00 525. 759.
5 sigma 1.00 885. 1042.
6 Intercept 1.00 1160. 1211.
7 lprior 1.01 815. 1067.
8 lp__ 1.01 736. 1256.The convergence of each of the submodels looks good. Accordingly, we
can proceed with further post-processing and interpretation of the
results. For instance, we could investigate the combined effect of
age and chl.
To summarize, the advantages of multiple imputation are obvious: One can apply it to all kinds of models, since model fitting functions do not need to know that the data sets were imputed, beforehand. Also, we do not need to worry about pooling across submodels when using fully Bayesian methods. The only drawback is the amount of time required for model fitting. Estimating Bayesian models is already quite slow with just a single data set and it only gets worse when working with multiple imputation.
Compatibility with other multiple imputation packages
brms offers built-in support for
mice mainly because I use the latter in some of my own
research projects. Nevertheless, brm_multiple supports all
kinds of multiple imputation packages as it also accepts a list
of data frames as input for its data argument. Thus, you
just need to extract the imputed data frames in the form of a list,
which can then be passed to brm_multiple. Most multiple
imputation packages have some built-in functionality for this task. When
using the mi package, for instance, you simply need to
call the mi::complete function to get the desired
output.
