---

First, using the Amelia package, we multiply impute a dataset with missing values:

library(purrr)
library(Amelia)

## Loading required package: Rcpp

## ## 
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.4, built: 2015-12-05)
## ## Copyright (C) 2005-2017 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##

data(africa)
a.out <- amelia(x = africa, cs = "country", ts = "year", logs = "gdp_pc")

## -- Imputation 1 --
## 
##   1  2
## 
## -- Imputation 2 --
## 
##   1  2  3
## 
## -- Imputation 3 --
## 
##   1  2
## 
## -- Imputation 4 --
## 
##   1  2
## 
## -- Imputation 5 --
## 
##   1  2

Then we can use the output object from the Amelia package to fit a list of models:

models <- map(a.out$imputations, ~ lm(gdp_pc ~ trade + civlib, data = .))

Amelia includes a mi.meld function that can be used to combine the results:

mi.meld(reduce(map(models, coef), rbind),
        reduce(map(models, ~ coef(summary(.))[, "Std. Error"]), rbind))

## $q.mi
##      (Intercept)    trade   civlib
## [1,]    115.9486 18.15496 -637.201
## 
## $se.mi
##      (Intercept)    trade   civlib
## [1,]    96.92541 1.239931 182.6399

Nicer output is available from several of the available multiple-imputation packages in R, including the mice and mitools.

library(mice)

## Loading required package: lattice

(foo <- summary(pool(as.mira(models))))

##                    est         se         t       df     Pr(>|t|)
## (Intercept)  115.94859  96.925405  1.196266 114.2037 0.2340703917
## trade         18.15496   1.239931 14.641908 114.7050 0.0000000000
## civlib      -637.20096 182.639877 -3.488838 112.7967 0.0006928665
##                  lo 95      hi 95 nmis        fmi      lambda
## (Intercept)  -76.05622  307.95339   NA 0.02320620 0.006248513
## trade         15.69883   20.61109   NA 0.01972270 0.002778270
## civlib      -999.05056 -275.35137   NA 0.03110488 0.014076359

Quantities of interest for MI

Using simulations to estimate uncertainty makes it pretty easy to incorporate additional uncertainty from imputed values.

library(smargins)

models.sm <- map(models,
                 ~smargins(., trade = quantile(trade,
                                               c(.25, .50, .75),
                                               na.rm = TRUE)))

summary(reduce(models.sm, rbind))

##     trade      mean       sd    median lower_2.5 upper_97.5
## 1 38.4075  629.5738 42.89750  629.5891  544.1707   713.8550
## 2 58.8350  634.5849 42.95185  635.4293  551.0451   719.7549
## 3 80.7875  629.5525 44.51840  629.5309  543.9320   716.1761
## 4 59.1900  999.1990 32.23231  999.1479  936.3776  1062.7405
## 5 81.1050 1011.3282 32.20677 1011.5555  948.0754  1073.4816
## 6 38.7000 1006.1459 32.75386 1006.0958  942.5650  1072.0525
## 7 59.5300 1397.5705 39.95786 1397.6203 1319.2201  1475.5956
## 8 38.4650 1410.9172 39.92553 1410.7654 1331.0978  1487.8401