Multiple Imputation for Missing Data

Multiple imputation (MI) creates several plausible completed datasets by filling in missing values, analyzes each separately, and combines results using Rubin's rules.

Why Multiple Imputation?

The MICE Algorithm

Multiple Imputation by Chained Equations (MICE) is the most popular MI method.

Steps

| Step | Action |

|------|--------|

| 1 | Initialize missing values with simple imputation (e.g., mean) |

| 2 | For each variable with missing values: fit regression on other variables |

| 3 | Draw imputed values from the predictive distribution |

| 4 | Repeat steps 2-3 for many cycles (typically 10-20) |

| 5 | After convergence, save the imputed dataset |

| 6 | Repeat steps 1-5 M times to create M datasets |

Predictive Mean Matching (PMM)

Steps for each missing value:

1. Fit a regression predicting from other variables using observed data

2. Predict for both observed and missing cases

3. For each missing case, find the observed case with the closest predicted value

4. Use the observed value as the imputation

Rubin's Rules

Combined Point Estimate

Combined Variance

Number of Imputations

| Missing Fraction | Recommended M |

|-----------------|---------------|

| < 10% | 5-10 |

| 10-30% | 20-40 |

| 30-50% | 40-100 |

| > 50% | > 100 |

Diagnostics

Convergence

Plot the mean and variance of imputed values across iterations. They should stabilize after 10-20 iterations.

Imputation Quality

| Diagnostic | What to Check |

|-----------|---------------|

| Density plots | Imputed and observed distributions should overlap |

| Scatter plots | Relationships between variables should be similar |

| Trace plots | MICE chains should mix and converge |

Assumptions

Python Implementation

import numpy as np

import pandas as pd

from sklearn.experimental import enable_iterative_imputer

from sklearn.impute import IterativeImputer

import matplotlib.pyplot as plt



np.random.seed(42)



# Simulate data

n = 500

X1 = np.random.randn(n)

X2 = 0.6 * X1 + np.random.randn(n) * 0.5

X3 = 0.3 * X1 + 0.5 * X2 + np.random.randn(n) * 0.7



# Create MAR missingness in X1

missing_prob = 1 / (1 + np.exp(-(-1 + 0.8*X2)))

R = np.random.binomial(1, missing_prob)

X1_obs = X1.copy()

X1_obs[R == 1] = np.nan



df = pd.DataFrame({'X1': X1_obs, 'X2': X2, 'X3': X3})

true_mean = X1.mean()

obs_mean = np.nanmean(X1_obs)

missing_pct = df['X1'].isna().mean()



print(f"True X1 mean: {true_mean:.3f}")

print(f"Observed X1 mean: {obs_mean:.3f}")

print(f"Missing: {missing_pct:.1%}")



# Multiple Imputation (M = 30)

M = 30

estimates = []

for m in range(M):

    imputer = IterativeImputer(random_state=m, max_iter=20)

    imputed = imputer.fit_transform(df)

    estimates.append(imputed[:, 0].mean())



# Rubin's rules

Q_bar = np.mean(estimates)

B = np.var(estimates, ddof=1)

# Approximate within-imputation variance

U_bar = np.var(X1) / (n * (1 - missing_pct))



T = U_bar + (1 + 1/M) * B

SE = np.sqrt(T)



print(f"\nMI estimate: {Q_bar:.3f} (SE: {SE:.3f})")

print(f"95% CI: [{Q_bar - 1.96*SE:.3f}, {Q_bar + 1.96*SE:.3f}]")



# Plot convergence

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

axes[0].plot(estimates, 'o-')

axes[0].axhline(y=true_mean, color='red', linestyle='--', label='True')

axes[0].set_xlabel('Imputation')

axes[0].set_ylabel('X1 Mean')

axes[0].set_title('Imputation Convergence')

axes[0].legend()



# Density comparison

axes[1].hist(X1, bins=30, alpha=0.5, density=True, label='Observed')

axes[1].hist([e for e in estimates], bins=30, alpha=0.3, density=True, label='Imputed means')

axes[1].legend()

axes[1].set_title('Density Comparison')

plt.tight_layout()

plt.show()

Worked Example

Key Takeaways

Related Topics

• See Missing Data for mechanisms (MCAR, MAR, MNAR)

• See Bootstrap Methods for resampling methods

• See Regression Diagnostics for model checking