Single imputation creates one completed dataset, which allows downstream analysis but fails to reflect uncertainty from missing values. Multiple imputation addresses this by generating several completed datasets, each representing a plausible set of replacements for missing values.

Multiple imputation is closely related to data augmentation: iterative simulation of missing values and parameters naturally produces multiple versions of $y_{\text{mis}}$.

The typical setting involves a regression model $p(y \mid X, \theta)$ with missing entries in $X$. A full Bayesian approach would require a joint model for $X$ given parameters $\psi$, producing inference for $(X_{\text{mis}}, \psi, \theta \mid X_{\text{obs}}, y)$. This can be difficult because it requires specifying a potentially complex distribution for the covariates. Practical missing-data methods simplify this by separating the imputation and analysis tasks.

Bayesian computation for missing data produces simulations $(y^{(s)}_{\text{mis}}, \theta^{(s)})$ for $s=1,\dots,S$. After obtaining these draws, two options exist:

  • Use the simulations directly to infer parameters, missing values, or predictions.
  • Provide $y^{(s)}_{\text{mis}}$ as multiple imputations for users who will analyze each completed dataset with standard methods.

In practice, the workflow is often split: first perform multiple imputation to clean the dataset, then conduct the final inference using the completed datasets.

Basic Steps of Multiple Imputation

  1. Jointly model $(X, y)$ and generate draws of $X^{(s)}_{\text{mis}}$ and $\theta^{(s)}$.
    After simulation, discard $\theta^{(s)}$ and retain completed datasets
    $X^{(s)} = (X_{\text{obs}}, X^{(s)}_{\text{mis}})$ for a small number of draws $s$.
  2. For each completed dataset $X^{(s)}$, fit the desired model $p(y \mid X^{(s)}, \theta)$
    treating the imputed values as known during this step.
  3. Combine the inferences from the $K$ completed datasets.
    In Bayesian analyses, simply mix the posterior draws across imputations.
    Analytic combining rules are also available.

Computation Using EM and Data Augmentation

Imputation typically begins with crude initial values, often based on simple assumptions such as MCAR. These initial values serve as starting points for EM or Gibbs sampling.

EM and Gibbs sampling both follow a repeated structure:
impute missing values, update parameters, and iterate.

Let $(y_{\text{obs}}, y_{\text{mis}})$ denote the observed and missing components; any latent variables (e.g., mixture indicators) are included in $y_{\text{mis}}$.

Correspondence to earlier notation:

  • Data: $y = (y_{\text{obs}}, I)$
  • Parameters: posterior mode of $\theta$
  • Missing components: averaged over $y_{\text{mis}}$

The EM algorithm focuses on the posterior mode. In exponential-family models, the expected complete-data log posterior requires only the expected sufficient statistics. Gibbs sampling replaces conditional expectations with simulated draws and targets the full posterior instead of a mode.

Nonignorable models

If the missingness mechanism is nonignorable, EM proceeds by explicitly conditioning on the inclusion indicator $I$. The expected complete-data log posterior is computed as a function of $(\theta, \phi)$, conditional on $(y_{\text{obs}}, I)$, averaging over $y_{\text{mis}}$ at current parameter values.

Monotone Missing Data Patterns

A dataset exhibits a monotone missingness pattern if variables can be ordered so that each block is at least as observed as the next. Many practical datasets approximate this pattern.

For monotone patterns, posterior mode computation is efficient:

  • Only one regression is needed per missingness pattern.
  • Values consistent with the monotone pattern can be imputed in a single block.
  • For nearly monotone datasets, apply EM in two parts:
    (a) handle non-monotone entries with standard E-steps,
    (b) apply efficient monotone-pattern updates to the remainder.

Inference With Multiple Imputations

Suppose $K$ completed datasets are produced.
For each dataset $k$, obtain:

  • parameter estimate $\hat{\theta}_k$
  • complete-data variance estimate $W_k$

The combined estimate isθˉK=1K∑k=1Kθ^k.\bar{\theta}_K = \frac{1}{K} \sum_{k=1}^K \hat{\theta}_k.

The variances include:

  • Within-imputation variance

WK=1K∑k=1KWk,W_K = \frac{1}{K} \sum_{k=1}^K W_k,

  • Between-imputation variance

BK=1K−1∑k=1K(θ^k−θˉK)2.B_K = \frac{1}{K – 1} \sum_{k=1}^K (\hat{\theta}_k – \bar{\theta}_K)^2.

The total variance isTK=WK+K+1KBK.T_K = W_K + \frac{K+1}{K} B_K.

Approximate interval estimates use a $t$ distribution with degrees of freedom:d.f.=(K−1)(1+KK+1WKBK)2.\text{d.f.} = (K – 1) \left( 1 + \frac{K}{K+1} \frac{W_K}{B_K} \right)^{2}.

If the fraction of missing information is low, results tend to be insensitive to missing-data assumptions. A common approach is to begin with an ignorable model and then examine sensitivity to plausible nonignorable mechanisms.