Multiple Imputation

The marginaleffects package offers convenience functions to compute and display predictions, contrasts, and marginal effects from models with multiple imputation from the mice package. The workflow follows Rubin’s rules (Rubin, 1987, p. 76) that uses the following steps:

Impute n data sets
Fit in each of the n imputed data sets
Marginal effects in each of the n data sets
Pool results

To highlight the workflow, we can consider 3 situations: linear regression, logistic regression, and multilevel models with lme4. To show the linear and logistic models, we are going to use the classic mtcars data set but we will artificially add missing data. For the multilevel models, we will use the sleepstudy data set from lme4.

Linear Regression with `mice`

Let’s first set up the data.

library(mice)
library(marginaleffects)
library(modelsummary)

dat <- mtcars
dat$am[c(2, 5, 9, 12)] <- NA
dat$mpg[c(3, 1, 8, 16)] <- NA
dat$hp[c(1, 10, 13, 18)] <- NA

The next steps are to use mice to create the imputed data sets. Here, we are asking for m = 20 imputations. Importantly, when mice creates the multiple imputation, it creates a list of data sets. So here, dat_mice has 20 nearly identical versions of the data in a list where any missing values were imputed.

dat_mice <- mice(dat, m = 20, printFlag = FALSE, .Random.seed = 1024)
dat_mice <- complete(dat_mice, "all")

To work with the list of data sets, we’ll create a function (in this case called fit_reg) that fits the model and computes the marginal effects.

fit_reg <- function(dat) {
    mod <- lm(mpg ~ hp, data = dat)
    out <- slopes(mod, newdata = dat)
    return(out)
}

Using this function, we can apply it to each data set in the dat_mice list using lapply(). From there, we can pool and get a summary.

mod_imputation <- lapply(dat_mice, fit_reg)
mod_imputation <- pool(mod_imputation)

summary(mod_imputation)
#>   term    contrast    estimate  std.error statistic       df     p.value
#> 1   hp mean(dY/dX) -0.06654044 0.01181385 -5.632407 1386.179 2.14888e-08

We can compare this with what the model would have looked like without any missing data. The estimates are very similar (within one standard error) and the p-value for the imputed models is slightly higher than the full model (as expected).

mod_missing <- lm(mpg ~ hp, data = dat)
mod_missing <- slopes(mod_missing)
mod_complete <- lm(mpg ~ hp, data = mtcars)
mod_complete <- slopes(mod_complete)

models <- list(
    "Listwise Deletion" = mod_missing,
    "Complete" = mod_complete,
    "Multiple Imputation" = mod_imputation)

modelsummary(models)

	Listwise Deletion	Complete	Multiple Imputation
hp	−0.063	−0.068	−0.067
	(0.011)	(0.010)	(0.012)
Num.Obs.	25	32
R2	0.573	0.602
R2 Adj.	0.555	0.589
AIC	143.6	181.2
BIC	147.3	185.6
Log.Lik.	−68.815	−87.619
F	30.881	45.460
RMSE	3.79	3.74

Categories and Contrasts: Problem and Solution

One particular problem arises in the cases of contrasts and categorical predictors. To see it, notice that when there are contrasts or categorical predictors, the tidy() method of marginaleffects identifies unique estimates using two columns called term and contrast:

mod <- lm(mpg ~ factor(cyl), data = dat)
mfx <- slopes(mod)
tidy(mfx)
#> # A tibble: 2 × 9
#>   type     term  contrast       estim…¹ std.e…² stati…³  p.value conf.…⁴ conf.…⁵
#>   <chr>    <chr> <chr>            <dbl>   <dbl>   <dbl>    <dbl>   <dbl>   <dbl>
#> 1 response cyl   mean(6) - mea…   -7.81    1.67   -4.67 2.96e- 6   -11.1   -4.54
#> 2 response cyl   mean(8) - mea…  -11.9     1.38   -8.64 5.55e-18   -14.6   -9.19
#> # … with abbreviated variable names ¹estimate, ²std.error, ³statistic,
#> #   ⁴conf.low, ⁵conf.high

This poses problems because the mice::pool function merges estimates based only on the term column. This means that our original procedure will erroneously combine different contrast levels. For example:

fit_reg <- function(dat) {
    mod <- lm(mpg ~ factor(cyl), data = dat)
    out <- slopes(mod, newdata = dat)
    return(out)
}
mod_imputation <- lapply(dat_mice, fit_reg)
mod_imputation <- pool(mod_imputation)

summary(mod_imputation)
#>   term          contrast   estimate std.error statistic        df      p.value
#> 1  cyl mean(6) - mean(4)  -7.201169  1.733322 -4.154549 1436.2230 3.451518e-05
#> 2  cyl mean(8) - mean(4) -11.601169  1.490108 -7.785456  669.6921 2.644828e-14

One hack to work around this limitation is to assign a custom class to the object and to create a custom tidy method that combines the term and contrast columns:

fit_reg <- function(dat) {
    mod <- lm(mpg ~ factor(cyl), data = dat)
    out <- slopes(mod, newdata = dat)
    # the next line assigns a custom class
    class(out) <- c("custom", class(out))
    return(out)
}

