This tutorial covers techniques of multiple imputation. Multiple imputation is a strategy for dealing with missing data. Whereas we typically (i.e., automatically) deal with missing data through casewise deletion of any observations that have missing values on key variables, imputation attempts to replace missing values with an estimated value. In single imputation, we guess that missing value one time (perhaps based on the means of observed values, or a random sampling of those values). In multiple imputation, we instead draw multiple values for each missing value, effectively building multiple datasets, each of which replaces the missing data in a different way. There are numerous algorithms for this, each of which builds those multiple datasets in different ways. We're not going to discuss the details here, but instead focus on executing multiple imputation in R. The main challenge of multiple imputation is not the analysis (it simply proceeds as usual on each imputed dataset) but instead the aggregation of those separate analyses. The examples below discuss how to do this.
To get a basic feel for the process, let's imagine that we're trying to calculate the mean of a vector of values that contains missing values. We can impute the missing values by drawing from the observed values, repeat the process several times, and then average across the estimated means to get an estimate of the mean with a measure of uncertainty that accounts for the uncertainty due to imputation. Let's create a vector of ten values, seven of which we observe and three of which are missing, and imagine that they are random draws from the population whose mean we're trying to estimate:
set.seed(10)
x <- c(sample(1:10, 7, TRUE), rep(NA, 3))
x
## [1] 6 4 5 7 1 3 3 NA NA NA
We can find the mean using case deletion:
mean(x, na.rm = TRUE)
## [1] 4.143
Our estimate of the sample standard error is then:
sd(x, na.rm = TRUE)/sqrt(sum(!is.na(x)))
## [1] 0.7693
Now let's impute several times to generate a list of imputed vectors:
imp <- replicate(15, c(x[!is.na(x)], sample(x[!is.na(x)], 3, TRUE)), simplify = FALSE)
imp
## [[1]]
## [1] 6 4 5 7 1 3 3 4 1 7
##
## [[2]]
## [1] 6 4 5 7 1 3 3 1 7 6
##
## [[3]]
## [1] 6 4 5 7 1 3 3 1 5 7
##
## [[4]]
## [1] 6 4 5 7 1 3 3 6 4 5
##
## [[5]]
## [1] 6 4 5 7 1 3 3 3 3 1
##
## [[6]]
## [1] 6 4 5 7 1 3 3 3 5 5
##
## [[7]]
## [1] 6 4 5 7 1 3 3 1 3 4
##
## [[8]]
## [1] 6 4 5 7 1 3 3 3 5 7
##
## [[9]]
## [1] 6 4 5 7 1 3 3 6 4 3
##
## [[10]]
## [1] 6 4 5 7 1 3 3 5 3 3
##
## [[11]]
## [1] 6 4 5 7 1 3 3 3 1 7
##
## [[12]]
## [1] 6 4 5 7 1 3 3 4 4 6
##
## [[13]]
## [1] 6 4 5 7 1 3 3 3 4 4
##
## [[14]]
## [1] 6 4 5 7 1 3 3 6 7 6
##
## [[15]]
## [1] 6 4 5 7 1 3 3 3 5 3
The result is a list of five vectors. The first seven values of each is the same as our original data, but the three missing values have been replaced with different combinations of the observed values. To get our new estimated maen, we simply take the mean of each vector, and then average across them:
means <- sapply(imp, mean)
means
## [1] 4.1 4.3 4.2 4.4 3.6 4.2 3.7 4.4 4.2 4.0 4.0 4.3 4.0 4.8 4.0
grandm <- mean(means)
grandm
## [1] 4.147
The result is 4.147, about the same as our original estimate. To get the standard error of our multiple imputation estimate, we need to combine the standard errors of each of our estimates, so that means we need to start by getting the SEs of each imputed vector:
ses <- sapply(imp, sd)/sqrt(10)
Aggregating the standard errors is a bit complicated, but basically sums the mean of the SEs (i.e., the “within-imputation variance”) with the variance across the different estimated means (the “between-imputation variance”). To calculate the within-imputation variance, we simply average the SE estimates:
within <- mean(ses)
To calculate the between-imputation variance, we calculate the sum of squared deviations of each imputed mean from the grand mean estimate:
between <- sum((means - grandm)^2)/(length(imp) - 1)
Then we sum the within- and between-imputation variances (multiply the latter by a small correction):
grandvar <- within + ((1 + (1/length(imp))) * between)
grandse <- sqrt(grandvar)
grandse
## [1] 0.8387
The resulting standard error is interesting because we increase the precision of our estimate by using 10 rather than 7 values (and standard errors are proportionate to sample size), but is larger than our original standard error because we have to account for uncertainty due to imputation. Thus if our missing values are truly missing at random, we can get a better estimate that is actually representative of our original population. Most multiple imputation algorithms are, however, applied to multivariate data rather than a single data vector and thereby use additional information about the relationship between observed values and missingness to reach even more precise estimates of target parameters.
There are three main R packages that offer multiple imputation techniques. Several other packages - described in the OfficialStatistics Task View - supply other imputation techniques, but packages Amelia (by Gary King and collaborators), mi (by Andrew Gelman and collaborators), and mice (by Stef van Buuren and collaborators) provide more than enough to work with. Let's start by installing these packages:
install.packages(c("Amelia", "mi", "mice"), repos = "http://cran.r-project.org")
## Warning: packages 'Amelia', 'mi', 'mice' are in use and will not be
## installed
Now, let's consider an imputation situation where we plan to conduct a regression analysis predicting y by two covariates: x1 and x2 but we have missing data in x1 and x2. Let's start by creating the dataframe:
x1 <- runif(100, 0, 5)
x2 <- rnorm(100)
y <- x1 + x2 + rnorm(100)
mydf <- cbind.data.frame(x1, x2, y)
Now, let's randomly remove some of the observed values of the independent variables:
mydf$x1[sample(1:nrow(mydf), 20, FALSE)] <- NA
mydf$x2[sample(1:nrow(mydf), 10, FALSE)] <- NA
The result is the removal of thirty values, 20 from x1 and 10 from x2:
summary(mydf)
## x1 x2 y
## Min. :0.098 Min. :-2.321 Min. :-1.35
## 1st Qu.:1.138 1st Qu.:-0.866 1st Qu.: 1.17
## Median :2.341 Median : 0.095 Median : 2.39
## Mean :2.399 Mean :-0.038 Mean : 2.28
## 3rd Qu.:3.626 3rd Qu.: 0.724 3rd Qu.: 3.69
## Max. :4.919 Max. : 2.221 Max. : 6.26
## NA's :20 NA's :10
If we estimate the regression on these data, R will force casewise deletion of 28 cases:
lm <- lm(y ~ x1 + x2, data = mydf)
summary(lm)
##
## Call:
## lm(formula = y ~ x1 + x2, data = mydf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5930 -0.7222 0.0018 0.7140 2.4878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0259 0.2196 0.12 0.91
## x1 0.9483 0.0824 11.51 < 2e-16 ***
## x2 0.7487 0.1203 6.23 3.3e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.969 on 69 degrees of freedom
## (28 observations deleted due to missingness)
## Multiple R-squared: 0.706, Adjusted R-squared: 0.698
## F-statistic: 82.9 on 2 and 69 DF, p-value: <2e-16
We should thus be quite skeptical of our results given taht we're discarding a substantial portion of our observations (28%, in fact). Let's see how the various multiple imputation packages address this and affect our inference.
library(Amelia)
imp.amelia <- amelia(mydf)
## -- Imputation 1 --
##
## 1 2 3 4 5 6 7 8 9
##
## -- Imputation 2 --
##
## 1 2 3 4 5 6 7
##
## -- Imputation 3 --
##
## 1 2 3 4 5 6
##
## -- Imputation 4 --
##
## 1 2 3 4 5 6
##
## -- Imputation 5 --
##
## 1 2 3 4 5 6 7
Once we've run our multiple imputation, we can see where are missing data lie:
missmap(imp.amelia)
We can also run our regression model on each imputed dataset. We'll use the lapply function to do this quickly on each of the imputed dataframes:
lm.amelia.out <- lapply(imp.amelia$imputations, function(i) lm(y ~ x1 + x2,
data = i))
If we look at lm.amelia.out we'll see the results of the model run on each imputed dataframe separately:
lm.amelia.out
## $imp1
##
## Call:
## lm(formula = y ~ x1 + x2, data = i)
##
## Coefficients:
## (Intercept) x1 x2
## 0.247 0.854 0.707
##
##
## $imp2
##
## Call:
## lm(formula = y ~ x1 + x2, data = i)
##
## Coefficients:
## (Intercept) x1 x2
## 0.164 0.931 0.723
##
##
## $imp3
##
## Call:
## lm(formula = y ~ x1 + x2, data = i)
##
## Coefficients:
## (Intercept) x1 x2
## 0.0708 0.9480 0.8234
##
##
## $imp4
##
## Call:
## lm(formula = y ~ x1 + x2, data = i)
##
## Coefficients:
## (Intercept) x1 x2
## 0.0656 0.9402 0.6446
##
##
## $imp5
##
## Call:
## lm(formula = y ~ x1 + x2, data = i)
##
## Coefficients:
## (Intercept) x1 x2
## -0.064 0.956 0.820
To aggregate across the results is a little bit tricky because we have to extract the coefficients and standard errors from each model, format them in a particular way, and then feed that structure into the mi.meld function:
coefs.amelia <- do.call(rbind, lapply(lm.amelia.out, function(i) coef(summary(i))[,
1]))
ses.amelia <- do.call(rbind, lapply(lm.amelia.out, function(i) coef(summary(i))[,
2]))
mi.meld(coefs.amelia, ses.amelia)
## $q.mi
## (Intercept) x1 x2
## [1,] 0.09683 0.926 0.7436
##
## $se.mi
## (Intercept) x1 x2
## [1,] 0.2291 0.08222 0.1267
Now let's compare these results to those of our original model:
t(do.call(rbind, mi.meld(coefs.amelia, ses.amelia)))
## [,1] [,2]
## (Intercept) 0.09683 0.22908
## x1 0.92598 0.08222
## x2 0.74359 0.12674
coef(summary(lm))[, 1:2] # original results
## Estimate Std. Error
## (Intercept) 0.02587 0.21957
## x1 0.94835 0.08243
## x2 0.74874 0.12026
library(mi)
Let's start by visualizing the missing data:
mp.plot(mydf)
We can then see some summary information about the dataset and the nature of the missingness:
mi.info(mydf)
## names include order number.mis all.mis type collinear
## 1 x1 Yes 1 20 No positive-continuous No
## 2 x2 Yes 2 10 No continuous No
## 3 y Yes NA 0 No continuous No
With that information confirmed, it is incredibly issue to conduct our multiple imputation using the mi function:
imp.mi <- mi(mydf)
## Beginning Multiple Imputation ( Wed Nov 13 22:07:34 2013 ):
## Iteration 1
## Chain 1 : x1* x2*
## Chain 2 : x1* x2*
## Chain 3 : x1* x2*
## Iteration 2
## Chain 1 : x1* x2
## Chain 2 : x1* x2*
## Chain 3 : x1* x2*
## Iteration 3
## Chain 1 : x1* x2*
## Chain 2 : x1* x2
## Chain 3 : x1* x2
## Iteration 4
## Chain 1 : x1 x2
## Chain 2 : x1* x2
## Chain 3 : x1 x2
## Iteration 5
## Chain 1 : x1* x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2*
## Iteration 6
## Chain 1 : x1* x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2*
## Iteration 7
## Chain 1 : x1 x2*
## Chain 2 : x1* x2
## Chain 3 : x1 x2
## Iteration 8
## Chain 1 : x1 x2
## Chain 2 : x1* x2
## Chain 3 : x1* x2
## Iteration 9
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1* x2
## Iteration 10
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 11
## Chain 1 : x1* x2*
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 12
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 13
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 14
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## mi converged ( Wed Nov 13 22:07:37 2013 )
## Run 20 more iterations to mitigate the influence of the noise...
## Beginning Multiple Imputation ( Wed Nov 13 22:07:37 2013 ):
## Iteration 1
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 2
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 3
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 4
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 5
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 6
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 7
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 8
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 9
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 10
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 11
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 12
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 13
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 14
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 15
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 16
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 17
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 18
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 19
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## Iteration 20
## Chain 1 : x1 x2
## Chain 2 : x1 x2
## Chain 3 : x1 x2
## mi converged ( Wed Nov 13 22:07:41 2013 )
imp.mi
##
## Multiply imputed data set
##
## Call:
## .local(object = object, n.iter = ..3, R.hat = ..4, max.minutes = ..2,
## run.past.convergence = TRUE)
##
## Number of multiple imputations: 3
##
## Number and proportion of missing data per column:
## names type number.mis proportion
## 1 x1 positive-continuous 20 0.2
## 2 x2 continuous 10 0.1
## 3 y continuous 0 0.0
##
## Total Cases: 100
## Missing at least one item: 2
## Complete cases: 72
The results above report how many imputed datasets were produced and summarizes some of the results we saw above. For linear regression (and several other common models), the mi package includes functions that automatically run the model on each imputed dataset and aggregate the results:
lm.mi.out <- lm.mi(y ~ x1 + x2, imp.mi)
We can extract the results using the following:
coef.mi <- [email protected]
# or see them quickly with:
display(lm.mi.out)
## =======================================
## Separate Estimates for each Imputation
## =======================================
##
## ** Chain 1 **
## lm(formula = formula, data = mi.data[[i]])
## coef.est coef.se
## (Intercept) -0.01 0.18
## x1 0.96 0.06
## x2 0.75 0.09
## ---
## n = 100, k = 3
## residual sd = 0.89, R-Squared = 0.75
##
## ** Chain 2 **
## lm(formula = formula, data = mi.data[[i]])
## coef.est coef.se
## (Intercept) 0.03 0.18
## x1 0.94 0.07
## x2 0.67 0.09
## ---
## n = 100, k = 3
## residual sd = 0.91, R-Squared = 0.74
##
## ** Chain 3 **
## lm(formula = formula, data = mi.data[[i]])
## coef.est coef.se
## (Intercept) -0.02 0.20
## x1 0.96 0.07
## x2 0.69 0.10
## ---
## n = 100, k = 3
## residual sd = 0.96, R-Squared = 0.71
##
## =======================================
## Pooled Estimates
## =======================================
## lm.mi(formula = y ~ x1 + x2, mi.object = imp.mi)
## coef.est coef.se
## (Intercept) 0.00 0.19
## x1 0.95 0.07
## x2 0.71 0.10
## ---
Let's compare these results to our original model:
do.call(cbind, coef.mi) # multiply imputed results
## coefficients se
## (Intercept) 0.00123 0.18878
## x1 0.95311 0.06901
## x2 0.70687 0.10411
coef(summary(lm))[, 1:2] # original results
## Estimate Std. Error
## (Intercept) 0.02587 0.21957
## x1 0.94835 0.08243
## x2 0.74874 0.12026
library(mice)
To conduct the multiple imputation, we simply need to run the mice function:
imp.mice <- mice(mydf)
##
## iter imp variable
## 1 1 x1 x2
## 1 2 x1 x2
## 1 3 x1 x2
## 1 4 x1 x2
## 1 5 x1 x2
## 2 1 x1 x2
## 2 2 x1 x2
## 2 3 x1 x2
## 2 4 x1 x2
## 2 5 x1 x2
## 3 1 x1 x2
## 3 2 x1 x2
## 3 3 x1 x2
## 3 4 x1 x2
## 3 5 x1 x2
## 4 1 x1 x2
## 4 2 x1 x2
## 4 3 x1 x2
## 4 4 x1 x2
## 4 5 x1 x2
## 5 1 x1 x2
## 5 2 x1 x2
## 5 3 x1 x2
## 5 4 x1 x2
## 5 5 x1 x2
We can see some summary information about the imputation process:
summary(imp.mice)
## Multiply imputed data set
## Call:
## mice(data = mydf)
## Number of multiple imputations: 5
## Missing cells per column:
## x1 x2 y
## 20 10 0
## Imputation methods:
## x1 x2 y
## "pmm" "pmm" ""
## VisitSequence:
## x1 x2
## 1 2
## PredictorMatrix:
## x1 x2 y
## x1 0 1 1
## x2 1 0 1
## y 0 0 0
## Random generator seed value: NA
To run our regression we use the lm function wrapped in a with call, which estimates our model on each imputed dataframe:
lm.mice.out <- with(imp.mice, lm(y ~ x1 + x2))
summary(lm.mice.out)
##
## ## summary of imputation 1 :
##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7893 -0.7488 0.0955 0.7205 2.3768
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0636 0.1815 0.35 0.73
## x1 0.9528 0.0660 14.44 < 2e-16 ***
## x2 0.6548 0.0916 7.15 1.6e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.915 on 97 degrees of freedom
## Multiple R-squared: 0.735, Adjusted R-squared: 0.73
## F-statistic: 135 on 2 and 97 DF, p-value: <2e-16
##
##
## ## summary of imputation 2 :
##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4456 -0.6911 0.0203 0.6839 2.6271
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1618 0.1855 0.87 0.39
## x1 0.8849 0.0671 13.19 < 2e-16 ***
## x2 0.7424 0.0942 7.88 4.7e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.939 on 97 degrees of freedom
## Multiple R-squared: 0.721, Adjusted R-squared: 0.716
## F-statistic: 126 on 2 and 97 DF, p-value: <2e-16
##
##
## ## summary of imputation 3 :
##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5123 -0.6683 0.0049 0.6717 2.5072
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0800 0.1872 0.43 0.67
## x1 0.9402 0.0674 13.95 < 2e-16 ***
## x2 0.8150 0.0947 8.61 1.3e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.94 on 97 degrees of freedom
## Multiple R-squared: 0.721, Adjusted R-squared: 0.715
## F-statistic: 125 on 2 and 97 DF, p-value: <2e-16
##
##
## ## summary of imputation 4 :
##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5739 -0.7310 0.0152 0.6534 2.4748
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0787 0.1815 0.43 0.67
## x1 0.9438 0.0660 14.30 < 2e-16 ***
## x2 0.7833 0.0916 8.55 1.8e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.922 on 97 degrees of freedom
## Multiple R-squared: 0.732, Adjusted R-squared: 0.726
## F-statistic: 132 on 2 and 97 DF, p-value: <2e-16
##
##
## ## summary of imputation 5 :
##
## Call:
## lm(formula = y ~ x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6130 -0.6085 0.0085 0.6907 2.4719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0477 0.1858 0.26 0.8
## x1 0.9463 0.0672 14.07 < 2e-16 ***
## x2 0.7436 0.0893 8.33 5.4e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.91 on 97 degrees of freedom
## Multiple R-squared: 0.738, Adjusted R-squared: 0.733
## F-statistic: 137 on 2 and 97 DF, p-value: <2e-16
The results above are for each separate dataset. But, to pool them, we use pool:
pool.mice <- pool(lm.mice.out)
Let's compare these results to our original model:
summary(pool.mice) # multiply imputed results
## est se t df Pr(>|t|) lo 95 hi 95 nmis
## (Intercept) 0.08637 0.19056 0.4533 81.41 6.516e-01 -0.2927 0.4655 NA
## x1 0.93361 0.07327 12.7422 50.19 0.000e+00 0.7865 1.0808 20
## x2 0.74782 0.11340 6.5944 22.53 1.104e-06 0.5130 0.9827 10
## fmi lambda
## (Intercept) 0.08664 0.06447
## x1 0.20145 0.17025
## x2 0.38954 0.33765
coef(summary(lm))[, 1:2] # original results
## Estimate Std. Error
## (Intercept) 0.02587 0.21957
## x1 0.94835 0.08243
## x2 0.74874 0.12026
It is useful at this point to compare the coefficients from each of our multiple imputation methods. To do so, we'll pull out the coefficients from each of the three packages' results, our original observed results (with case deletion), and the results for the real data-generating process (before we introduced missingness). Amelia package results
s.amelia <- t(do.call(rbind, mi.meld(coefs.amelia, ses.amelia)))
mi package results
s.mi <- do.call(cbind, coef.mi) # multiply imputed results
mice package results
s.mice <- summary(pool.mice)[, 1:2] # multiply imputed results
Original results (case deletion)
s.orig <- coef(summary(lm))[, 1:2] # original results
Real results (before missingness was introduced)
s.real <- summary(lm(y ~ x1 + x2))$coef[, 1:2]
Let's print the coefficients together to compare them:
allout <- cbind(s.real[, 1], s.amelia[, 1], s.mi[, 1], s.mice[, 1], s.orig[,
1])
colnames(allout) <- c("Real Relationship", "Amelia", "MI", "mice", "Original")
allout
## Real Relationship Amelia MI mice Original
## (Intercept) 0.04502 0.09683 0.00123 0.08637 0.02587
## x1 0.95317 0.92598 0.95311 0.93361 0.94835
## x2 0.82900 0.74359 0.70687 0.74782 0.74874
All three of the multiple imputation models - despite vast differences in underlying approaches to imputation in the three packages - yield strikingly similar inference. This was a relatively basic and all of the packages offer a number of options for more complicated situations than what we examined here. While executing multiple imputation requires choosing a package and typing some potentially tedious code, the results are almost always going to be better than doing the easier thing of deleting cases and ignoring the consequences thereof.