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`

.

`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 <- marginaleffects(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 estimate std.error statistic df p.value
#> 1 hp -0.0653413 0.01133834 -5.762865 5372.835 8.729996e-09
```

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 <- marginaleffects(mod_missing)
mod_complete <- lm(mpg ~ hp, data = mtcars)
mod_complete <- marginaleffects(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.065 |

(0.011) | (0.010) | (0.011) | |

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 |

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 <- marginaleffects(mod)
tidy(mfx)
#> type term contrast estimate std.error statistic p.value conf.low
#> 1 response cyl 6 - 4 -7.811111 1.671348 -4.673540 2.960519e-06 -11.08689
#> 2 response cyl 8 - 4 -11.882906 1.375107 -8.641441 5.550854e-18 -14.57807
#> conf.high
#> 1 -4.535330
#> 2 -9.187746
```

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 <- marginaleffects(mod, newdata = dat)
return(out)
}
mod_imputation <- lapply(dat_mice, fit_reg)
mod_imputation <- pool(mod_imputation)
summary(mod_imputation)
#> term estimate std.error statistic df p.value
#> 1 cyl -8.983929 2.851044 -3.151101 77.89198 0.002308526
```

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 <- marginaleffects(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.marginaleffects(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 estimate std.error statistic df p.value
#> 1 cyl 6 - 4 -6.702857 1.782722 -3.759900 1381.464 1.770773e-04
#> 2 cyl 8 - 4 -11.265000 1.472236 -7.651627 1856.132 3.153033e-14
```

`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 <- marginaleffects(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 estimate std.error statistic df p.value
#> 1 mpg 0.04939145 0.01036178 4.766693 345.6361 2.763556e-06
```

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 <- marginaleffects(mod_missing)
mod_complete <- marginaleffects(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.010) | |

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 |

`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 <- marginaleffects(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 estimate std.error statistic df p.value
#> 1 Days 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 <- marginaleffects(mod_complete)
mod_missing <- marginaleffects(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 |