# this custom method will be called automatically for all objects produced by fit_reg()
tidy.custom <- function(x, ...) {
    out <- marginaleffects:::tidy.slopes(x, ...)
    out$term <- paste(out$term, out$contrast)
    return(out)
}

mod_imputation <- lapply(dat_mice, fit_reg)
mod_imputation <- pool(mod_imputation)

summary(mod_imputation)
#>                    term          contrast   estimate std.error statistic
#> 1 cyl mean(6) - mean(4) mean(6) - mean(4)  -7.201169  1.733322 -4.154549
#> 2 cyl mean(8) - mean(4) mean(8) - mean(4) -11.601169  1.490108 -7.785456
#>          df      p.value
#> 1 1436.2230 3.451518e-05
#> 2  669.6921 2.644828e-14

Logistic Regression with `mice`

For the logistic regression, we’ll work with the same dat_mice imputed data sets. We’ll update our function to run the logistic regression that we want and call it fit_logistic.

fit_logistic <- function(dat) {
    mod <- glm(am ~ mpg, data = dat, family = binomial)
    out <- slopes(mod, newdata = dat)
    return(out)
}

Using this function, we can apply it to each data set in the dat_mice list using lapply(). From there, we can pool and get a summary.

mod_imputation <- lapply(dat_mice, fit_logistic)
mod_imputation <- pool(mod_imputation)

summary(mod_imputation)
#>   term    contrast   estimate  std.error statistic      df      p.value
#> 1  mpg mean(dY/dX) 0.04861021 0.01088347  4.466427 204.956 1.313484e-05

Again, we can compare this with what the model would have looked like without any missing data. The estimates, again, are very similar (within one standard error) and the p-value for the imputed models is slightly higher than the full model (as expected).

mod_missing <- glm(am ~ mpg, data = dat, family = binomial)
mod_complete <- glm(am ~ mpg, data = mtcars, family = binomial)
mod_missing <- slopes(mod_missing)
mod_complete <- slopes(mod_complete)

models <- list(
    "Listwise Deletion" = mod_missing,
    "Complete" = mod_complete,
    "Multiple Imputation" = mod_imputation)

modelsummary(models)

	Listwise Deletion	Complete	Multiple Imputation
mpg	0.044	0.046	0.049
	(0.009)	(0.009)	(0.011)
Num.Obs.	24	32
AIC	24.0	33.7
BIC	26.4	36.6
Log.Lik.	−10.023	−14.838
F	5.702	7.148
RMSE	0.37	0.39

Multilevel Modeling with `lme4`

Our last example with use data from lme4 known as sleepstudy. Let’s first set up the data. We randomly create missing in the outcome variable known as Reaction.

library(lme4)
data("sleepstudy")

set.seed(1234)

dat2 <- sleepstudy
dat2$Reaction[sample(1:180, 10)] <- NA

As before, the next steps are to use mice to create the imputed data sets.

dat_mice2 <- mice(dat2, m = 20, printFlag = FALSE, .Random.seed = 1024)
dat_mice2 <- complete(dat_mice2, "all")

To work with the list of data sets, we’ll create a function (in this case called fit_reg) that fits the model and computes the marginal effects.

fit_mlm <- function(dat) {
    mod <- lmer(Reaction ~ Days + (1 + Days|Subject), data = dat)
    out <- slopes(mod, newdata = dat)
    return(out)
}

Using this function, we can apply it to each data set in the dat_mice list using lapply(). From there, we can pool and get a summary.

mod_imputation <- lapply(dat_mice2, fit_mlm)
mod_imputation <- pool(mod_imputation)

summary(mod_imputation)
#>   term    contrast estimate std.error statistic       df      p.value
#> 1 Days mean(dY/dX) 10.25971  1.563251  6.563063 172.3149 5.950207e-10

Like the previous models, we can compare this with what the model would have looked like without any missing data. The estimates are very similar (within one standard error) and the p-value for the imputed models is slightly higher than the full model (as expected).

mod_complete <- lmer(Reaction ~ Days + (1 + Days|Subject), data = sleepstudy)
mod_missing <- lmer(Reaction ~ Days + (1 + Days|Subject), data = dat2)
mod_complete <- slopes(mod_complete)
mod_missing <- slopes(mod_missing)

models <- list(
    "Listwise Deletion" = mod_missing,
    "Complete" = mod_complete,
    "Multiple Imputation" = mod_imputation)

modelsummary(models)

	Listwise Deletion	Complete	Multiple Imputation
Days	10.407	10.467	10.260
	(1.622)	(1.546)	(1.563)
Num.Obs.	170	180	180
Num.Imp.			20
R2 Marg.	0.274	0.279
R2 Cond.	0.796	0.799
AIC	1664.9	1755.6
BIC	1683.7	1774.8
ICC	0.7	0.7
RMSE	23.55	23.44

Linear Regression with mice

Categories and Contrasts: Problem and Solution

Logistic Regression with mice

Multilevel Modeling with lme4

Linear Regression with `mice`

Logistic Regression with `mice`

Multilevel Modeling with `lme4`