Missing data and imputation

Introduction

The aim of this practical is to build understanding of missing data, i.e., how different missingness mechanisms affect your estimates, why naive fixes like mean imputation and listwise deletion don’t solve the problem, and how proper (multiple) imputation does.

In Part 1, we manually generate missingness under three mechanisms in an income dataset, and show why mean imputation fails. In Part 2, we move to a new dataset and work through a sequence of increasingly better imputation approaches, ending with multiple imputation using the mice package, and compare how each one affects a regression estimate.

library(MASS)
library(tidyverse)
library(mice)
library(ggmice)
library(ISLR)

Make sure to load MASS before tidyverse, otherwise the function MASS::select() will overwrite dplyr::select().

Part 1: Missingness mechanisms

Before we start, set a seed for reproducibility, we use 45.

set.seed(45)

Income dataset

In this part, we will manually create missingness to get an idea how different missingness mechanisms affect the outcomes of analyses. To this end, we start with a complete dataset, income.rds which can be downloaded here. Load the data in R after downloading as follows:

income <- readRDS("data/income.rds")

Generating missing values

The general expression of the missing data model is \(Pr(R = 0|Y_{\text{obs}}, Y_{\text{mis}}, \psi)\). Let’s go through this formula step by step.

• \(R\) is a matrix with response indicators, which contains the locations of the missing values. For each person and each variable, \(R_{i,j}\) indicates whether person \(i\) has a missing value (\(R_{i,j} = 0\)) or an observed value (\(R_{i,j} = 1\)) on variable \(j\). In our case, we have \(R_{\text{age}}\), \(R_{\text{gender}}\) and \(R_{\text{income}}\).

• \(Y_{\text{obs}}\) denotes the observed data.

• \(Y_{\text{mis}}\) denotes the missing data. Note, we don’t have these data, but they might influence whether or not a person has missings on a variable.

• \(\psi\) denotes the parameters for the missing data models. These parameters relate the observed and missing values to the response indicator.

Now we have a binary response indicator, and some additional variables, we can build a classification model to predict whether or not someone has a missing value on some variable. Normally, we would use the data to evaluate whether the missingness is related to someone’s values on other variables. However, we can also reverse the process, and impose a missingness structure on the observed data that depends on the variables.

In this practical, we will create missingness in the variable income, and we will do that using a logistic regression model. The logistic regression model predicts the probability of a missing value, given a set of indicators.

Note that a logistic regression model for the missingness can be defined as \[ \Pr(R_{\text{income}} = 0 | Y_{\text{obs}}, Y_{\text{mis}}, \psi) = \frac{\exp\{\psi_0 + \psi_1 * \text{age} + \psi_2 * (\text{gender}=\text{male}) + \psi_3 * \text{income}\}}{1 + \exp\{\psi_0 + \psi_1 * \text{age} + \psi_2 * (\text{gender} = \text{male}) + \psi_3 * \text{income}\}}. \] In this equation, \(\psi_0\) is a parameter for the baseline probability of missingness in \(\text{income}\), \(\psi_1\) and \(\psi_2\) relate the observed data to the missingness in \(\text{income}\), and \(\psi_3\) relates unobserved information to the missingness in \(\text{income}\).

If we would translate this equation into R-code, we would have

prob_mis <- function(psi0, psi1, psi2, psi3, age, gender, income) {
  exp(psi0 + psi1 * age + psi2 * (gender == "Male") + psi3 * income) /
    (1 + exp(psi0 + psi1 * age + psi2 * (gender == "Male") + psi3 * income))
}

Not Data-Dependent (NDD)

The data are considered to be NDD if \(Pr(R = 0|Y_{obs}, Y_{mis}, \psi) = Pr(R=0|\psi)\), indicating that the missingness probability is unrelated to the data \(Y\) and only depends on some parameters \(\psi\).

1. Use the function above to add a variable with the missingness probabilities misprob_ndd to the income data, in which you specify all missingness parameters to be 0 (i.e., \(\psi_0 = \psi_1 = \psi_2 = \psi_3\) = 0), such that the probability of having a missing value equals \(0.5\) for everyone.

Tip: use the mutate() function to create the new variable.

Using these probabilities, we can randomly draw the observations that will have a missing value on \(\text{income}\) from a binomial distribution.

2. Use rbinom() to create a missingness indicator R_ndd using the misprob_ndd probabilities, and add this indicator to the income data.

NB: The function prob_mis gives the probability of having a missing, while the missingness indicator must give zero if income is not observed, so use 1 - rbinom() to create a proper missingness indicator.

