A Bayesian Perspective on Missing Data Imputation

Missing data is a quite prevalent issue among many applied data sets. This lecture is a light introduction to one way we can perform missing data imputation from a Bayesian perspective. For more information, the book Statistical Analysis with Missing Data by Roderick Little and Donald Rubin is a great resource.

Diabetes example

The motivating example comes from Hoff’s book, page 115. Consider two measurements on 200 Arizona women with Pima Indian heritage:

• skin, the skin fold thickness, and
• bmi, the body mass index.

Our goal is to study the relationship between these two variables, accounting for the missing data present in the table.

1
2
3

library(tidyverse)
library(ggpubr)
library(rjags)

1
2
3
4
5
6
7
8
9

diab <- data.table::fread(
  "http://www2.stat.duke.edu/~pdh10/FCBS/Exercises/diabetes_200_miss.dat",
  header = T
) %>%
  select(skin, bmi)

diab %>%
  head() %>%
  knitr::kable()

Exploratory data analysis

As shown in the table above, there are obviously missing values coded as NA’s in our data set¹.

We can use the summary() function in R to get the 5-number summary, mean, and number of missing values for the columns.

Histograms² are helpful in determining the distributions of the variables. Both are roughly symmetric, and there’s one extreme observation for skin fold thickness.

Figure 1: Observed BMI measurements and skin fold thicknesses.

We can also use a scatter plot to show the relationship between the two variables.

We observe a relatively strong positive relationship except for the one outlier, and we might want to quantify this in terms of how correlated the two variables are.

Finally, when we have missing values in our data set, it can be helpful to visualize the pattern of missingness. The aggr() function from the VIM R package shows proportions of missing values for each variable, as well as how frequently missingness occurrs for various variable combinations.

1

VIM::aggr(diab)

There’s about 13% missingness for the skin variable, and about 11% for the bmi variable. The blue color on the right indicates presence of data, whereas the red indicates missing values. We have both values for the large majority of the data, but in a few entries we are missing one or both of the values.

Different strategies might be applied for different types of missingness. For example, if we always have values in one variable but have missing entries in another variable, or if the two variables are always missing together.

Bayesian models

Our two variables appear correlated with each other, so we specify a bivariate normal data distribution for the sampling model:

where \(Y_{i1}\) is the \(i^{th}\) person’s skin, and \(Y_{i2}\) is the \(i^{th}\) person’s BMI. The two parameters of the MVN distribution are its mean vector and covariance matrix. The multivariate normal is just an extension of the normal distribution that accounts for correlation between components of the random variables³.

For the prior model, we have essentially six parameters. \(\theta_1\) and \(\theta_2\) are means, so we specify a normal distribution with a large variance:

For the covariance matrix of the multivariate normal distribution, there’s a conjugate prior called the inverse Wishart distribution, which is a distribution over matrices that has two parameters.

Posterior distributions for our parameters are given by Bayes’ Theorem:

\begin{aligned} p(\theta, \Sigma \mid y) &\propto p(y \mid \theta, \Sigma) p(\theta, \Sigma) \\ &\propto \left[ \prod_{i=1}^n p(y_i \mid \theta, \Sigma) \right] p(\theta) p(\Sigma) \end{aligned}

A problem arises when we evaluate the sampling model likelihood, since this depends on all of the data, but we have missing values. For example, if we were missing \(Y_{31}\), then we can’t evaluate \(p(y_3 \mid \theta, \Sigma)\).

Simple approaches

There are many ways to handle missing data. The simplest approach is to perform a complete case analysis, where we delete entries with a missing value for at least one variable. If the data are missing completely at random (MCAR), then we get unbiased parameter estimate. However:

• Credible intervals will usually be too wide, since the sample size has effectively been reduced (drastically).
• If missing values occur for other reasons, this yields biased estimates.

The next step up is instead of dropping all the missing entries, we replace missing values of a feature with the mean or median of available values for the feature. This is called single mean/median imputation, and is often recommended in many machine learning tasks if the data are MCAR as it can still lead to unbiased parameter estimates.

There are also some drawbacks with this method, especially if you’re interested in studying the variability. Single mean/median imputation artificially reduces the variance of features, which results in credible intervals which are too narrow. As shown in Figure 2, the peak at the mean is much higher relative to the surrounding values. The mean of the distribution can still be estimated appropriately, but the amount of spread is reduced.

1
2
3
4
5
6
7

diab %>%
  select(bmi) %>%
  mutate(Imputed = bmi) %>%
  replace_na(list(Imputed = mean(.$Imputed, na.rm = T))) %>%
  select(`No imputation` = bmi, `With mean imputation` = Imputed) %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "BMI") %>%
  gghistogram(x = "BMI", bins = 10, facet.by = "Variable")

Figure 2: Mean imputation with the BMI feature.

