30 Missing Data

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

Impute \(M\) data sets.
Fit a model in each of the \(M\) imputed data sets.
Compute marginal effects in each of the \(M\) data sets.
Pool results.

To highlight the workflow, we consider a simple linear regression model, although the same workflow should work with any model type that is fit using a formula interface and a data argument.

marginaleffects directly supports the mice and Amelia imputation packages, as well as any other package that can return a list of imputed data frames. This is demonstrated below using the iris dataset, in which we insert missing observations randomly and then impute missing values using several packages.

library(marginaleffects)
set.seed(1024)

dat <- iris
dat$Sepal.Length[sample(seq_len(nrow(iris)), 40)] <- NA
dat$Sepal.Width[sample(seq_len(nrow(iris)), 40)] <- NA
dat$Species[sample(seq_len(nrow(iris)), 40)] <- NA

30.1 `mice`

First, we impute the dataset using the mice package:

library(mice)

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

Then, we use the standard mice syntax to produce an object of class mira with all the models:

mod_mice <- with(dat_mice, lm(Petal.Width ~ Sepal.Length * Sepal.Width + Species))

Finally, we feed the mira object to a marginaleffects function:

mfx_mice <- avg_slopes(mod_mice, by = "Species")
mfx_mice
#> 
#>          Term            Contrast    Species Estimate Std. Error      t Pr(>|t|)    S   2.5 % 97.5 %   Df
#>  Sepal.Length dY/dX               setosa       0.0684     0.0560  1.222  0.22747  2.1 -0.0440  0.181 49.9
#>  Sepal.Length dY/dX               versicolor   0.0540     0.0558  0.968  0.33850  1.6 -0.0585  0.166 42.6
#>  Sepal.Length dY/dX               virginica    0.0582     0.0512  1.137  0.26149  1.9 -0.0449  0.161 44.8
#>  Sepal.Width  dY/dX               setosa       0.1890     0.0836  2.260  0.02632  5.2  0.0228  0.355 87.0
#>  Sepal.Width  dY/dX               versicolor   0.2092     0.0772  2.710  0.00977  6.7  0.0533  0.365 41.1
#>  Sepal.Width  dY/dX               virginica    0.2242     0.1041  2.155  0.03896  4.7  0.0121  0.436 31.3
#>  Species      versicolor - setosa setosa       1.1399     0.0977 11.668  < 0.001 49.7  0.9435  1.336 48.6
#>  Species      virginica - setosa  setosa       1.7408     0.1108 15.709  < 0.001 67.5  1.5182  1.963 50.3
#>  Species      versicolor - setosa versicolor   1.1399     0.0977 11.668  < 0.001 49.7  0.9435  1.336 48.6
#>  Species      virginica - setosa  versicolor   1.7408     0.1108 15.709  < 0.001 67.5  1.5182  1.963 50.3
#>  Species      versicolor - setosa virginica    1.1399     0.0977 11.668  < 0.001 49.7  0.9435  1.336 48.6
#>  Species      virginica - setosa  virginica    1.7408     0.1108 15.709  < 0.001 67.5  1.5182  1.963 50.3
#> 
#> Type:  response

30.2 `Amelia`

With Amelia, the workflow is essentially the same. First, we impute using Amelia:

library(Amelia)

dat_amelia <- amelia(dat, m = 20, noms = "Species", p2s = 0)

Then, we use Amelia syntax to produce an object of class amest with all the models:

mod_amelia <- with(dat_amelia, lm(Petal.Width ~ Sepal.Length * Sepal.Width + Species))

Finally, we feed the amest object to a marginaleffects function:

mfx_amelia <- avg_slopes(mod_amelia, by = "Species")
mfx_amelia
#> 
#>          Term            Contrast    Species Estimate Std. Error      t Pr(>|t|)    S  2.5 % 97.5 %   Df
#>  Sepal.Length dY/dX               setosa       0.3878     0.0907  4.278  < 0.001 11.8  0.200 0.5753 23.1
#>  Sepal.Length dY/dX               versicolor   0.3231     0.0802  4.030  < 0.001 11.4  0.159 0.4872 28.8
#>  Sepal.Length dY/dX               virginica    0.3467     0.0799  4.340  < 0.001 12.1  0.182 0.5118 23.5
#>  Sepal.Width  dY/dX               setosa      -0.2079     0.1491 -1.395  0.17399  2.5 -0.513 0.0973 28.4
#>  Sepal.Width  dY/dX               versicolor  -0.1157     0.1168 -0.991  0.33068  1.6 -0.355 0.1239 26.9
#>  Sepal.Width  dY/dX               virginica   -0.0452     0.1272 -0.355  0.72426  0.5 -0.303 0.2122 38.9
#>  Species      versicolor - setosa setosa       0.6127     0.1731  3.541  0.00217  8.8  0.251 0.9748 19.1
#>  Species      virginica - setosa  setosa       1.0364     0.2004  5.171  < 0.001 13.8  0.615 1.4582 17.6
#>  Species      versicolor - setosa versicolor   0.6127     0.1731  3.541  0.00217  8.8  0.251 0.9748 19.1
#>  Species      virginica - setosa  versicolor   1.0364     0.2004  5.171  < 0.001 13.8  0.615 1.4582 17.6
#>  Species      versicolor - setosa virginica    0.6127     0.1731  3.541  0.00217  8.8  0.251 0.9748 19.1
#>  Species      virginica - setosa  virginica    1.0364     0.2004  5.171  < 0.001 13.8  0.615 1.4582 17.6
#> 
#> Type:  response