3. Use the function ifelse() to add a new variable income_ndd to the income data, that equals NA if R_ndd equals 0, and the observed income value otherwise.

4. Verify that approximately \(50\%\) of the observations have a missing value on income_ndd.

Seen Data-Dependent (SDD)

So far, there was a baseline probability to have missing values, but this was the same for all observations (hence the term, Not Data-Dependent). Data are considered to be Seen Data-Dependent if the probability of having a missing value depends on the observed variables. Formally, we can define this as \(Pr(R = 0|Y_{obs}, Y_{mis}, \psi) = Pr(R=0|Y_{obs},\psi)\).

In this example, we let the missingness in \(\text{income}\) depend on the observed values of \(\text{age}\) and \(\text{gender}\), in such a way that older people have a higher probability of having a missing, and that males have a higher probability of having a missing than females.

5. Use the prob_mis() function to create misprob_sdd with psi0 = -3, psi1 = 0.05, psi2 = 0.5 and psi3 = 0. Then, as before, use rbinom() to create a missingness indicator R_sdd, and ifelse() to create income_sdd (NA where R_sdd equals 0).

6. Create a plot that shows how the probability of having a missing value (on the y-axis) depends on \(\text{age}\) (on the x-axis) and \(\text{gender}\) (as color aesthetic).

Unseen Data-Dependent (UDD)

We now take the missingness one step further again, by making it dependent on unobserved data. In practice, it is often the case that those who have a higher income are less willing to share their income, which is an example of unseen data-dependent missingness. Namely, we cannot observe whether those who have a missing on income have a higher income, because this is exactly the information that we miss. The missingness could also depend on variables that we have not observed, for example because people who are self-employed might be less willing to answer survey questions about their income than people who have an employer. In both situations, we call the missingness Unseen Data-Dependent. Formally, this is defined as \(Pr(R = 0|Y_{obs}, Y_{mis}, \psi) = Pr(R=0|Y_{obs},Y_{mis},\psi)\).

Because we have no unobserved variables in this example, we let the missingness in \(\text{income}\) depend on the income values themselves, such that higher incomes have a higher probability of being missing.

7. Use the prob_mis() function to create misprob_udd with psi0 = -7, psi1 = 0, psi2 = 0 and psi3 = 0.003. Then, as before, use rbinom() to create a missingness indicator R_udd, and ifelse() to create income_udd (NA where R_udd equals 0).

8. Create a plot that shows how the probability of having a missing value misprob_udd (on the y-axis) depends on income (on the x-axis).

Note: In practice, we would never be able to do this, because we would not have the actual incomes for those who refuse to provide their income values.

Visualizing missing data

9. Create a boxplot of the observed values and missing values of income under each of the three missing data models. What do you notice?

Mean imputation

We will now show why mean imputation is, in general, a bad idea.

10. Use mutate() in combination with ifelse() and mean(x, na.rm = TRUE) to create imp_income_ndd, imp_income_sdd and imp_income_udd, in which missings on income_ndd, income_sdd and income_udd are replaced with the means of the respective variables. Then compare the means and standard deviations of the complete, observed income variable with each of the imputed variables. What differences do you see?

11. Create a plot in which you visualize the density or histogram of observed income, imp_income_ndd, imp_income_sdd, imp_income_udd. What do you see? Do you think it is a good idea to do mean imputation?

Mean imputation clearly distorts the distribution. In Part 2, we work through better alternatives, ending with multiple imputation.

Part 2: Ad hoc and multiple imputation with mice

Before we start, set a seed for reproducibility, we use 123, and also use options(scipen = 999) to suppress scientific notations, making it easier to compare and interpret results later in the session.

set.seed(123)
options(scipen = 999)

Data processing

For this part of the practical, you can directly copy and paste the code as given.

From the ISLR package, we will use the College dataset. Check out help(College) to find out more about this dataset. We will only use the following four variables:

• Outstate: Out-of-state tuition
• PhD: Pct. of faculty with Ph.D.’s
• Terminal: Pct. of faculty with terminal degree
• Expend: Instructional expenditure per student

Hence, we create a new data frame containing only these four variables:

data_complete <- 
  College |> 
  select(Outstate, PhD, Terminal, Expend) 

