Imputing Missing Data

Listwise deletion

We will only use rows for which all the variables used in the model are non-missing. This is what lm() does by default when it encounters missing values.

mod_listwise <-
  lm(infantMortality ~ log(GDPperCapita) + contraception + educationFemale,
     data = UN)
summary(mod_listwise)

## 
## Call:
## lm(formula = infantMortality ~ log(GDPperCapita) + contraception + 
##     educationFemale, data = UN)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.251 -11.145   0.634   7.925  66.125 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        148.245     11.913  12.444   <2e-16 ***
## log(GDPperCapita)   -7.580      2.365  -3.205   0.0022 ** 
## contraception      -46.114     17.399  -2.650   0.0103 *  
## educationFemale     -2.491      1.386  -1.798   0.0775 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.9 on 58 degrees of freedom
##   (145 observations deleted due to missingness)
## Multiple R-squared:  0.7466, Adjusted R-squared:  0.7335 
## F-statistic: 56.97 on 3 and 58 DF,  p-value: < 2.2e-16

Indicator Variables for Missing Values

For each variable, we will add a dummy variable that is equal 1 if that variable is missing and 0 if not. Then for each variable we will replace its missing values with 0 (or its mean). To assist in that, we will define the function fill_na which fills in missing values in x with whatever value is given in fill.

fill_na <- function(x, fill = 0) {
  x[is.na(x)] <- fill
  x
}

Now we generate a new dataset with those dummy variables and with the missing values replaced by their means.

UN_with_dummies <-
  UN %>%
  mutate(GDPperCapita_na = as.integer(is.na(GDPperCapita)),
         GDPperCapita = fill_na(GDPperCapita, mean(GDPperCapita, na.rm = TRUE)),
         contraception_na = as.integer(is.na(contraception)),
         contraception = fill_na(contraception, mean(contraception, na.rm = TRUE)),
         educationFemale_na = as.integer(is.na(educationFemale)),
         educationFemale = fill_na(educationFemale, mean(educationFemale, na.rm = TRUE)))
mod_dummies <- lm(infantMortality ~ log(GDPperCapita) + GDPperCapita_na + contraception +
     contraception_na + educationFemale + educationFemale, data = UN_with_dummies)
mod_dummies

## 
## Call:
## lm(formula = infantMortality ~ log(GDPperCapita) + GDPperCapita_na + 
##     contraception + contraception_na + educationFemale + educationFemale, 
##     data = UN_with_dummies)
## 
## Coefficients:
##       (Intercept)  log(GDPperCapita)    GDPperCapita_na  
##          166.4830           -12.5506             3.7953  
##     contraception   contraception_na    educationFemale  
##          -55.0697             3.5130            -0.4682

Unconditional Mean Imputation

Another ad hoc method is to replace missing values with their means. Since we filled in the missing values in UN_with_dummies with the mean values, we will use that, but not include the dummy variables.

These bivariate plots illustrate how unconditional imputation means work, and give an indication of how they will affect results.

ggplot(filter(UN_with_dummies, !is.na(infantMortality),
              !is.na(educationFemale)),
       aes(y = infantMortality, x = educationFemale,
           colour = as.factor(educationFemale_na))) +
  geom_point() +
  scale_colour_discrete("Imputed") +
  theme_local()

ggplot(na.omit(select(UN_with_dummies, infantMortality,
                      contraception, contraception_na)),
       aes(y = infantMortality, x = contraception,
           colour = as.factor(contraception_na))) +
  geom_point() +
  scale_colour_discrete("Imputed") +
  theme_local()

mod_means <- lm(infantMortality ~ log(GDPperCapita) + contraception +
      educationFemale, data = UN_with_dummies)
mod_means

## 
## Call:
## lm(formula = infantMortality ~ log(GDPperCapita) + contraception + 
##     educationFemale, data = UN_with_dummies)
## 
## Coefficients:
##       (Intercept)  log(GDPperCapita)      contraception  
##          166.5511           -12.2072           -56.1593  
##   educationFemale  
##           -0.5595

Multiple Imputation

