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 |