Skip to contents

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:

  1. Impute n data sets
  2. Fit in each of the n imputed data sets
  3. Marginal effects in each of the n data sets
  4. 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