Our preferred method is multiple imputation. For multiple imputation we will use the implementation in Amelia.¹ Multiple imputation works by creates \(m\) datasets with the missing values filled in with imputations. The analyst then runs the analysis on each of the datasets, and there are (fairly easy) methods to combine the estimates and standard errors from the individual analyses into an overall point estimate and standard error. The function amelia creates these imputations.

• m = 5 is the number of observations. Usually 5 - 10 are sufficient. See Fox, p. 564.
• amelia models the distribution of the data as multivariate normal so it is useful to transform the data to make it more normal. Use the arguments logs to indicate variables to log transform, and lgstc to logistic transform variables.
• Additionally, you need to tell amelia which variables are nominal (noms =) so it creates dummy variables; and which variables are id variables (idvars =) so it can ignore them in the imputation process, but still include them in the final datasets.

UN_mi <-
  amelia(data.frame(UN),
         m = 5,
         logs = c("GDPperCapita"),
         lgstc = c("economicActivityMale", "economicActivityFemale",
                   "illiteracyMale", "illiteracyFemale", "contraception"),
         noms = c("region"),
         idvars = "country")

Note I had to use data.frame because at the moment amelia has bugs when using dplyr data_frame objects.

Note several things about the imputation process.

First, the number of imputations is small, m = 5. There is not much improvement in the performance of multiple imputation from using many imputations. Usually 5–10 is more than sufficient. See Fox, p. 564.

Second, we use the dependent variable as well as variables not used in the regression in the imputation. (Fox, p. 567) This is not an issue. Concerns like endogeneity and multicollinearity are not concerns in this context; we are not interested in estimating unbiased coefficients. We are concerned with predicting the values of these variables.

The multiple imputation process has created m = 5 separate datasets, which are stored as a list of data frames in the imputations element of the object:

str(UN_mi$imputations)

Each of these datasets is complete, meaning it has no missing observations. But for values which were missing in the original data, they are filled with different imputations. For example, compare this subset in the orginal data and two imputations.

countries <- c("Afghanistan", "Angola", "Algeria", "Antigua")
UN %>% filter(country %in% countries) %>% select(country, contraception)

## Source: local data frame [4 x 2]
## 
##       country contraception
##         (chr)         (dbl)
## 1 Afghanistan            NA
## 2     Algeria          0.52
## 3      Angola            NA
## 4     Antigua          0.53

UN_mi$imputations[[1]] %>% filter(country %in% countries) %>% select(country, contraception)

##       country contraception
## 1 Afghanistan    0.07640913
## 2     Algeria    0.52000000
## 3      Angola    0.14690295
## 4     Antigua    0.53000000

UN_mi$imputations[[2]] %>% filter(country %in% countries) %>% select(country, contraception)

##       country contraception
## 1 Afghanistan    0.06033049
## 2     Algeria    0.52000000
## 3      Angola    0.28073275
## 4     Antigua    0.53000000

The following creates plots of each variable, overplotting the values from each imputation against the observation number. Points which are black were observed in the original data; points which are gray are imputed.

all_imputations <-
  data_frame(i = seq_len(UN_mi$m)) %>%
    rowwise() %>%
    do({
      mutate(UN_mi$imputations[[.$i]],
             .obsnum = row_number(),
             .imputation = .$i)
    })

imputation_plot <- function(yvar) {
  ggplot(all_imputations, aes_string(x = ".obsnum",
                                     y = yvar)) +
  geom_point(alpha = 0.2) +
  theme_minimal() +
  xlab("Observation number") +
  theme_local()
}

imputation_plot("educationFemale")

imputation_plot("contraception")

imputation_plot("log(GDPperCapita)")

Let’s see how one of these imputed datasets looks and compare the imputed and non-imputed values as an informal way of judging its plausibility.

ggplot(UN_mi$imputations[[1]] %>%
         mutate(missing = UN_mi$missMatrix[ , "educationFemale"]),
       aes(x = educationFemale, y = infantMortality, colour = missing)) +
  geom_point() +
  theme_local()

ggplot(UN_mi$imputations[[1]] %>%
         mutate(missing = UN_mi$missMatrix[ , "contraception"]),
       aes(x = contraception, y = infantMortality, colour = missing)) +
  geom_point() +
  theme_local()

These look much more plausible than those of the unconditionally imputed mean.

