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:
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 |