30.3 Other imputation packages: `missRanger`, or lists of imputed data frames.

Several R packages can impute missing data. Indeed, the Missing Data CRAN View lists at least a dozen alternatives. Since user interfaces change a lot from package to package, marginaleffects supports a single workflow that can be used, with some adaptation, with all imputation packages:

Use an external package to create a list of imputed data frames.
Apply the datalist2mids() function from the miceadds package to convert the list of imputed data frames to a mids object.
Use the with() function to fit models to create mira object, as illustrated in the mice and Amelia sections above.
Pass the mira object to a marginaleffects function.

Consider the imputation package missRanger, which generates a list of imputed datasets:

library(miceadds)
library(missRanger)

## convert lists of imputed datasets to `mids` objects
dat_missRanger <- replicate(20, missRanger(dat, verbose = 0), simplify = FALSE)
mids_missRanger <- datlist2mids(dat_missRanger)

## fit models
mod_missRanger <- with(mids_missRanger, lm(Petal.Width ~ Sepal.Length * Sepal.Width + Species))

## `missRanger` slopes
mfx_missRanger <- avg_slopes(mod_missRanger, by = "Species")
mfx_missRanger
#> 
#>          Term            Contrast    Species Estimate Std. Error     t Pr(>|t|)     S   2.5 % 97.5 %  Df
#>  Sepal.Length dY/dX               setosa       0.0586     0.0434  1.35  0.17907   2.5 -0.0272  0.144 142
#>  Sepal.Length dY/dX               versicolor   0.0615     0.0392  1.57  0.11884   3.1 -0.0160  0.139 142
#>  Sepal.Length dY/dX               virginica    0.0605     0.0373  1.62  0.10724   3.2 -0.0133  0.134 142
#>  Sepal.Width  dY/dX               setosa       0.2302     0.0725  3.17  0.00185   9.1  0.0868  0.374 142
#>  Sepal.Width  dY/dX               versicolor   0.2260     0.0561  4.03  < 0.001  13.4  0.1152  0.337 142
#>  Sepal.Width  dY/dX               virginica    0.2229     0.0688  3.24  0.00148   9.4  0.0869  0.359 142
#>  Species      versicolor - setosa setosa       1.1629     0.0725 16.04  < 0.001 109.4  1.0196  1.306 142
#>  Species      virginica - setosa  setosa       1.7867     0.0849 21.03  < 0.001 148.5  1.6188  1.955 142
#>  Species      versicolor - setosa versicolor   1.1629     0.0725 16.04  < 0.001 109.4  1.0196  1.306 142
#>  Species      virginica - setosa  versicolor   1.7867     0.0849 21.03  < 0.001 148.5  1.6188  1.955 142
#>  Species      versicolor - setosa virginica    1.1629     0.0725 16.04  < 0.001 109.4  1.0196  1.306 142
#>  Species      virginica - setosa  virginica    1.7867     0.0849 21.03  < 0.001 148.5  1.6188  1.955 142
#> 
#> Type:  response

30.4 Comparing results with different imputation software

We can use the modelsummary package to compare the results with listwise deletion to the results using different imputations software:

library(modelsummary)

## listwise deletion slopes
mod_lwd <- lm(Petal.Width ~ Sepal.Length * Sepal.Width + Species, data = dat)
mfx_lwd <- avg_slopes(mod_lwd, by = "Species")

## regression table
models <- list(
    "LWD" = mfx_lwd,
    "mice" = mfx_mice,
    "Amelia" = mfx_amelia,
    "missRanger" = mfx_missRanger)
modelsummary(models, shape = term : contrast + Species ~ model)

	Species	LWD	mice	Amelia	missRanger
Sepal.Length dY/dX	setosa	0.033	0.068	0.388	0.059
		(0.061)	(0.056)	(0.091)	(0.043)
	versicolor	0.050	0.054	0.323	0.061
		(0.061)	(0.056)	(0.080)	(0.039)
	virginica	0.043	0.058	0.347	0.061
		(0.058)	(0.051)	(0.080)	(0.037)
Sepal.Width dY/dX	setosa	0.274	0.189	-0.208	0.230
		(0.091)	(0.084)	(0.149)	(0.073)
	versicolor	0.255	0.209	-0.116	0.226
		(0.074)	(0.077)	(0.117)	(0.056)
	virginica	0.234	0.224	-0.045	0.223
		(0.083)	(0.104)	(0.127)	(0.069)