The plots we just made suggest Amelia is working well, but let’s look at a couple of diagnostics that Amelia suggests: The first is to drop each observation of a variable, impute it, and plot the values against the original values. If the imputed values resemble the original values, then the imputation procedure has worked well.

Amelia suggests several diagnostics for evaluating The first is “overimputation”. This is similar to leave-one-out cross validation. Each observation is removed and treated as a missing value to be imputated. Amelia plots 90% confidence intervals for each observation against the true values. If the imputation is reasonable, then these confidence intervals should have 90% coverage of the true value. The following code runs the overimputation diagnostic for the variables of interest (the ones that will be used in the regression). The overimputation diagnostics do not suggest any issues with the imputation of these variables.

for (var in c("GDPperCapita", "contraception", "illiteracyFemale")) {
  overimpute(UN_mi, var, main = var)
}

The second method is to plot the marginal density of the imputed values against observed values. If these appear too different, then there may be issues with imputation. However, there is no particular criteria for what is too different, and there may be reasons why these distributions diverge; see the Amelia vignette for an example of what would be a big difference. But if these distributions are different, then you should revisit the data to see whether imputation makes sense. The following code performs this for each variable of interest. The plots do not suggest any differences to be concerned about.

for (var in c("GDPperCapita", "contraception", "illiteracyFemale")) {
  compare.density(UN_mi, var, main = var)
}

Now that we are satisfied that the imputations are reasonable, we run the analysis on each of the imputated datasets. The following code runs the regression on each imputed dataset saving the coefficients to b_out, and the standard errors of the coefficients to se_out. These are originally saved to a list, but rbind converts it to a matrix with each imputations as rows and coefficients as columns.

Now estimate the model using each imputed dataset. We save the results to b_out and se_out.

b_out <- list()
se_out <- list()
for (i in seq_along(UN_mi$imputations)) {
  mod <- lm(infantMortality ~ log(GDPperCapita) + contraception + educationFemale, data = UN_mi$imputations[[i]])
  b_out[[i]] <- coef(mod)
  se_out[[i]] <- sqrt(diag(vcov(mod)))
}
b_out <- do.call(rbind, b_out)
se_out <- do.call(rbind, se_out)

As described in Fox or Gelman, the coefficients from the imputations can be combined to calculate a single point estimate and standard error. The function mi.meld does this, and returns a list with point estimates in the element q.mi and standard errors in se.mi:

mod_mi_res <- mi.meld(q = b_out, se = se_out)
mod_mi_res

## $q.mi
##      (Intercept) log(GDPperCapita) contraception educationFemale
## [1,]    134.5969         -2.699894     -45.31222       -5.023085
## 
## $se.mi
##      (Intercept) log(GDPperCapita) contraception educationFemale
## [1,]    13.43238          3.431773      14.32526        1.894144

Model Comparison

Now that we’ve estimated this model using various methods to handle missing data, let’s compare the results. First, we combine these into a single dataset. For the results returned by Amelia, since there is no tidy function defined for Amelia objects we need to manually create a data frame consistent with those returned by tidy.

model_comp <-
  bind_rows(tidy(mod_listwise) %>%
              mutate(model = "listwise"),
            tidy(mod_dummies) %>%
              mutate(model = "dummies"),
            tidy(mod_means) %>%
              mutate(model = "means"),
            data_frame(term = colnames(mod_mi_res$q.mi),
                       estimate = as.numeric(mod_mi_res$q.mi),
                       std.error = as.numeric(mod_mi_res$se.mi),
                       model = "mi"))

plot_mi_models <- function(data, variable) {
  ggplot(data %>%
         filter(term %in% variable),
         aes(x = model, y = estimate,
             ymin = estimate - 2 * std.error,
             ymax = estimate + 2 * std.error)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, colour = "red", alpha = 0.2) +
  coord_flip() +
  theme_minimal()
}

plot_mi_models(model_comp, "log(GDPperCapita)")

plot_mi_models(model_comp, "contraception")

plot_mi_models(model_comp, "educationFemale")

These plots show that the various methods produce different answers both in terms of point-estimates and standard errors. Though, we would have to use simulations to determine why the imputation methods are better than the alternatives.