Another drawback is it doesn’t account for relationship between variables, thus reduces correlation. For example, the correlation between the two variables is 0.665 if we only look at the complete cases. If we use single mean imputation, the correlation decreases to 0.589.

Note that if you don’t care about assessing variability in your estimates, which is often the case for classification tasks, then mean imputation can work quite well.

A Bayesian approach

A popular approach to address the issues of single mean imputation is called multiple imputation by chained equations (MICE), and proceeds as follows.

1. Specify separate imputation models for each feature with missing values, which are conditional on all other features. This allows us to exploit relationship between variables.
2. Use Gibbs sampling with the imputation models to generate multiple imputed data sets. This has the effect of allowing us to more appropriately assess uncertainty in the missing values.
3. Analyze each imputed data set separately, i.e. by sampling from the posterior using JAGS.
4. Combine the results across imputed data sets at the end.

Recall that in the Gibbs sampler we generate values at step \(s+1\) using the full conditional distributions given values of other parameters at step \(s\). In our case since missing values are unknown (just like parameters), we can replace the parameters in the Gibbs sampler with our missing values.

Instead of deriving full conditional distributions, MICE takes a shortcut for missing data by directly specifying full conditional distributions. One hopes that there is a proper posterior distribution which would yield these full conditional distributions.

A common choice is to specify a Bayesian linear regression model, in our to exploit relationships between variables. Let \(Y_{i1}\) denote the skin variable and \(Y_{i2}\) the BMI:

\begin{gathered} Y_{i1} \mid Y_{i2}, \beta^{(2)}, {\sigma^2}^{(2)} \sim \mathcal{N}(\beta_1^{(2)} + \beta_2^{(2)} Y_{i2}, {\sigma^2}^{(2)}) \\ Y_{i2} \mid Y_{i1}, \beta^{(1)}, {\sigma^2}^{(1)} \sim \mathcal{N}(\beta_1^{(1)} + \beta_2^{(1)} Y_{i1}, {\sigma^2}^{(1)}) \end{gathered}

where the two mean terms are linear regression of BMI on skin and vice versa. Then we specify appropriate priors on \(\beta^{(1)}\), \(\beta^{(2)}\), \({\sigma^2}^{(1)}\), and \({\sigma^2}^{(2)}\). This gives us a Bayesian model to estimate the missing skin values \(Y_{i1}\) and missing BMI values \(Y_{i2}\).

What’s a bit different from before is that we are only generating a small number of samples, because we are filling in missing data instead of parameter values. This is usually only run for 5-10 iterations.

Implementation in R

This imputation procedure can be done automatically in R using the mice() function in the mice package. Below we run \(D=5\) iterations to obtain \(D\) separate, imputed data sets:

1
2
3
4
5
6
7
8
9

set.seed(42)
D <- 5
diab_mice <- mice::mice(
  diab,            # dataset with missing values
  m = D,           # number of imputed datasets
  method = "norm"  # use Bayesian linear regression discussed above
)

Y_imp <- map(seq(D), ~mice::complete(diab_mice, .x))

The Y_imp variable is a list of five data frames, each with missing values filled in by the corresponding values in diab_mice$imp.

Just like with other MCMC methods, we can assess convergence by looking at trace plots. With only five iterations, we want to see them randomly dispersed among each other, e.g. it would be undesirable if the red curve was always above the other curves.

Figure 3: Mean and standard deviation of imputed values across the five generated data sets.

Drawing posterior samples

We can then treat each of the five imputed data sets separately, and use JAGS⁴ to generate separate posterior samples for the parameters of interest. The table below has 95% credible intervals for the specified parameters:

These values are fairly similar, mainly because the sample size is large, and there is a relatively small proportion of missing values.

Combining results

Generally, we do not care about posterior inference of parameters for each imputed data set. Instead, we would rather pool our results together to account for imputation variability.

In frequentist inference, this is done using Rubin's rules; in Bayesian inference, we just aggregate all of our posterior samples together across the imputed data sets. The credible intervals after aggregating posterior samples are given below.

From the 95% CI of \(\rho\), we conclude there is fairly strong positive correlation between skin and BMI.

Concluding remarks

For many inferential tasks, simple imputation methods can perform well! Advanced methods are necessary if one is focused on quantifying uncertainty in estimates correctly, particularly for data sets with low sample sizes or a large amount of missingness.

Missing data is very application-specific, so you should really think about what the variables represent, i.e. use domain knowledge and explore the data before considering how to impute. Missing data often does not occur completely randomly, e.g.:

• Non-response to survey questions – it might be a sensitive question, or maybe the survey is too long.
• Patient dropout in a longitudinal study – maybe the patients died or moved to somewhere else.

Finally, never, ever replace all of the missing values with zeros (or other arbitrary values)!