Species versicolor - setosa	setosa	1.157	1.140	0.613	1.163
		(0.097)	(0.098)	(0.173)	(0.073)
	versicolor	1.157	1.140	0.613	1.163
		(0.097)	(0.098)	(0.173)	(0.073)
	virginica	1.157	1.140	0.613	1.163
		(0.097)	(0.098)	(0.173)	(0.073)
Species virginica - setosa	setosa	1.839	1.741	1.036	1.787
		(0.123)	(0.111)	(0.200)	(0.085)
	versicolor	1.839	1.741	1.036	1.787
		(0.123)	(0.111)	(0.200)	(0.085)
	virginica	1.839	1.741	1.036	1.787
		(0.123)	(0.111)	(0.200)	(0.085)
Num.Obs.		60	150	150	150
Num.Imp.			20	20	20
R2		0.953	0.930	0.853	0.948
R2 Adj.		0.949	0.928	0.848	0.946
AIC		-34.0
BIC		-19.3
Log.Lik.		23.997
F		220.780
RMSE		0.16

30.5 Passing new data arguments: `newdata`, `wts`, `by`, etc.

Sometimes we want to pass arguments changing or specifying the data on which we will do our analysis using marginaleffects. This can be for reasons such as wanting to specify the values or weights at which we evaluate e.g. avg_slopes(), or due to the underlying models not robustly preserving all the original data columns (such as fixest objects not saving their data in the fit object making it potentially challenging to retrieve, and even if retrievable it will not include the weights used during fitting as a column as wts expects when given a string).

If we are not using multiple imputation, or if we want to just pass a single dataset to the several fitted models after multiple imputation, we can pass a single dataset to the newdata argument. However, if we wish to supply each model in our list resulting after multiple imputation with a /different/ dataset on which to calculate results, we cannot use newdata. Instead, in this case it can be useful to revert to a more manual (but still very easy) approach. Here is an example calculating avg_slopes using a different set of weights for each of the fixest models which we fit after multiple imputation.

set.seed(1024)
library(mice)
library(fixest)
library(marginaleffects)

dat <- mtcars

## insert missing values
dat$hp[sample(seq_len(nrow(mtcars)), 10)] <- NA
dat$mpg[sample(seq_len(nrow(mtcars)), 10)] <- NA
dat$gear[sample(seq_len(nrow(mtcars)), 10)] <- NA

## multiple imputation
dat <- mice(dat, m = 5, method = "sample", printFlag = FALSE)
dat <- complete(dat, action = "all")

## fit models
mod <- lapply(dat, \(x) 
    feglm(am ~ mpg * cyl + hp,
        weight = ~gear,
        family = binomial,
        data = x))

## slopes without weights
lapply(seq_along(mod), \(i) 
    avg_slopes(mod[[i]], newdata = dat[[i]])) |>
    mice::pool()
#> Class: mipo    m = 5 
#>   term contrast m     estimate         ubar            b            t dfcom       df      riv    lambda       fmi
#> 1  cyl    dY/dX 5 -0.134280454 7.097466e-04 2.347331e-03 3.526544e-03    27 2.803719 3.968737 0.7987416 0.8680966
#> 2   hp    dY/dX 5  0.001649773 5.709036e-07 1.375452e-06 2.221446e-06    27 3.419709 2.891105 0.7430036 0.8230684
#> 3  mpg    dY/dX 5  0.006082804 1.080647e-04 2.722234e-04 4.347329e-04    27 3.324494 3.022893 0.7514227 0.8300305

## slopes with weights
lapply(seq_along(mod), \(i) 
    avg_slopes(mod[[i]], newdata = dat[[i]], wts = "gear")) |>
    mice::pool()
#> Class: mipo    m = 5 
#>   term contrast m     estimate         ubar            b            t dfcom       df      riv    lambda       fmi
#> 1  cyl    dY/dX 5 -0.135839444 7.281041e-04 2.481021e-03 3.705329e-03    27 2.752189 4.089010 0.8034981 0.8718206
#> 2   hp    dY/dX 5  0.001671173 5.697747e-07 1.424648e-06 2.279352e-06    27 3.340200 3.000446 0.7500278 0.8288809
#> 3  mpg    dY/dX 5  0.006251144 1.056103e-04 2.705239e-04 4.302390e-04    27 3.289588 3.073835 0.7545310 0.8325867

30.1 mice

30.2 Amelia

30.3 Other imputation packages: missRanger, or lists of imputed data frames.

30.4 Comparing results with different imputation software

30.5 Passing new data arguments: newdata, wts, by, etc.

30.1 `mice`

30.2 `Amelia`

30.3 Other imputation packages: `missRanger`, or lists of imputed data frames.

30.5 Passing new data arguments: `newdata`, `wts`, `by`, etc.