summary(data_complete)

    Outstate          PhD            Terminal         Expend     
 Min.   : 2340   Min.   :  8.00   Min.   : 24.0   Min.   : 3186  
 1st Qu.: 7320   1st Qu.: 62.00   1st Qu.: 71.0   1st Qu.: 6751  
 Median : 9990   Median : 75.00   Median : 82.0   Median : 8377  
 Mean   :10441   Mean   : 72.66   Mean   : 79.7   Mean   : 9660  
 3rd Qu.:12925   3rd Qu.: 85.00   3rd Qu.: 92.0   3rd Qu.:10830  
 Max.   :21700   Max.   :103.00   Max.   :100.0   Max.   :56233  

During the practical, we focus on how different imputation approaches affect the linear relationship between Outstate and Expend. Therefore, we first inspect how this relationship looks like on the complete dataset.

lm_complete <- lm(Expend ~ Outstate, 
                  data = data_complete)

summary(lm_complete)

Call:
lm(formula = Expend ~ Outstate, data = data_complete)

Residuals:
   Min     1Q Median     3Q    Max 
 -6131  -2022   -637   1027  39273 

Coefficients:
             Estimate Std. Error t value            Pr(>|t|)    
(Intercept) 542.86969  385.92915   1.407                0.16    
Outstate      0.87325    0.03449  25.315 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3866 on 775 degrees of freedom
Multiple R-squared:  0.4526,    Adjusted R-squared:  0.4519 
F-statistic: 640.9 on 1 and 775 DF,  p-value: < 0.00000000000000022

We now generate missingness in the variables Expend and Terminal under a Seen Data-Dependent (SDD) mechanism, the same mechanism you generated by hand in Part 1, assuming that the missingness in these variables are dependent on the values of Outstate and PhD respectively.

logistic <- function(x) (exp(x) / (1 + exp(x))) # logistic function

N <- nrow(data_complete) # number of observations

misprob_out <- logistic(scale(data_complete$Outstate)) # predictor Outstate
misprob_phd <- logistic(scale(data_complete$PhD)) # predictor PhD

data_missing <- 
  data_complete |> 
  mutate(
    R_out = 1 - rbinom(N, 1, misprob_out), # reponse indicator outstate
    R_phd = 1 - rbinom(N, 1, misprob_phd), # response indicator PhD
    m_Expend = ifelse(R_out == 0, NA, Expend),    # Missing in expend
    m_Terminal = ifelse(R_phd == 0, NA, Terminal) # and in terminal
  ) |> 
  select(Outstate, m_Expend, m_Terminal) # Select variables that we will use

Inspection of missingness

12. Inspect the missingness pattern of the dataset using the function plot_pattern() from the ggmice package.

Ad hoc imputation methods

Listwise deletion

13. Create the object lm_listwise by fitting a regression model with m_Expend as the outcome, predicted by Outstate, using the data containing missing values. Then inspect it with summary(). How many observations are excluded? Is there a significant effect of Outstate?

Hint: By default, lm() removes the observations that have missings on one of the variables in the model.

Regression imputation

14. Use mice() to impute data_missing with regression imputation (method = "norm.predict"), using a single imputation (m = 1) and a single iteration (maxit = 1), name the output imp_regression, then extract the complete data with complete() and confirm there are no more missing values.

NB: m = 1 implies that only a single imputed data set is created, whereas maxit = 1 implies that mice() does not iterate to converge to a distribution.

The with(data, expr) function allows to fit lm() and glm() models with imputed data sets (mids-objects, to be precise), without requiring any further data handling. Check how the function works by running ?with.mids().

15. Use with() to fit lm(m_Expend ~ Outstate) on imp_regression, assign it to lm_regression, and inspect the results. What happens to the regression coefficient and the standard error compared to listwise deletion?

16. Use ggmice() to create a scatterplot of m_Expend (y) versus Outstate (x) for the regression-imputed data. What do you notice?

Hint: ggmice() works much the same way as ggplot(), but accepts mids objects. Check ?ggmice for help.

Multiple imputation

We skip straight to the proper solution: multiple imputation, which, unlike the ad hoc methods above, accounts for the uncertainty in both the imputed values themselves and in the regression coefficients used to generate them.

17. Impute the missing data by calling mice() with default settings, assign the result to imp_multiple, and inspect the output object. What information does it give you?

18. Perform the linear regression on the imputed data using with() and assign it to lm_multiple. What do you notice when you inspect this object?

19. Use pool() to combine the regression analyses into pool_multiple. Inspect pool_multiple and summary(pool_multiple). What information do these objects provide?

Unlike regression imputation, multiple imputation carries the uncertainty from both the imputed values and the estimated coefficients through into the pooled standard error, which is exactly what the lambda and fmi figures above quantify.

20. Now compare the results from the different methods applied (complete data, listwise deletion, regression imputation, multiple imputation). What is your conclusion?