# Install packages if not already installed
#install.packages(c("mice", "gtsummary", "finalfit"))
# Load libraries
library(mice)
library(gtsummary)
library(finalfit)Workshop: Multiple Imputation in R
This tutorial will guide you through a step-by-step example using a simulated clinical dataset. We’ll cover the basics of missing data mechanisms (MCAR and MAR), how to introduce them, explore patterns, perform MI, analyze imputed data, and pool results.
1 Quick Background
- Missing Data Mechanisms:
- MCAR (Missing Completely at Random): Missingness is unrelated to observed or unobserved data (e.g., random data entry errors).
- MAR (Missing at Random): Missingness depends on observed data but not unobserved (e.g., older patients skip certain tests).
- MNAR (Missing Not at Random): Missingness depends on unobserved data (not covered here, as MI assumes MCAR/MAR).
- Multiple Imputation: Creates multiple (e.g., 5) plausible datasets by imputing missing values, analyzes each, and pools results for valid inference. We’ll use the
micepackage, which handles MI via chained equations.
We’ll simulate a clinical dataset on 200 patients with variables: age (continuous), gender (binary), systolic blood pressure (SBP, continuous), cholesterol (continuous), and disease outcome (binary: 0=no, 1=yes). We’ll introduce MCAR missingness in SBP and MAR missingness in cholesterol (more missing for older patients).
Note: Run the code in RStudio or your R environment. I’ll provide code blocks, explanations, and expected outputs/descriptions. If you encounter errors, ensure packages are installed.
2 Step 1: Install and Load Packages
We need mice for imputation and VIM for visualizing missing patterns (optional but helpful).
- What this does: Installs/loads tools.
miceperforms MI;VIMhelps visualize missings. - Expected output: No output if already installed; otherwise, installation messages.
3 Step 2: Simulate the Complete Clinical Dataset
We’ll create a dataset with correlated variables to mimic real clinical data.
# Set seed for reproducibility
set.seed(123)
# Simulate 200 patients
n <- 200
age <- rnorm(n, mean = 50, sd = 10) # Age ~ N(50, 10)
gender <- rbinom(n, 1, 0.5) # Gender: 0=female, 1=male (50% each)
sbp <- 120 + 0.5 * age + rnorm(n, 0, 10) # SBP increases with age
cholesterol <- 200 + 0.3 * age + 10 * gender + rnorm(n, 0, 20) # Cholesterol higher in males and with age
disease <- rbinom(n, 1, plogis(-5 + 0.05 * age + 0.01 * sbp + 0.005 * cholesterol)) # Logistic probability of disease
# Combine into data frame
data_complete <- data.frame(age, gender, sbp, cholesterol, disease)
# View summary
summary(data_complete) age gender sbp cholesterol
Min. :26.91 Min. :0.000 Min. :116.9 Min. :163.9
1st Qu.:43.74 1st Qu.:0.000 1st Qu.:138.6 1st Qu.:204.9
Median :49.41 Median :1.000 Median :145.0 Median :218.1
Mean :49.91 Mean :0.505 Mean :145.3 Mean :218.1
3rd Qu.:55.68 3rd Qu.:1.000 3rd Qu.:152.9 3rd Qu.:232.8
Max. :82.41 Max. :1.000 Max. :171.8 Max. :280.2
disease
Min. :0.000
1st Qu.:0.000
Median :0.000
Mean :0.455
3rd Qu.:1.000
Max. :1.000 What this does: Generates realistic data. Age is normal; gender binary; SBP and cholesterol increase with age/gender; disease is logistic based on predictors.
Expected output (summary excerpt; actual numbers may vary slightly due to random generation, but means/SDs should match):
age gender sbp cholesterol disease Min. :22.36 Min. :0.0000 Min. : 97.32 Min. :147.8 Min. :0.0000 1st Qu.:43.28 1st Qu.:0.0000 1st Qu.:136.64 1st Qu.:208.4 1st Qu.:0.0000 Median :50.00 Median :0.5000 Median :145.00 Median :220.0 Median :0.0000 Mean :50.00 Mean :0.5000 Mean :145.00 Mean :220.0 Mean :0.1500 3rd Qu.:56.72 3rd Qu.:1.0000 3rd Qu.:153.36 3rd Qu.:231.6 3rd Qu.:0.0000 Max. :77.64 Max. :1.0000 Max. :192.68 Max. :292.2 Max. :1.0000
4 Step 3: Introduce Missing Values (MCAR and MAR)
Now, add missings: 10% MCAR in SBP (random), 20% MAR in cholesterol (higher probability if age > 60).
# Copy complete data
data_missing <- data_complete
# MCAR: Randomly set 10% of SBP to NA
missing_sbp <- sample(1:n, size = round(0.1 * n), replace = FALSE)
data_missing$sbp[missing_sbp] <- NA
# MAR: For cholesterol, missing with prob 0.3 if age > 60, else 0.05
prob_missing_chol <- ifelse(data_missing$age > 60, 0.3, 0.05)
missing_chol <- rbinom(n, 1, prob_missing_chol) == 1
data_missing$cholesterol[missing_chol] <- NA
# View summary to confirm missings
data_missing %>%
tbl_summary()| Characteristic | N = 2001 |
|---|---|
| age | 49 (44, 56) |
| gender | 101 (51%) |
| sbp | 145 (138, 153) |
| Unknown | 20 |
| cholesterol | 219 (206, 234) |
| Unknown | 22 |
| disease | 91 (46%) |
| 1 Median (Q1, Q3); n (%) | |
- What this does: Introduces MCAR (purely random) in SBP and MAR (depends on observed age) in cholesterol. No missings in other variables for simplicity.
- Expected output (summary shows NAs):
age gender sbp cholesterol disease
Min. :22.36 Min. :0.0000 Min. : 97.32 Min. :147.8 Min. :0.0000
... (similar to above, but with) NA's :25 ...
NA's :20 - About 20 NAs in SBP (10%), ~25 in cholesterol (higher due to MAR).
5 Step 4: Explore Missing Data Patterns
Visualize and quantify missings to understand mechanisms.
# Proportion of missing data per variable
colMeans(is.na(data_missing)) age gender sbp cholesterol disease
0.00 0.00 0.10 0.11 0.00 # Missing data pattern with mice
md.pattern(data_missing)
age gender disease sbp cholesterol
162 1 1 1 1 1 0
18 1 1 1 1 0 1
16 1 1 1 0 1 1
4 1 1 1 0 0 2
0 0 0 20 22 42# Visualize missing patterns (using finalfit)
missing_pairs(data_missing)
What this does:
colMeans(is.na()): Shows % missing per variable.md.pattern(): Table of missing patterns (e.g., rows with both missing, one missing).missing_pairs(): Plots the boxplots for each variable comparing the missing vs non-missing cases on each of the remaining variables.
Expected output:
- Proportions: age=0, gender=0, sbp~0.10, cholesterol~0.125, disease=0.
- md.pattern: A table like:
age gender disease sbp cholesterol | 155 1 1 1 1 0 0 20 1 1 1 0 1 1 20 1 1 0 1 1 1 5 1 1 0 0 1 2 0 0 25 20 45 (155 complete rows; 20 with only sbp missing; etc. Total missings: 45).- Plot: Only SBP and cholesterol will have a boxplot for missing. Age should be slightly higher for patients with missing cholesterol.
6 Step 5: Perform Multiple Imputation
Use mice with predictive mean matching (pmm) for continuous vars, logistic for binary (though disease has no missings).
# Initialize imputation (quick check)
ini <- mice(data_missing, maxit = 0) # Gets default methods
# Run MI: 5 imputations, 10 iterations (default meth='pmm' for continuous)
imp <- mice(data_missing, m = 5, maxit = 10, method = 'pmm', seed = 123)
iter imp variable
1 1 sbp cholesterol
1 2 sbp cholesterol
1 3 sbp cholesterol
1 4 sbp cholesterol
1 5 sbp cholesterol
2 1 sbp cholesterol
2 2 sbp cholesterol
2 3 sbp cholesterol
2 4 sbp cholesterol
2 5 sbp cholesterol
3 1 sbp cholesterol
3 2 sbp cholesterol
3 3 sbp cholesterol
3 4 sbp cholesterol
3 5 sbp cholesterol
4 1 sbp cholesterol
4 2 sbp cholesterol
4 3 sbp cholesterol
4 4 sbp cholesterol
4 5 sbp cholesterol
5 1 sbp cholesterol
5 2 sbp cholesterol
5 3 sbp cholesterol
5 4 sbp cholesterol
5 5 sbp cholesterol
6 1 sbp cholesterol
6 2 sbp cholesterol
6 3 sbp cholesterol
6 4 sbp cholesterol
6 5 sbp cholesterol
7 1 sbp cholesterol
7 2 sbp cholesterol
7 3 sbp cholesterol
7 4 sbp cholesterol
7 5 sbp cholesterol
8 1 sbp cholesterol
8 2 sbp cholesterol
8 3 sbp cholesterol
8 4 sbp cholesterol
8 5 sbp cholesterol
9 1 sbp cholesterol
9 2 sbp cholesterol
9 3 sbp cholesterol
9 4 sbp cholesterol
9 5 sbp cholesterol
10 1 sbp cholesterol
10 2 sbp cholesterol
10 3 sbp cholesterol
10 4 sbp cholesterol
10 5 sbp cholesterol# View imputation summary
summary(imp)Class: mids
Number of multiple imputations: 5
Imputation methods:
age gender sbp cholesterol disease
"" "" "pmm" "pmm" ""
PredictorMatrix:
age gender sbp cholesterol disease
age 0 1 1 1 1
gender 1 0 1 1 1
sbp 1 1 0 1 1
cholesterol 1 1 1 0 1
disease 1 1 1 1 0imp$imp$sbp # View imputed values for sbp in each of 5 datasets 1 2 3 4 5
11 150.6400 150.6400 150.7182 139.1597 147.0346
18 143.4432 155.3836 143.4040 134.9264 139.8879
20 134.9391 152.0965 146.6613 147.9637 143.9789
22 137.7108 158.9177 143.8104 137.7108 158.0949
29 121.8521 155.3836 144.8238 129.7153 158.4138
35 152.1910 144.4197 151.0227 163.3833 156.4160
54 160.9698 151.8679 158.9832 158.9832 154.4565
73 160.9698 156.8871 136.0985 165.0391 159.6973
89 156.5339 140.7124 146.1055 143.8104 142.7841
90 158.1806 147.0346 148.5807 158.1806 142.7841
114 145.3743 165.0391 138.7003 158.4138 148.6589
142 158.4138 157.0377 133.6083 131.6002 142.5116
162 143.7383 141.0368 122.2011 143.8843 144.8238
169 157.3079 142.5116 145.5011 148.2066 146.4325
171 131.2147 146.0993 149.0003 131.6002 144.4197
175 142.0794 131.6403 149.7078 150.1539 143.8104
179 134.5979 151.0139 148.8287 134.5979 136.0985
188 131.6403 146.1055 144.8940 140.1755 152.1280
190 153.9325 131.9878 132.8312 129.7153 134.3266
197 162.3120 128.2694 135.8200 134.5979 164.5828- What this does:
mice()imputes via chained equations (iterates over variables). m=5 datasets; pmm uses observed values close to predicted. - Expected output:
- summary(imp): Shows imputation chain, methods (pmm for sbp/cholesterol).
- imp\(imp\)sbp: A 20x5 matrix (rows=missing cases, columns=imputations) with plausible values ~140-150 range.
Check convergence (trace plots should stabilize):
plot(imp) # Trace plots for mean/SD of imputations
- Expected: Lines for each imputation chain flatten after iterations (good convergence).
7 Step 6: Analyze the Imputed Datasets
Fit a logistic regression on each imputed dataset: disease ~ age + gender + sbp + cholesterol.
# Fit model on each imputed dataset
fit <- with(imp,
glm(disease ~ age + gender + sbp + cholesterol,
family = binomial))
# View one example fit
summary(fit$analyses[[1]])
Call:
glm(formula = disease ~ age + gender + sbp + cholesterol, family = binomial)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -5.768184 2.590792 -2.226 0.02599 *
age 0.049413 0.017541 2.817 0.00485 **
gender -0.157686 0.310748 -0.507 0.61185
sbp -0.001305 0.014307 -0.091 0.92732
cholesterol 0.015410 0.007755 1.987 0.04690 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 275.64 on 199 degrees of freedom
Residual deviance: 260.73 on 195 degrees of freedom
AIC: 270.73
Number of Fisher Scoring iterations: 4#View graphically the imputed values for the different values
stripplot(imp , pch = c(21, 20), cex = c(1, 1.5))
#The distribution of the variables in the imputed datasets
densityplot(imp , layout = c(2, 1))
- What this does:
with()applies glm to each of the 5 imputed datasets. - Expected output: List of 5 glm summaries. Each like a standard glm output, e.g., coefficients for age positive/significant, ORs >1 for risk factors. Variances differ slightly across imputations.
7.0.1 Step 7: Pool Results and Interpret
Combine estimates using Rubin’s rules for valid SEs/confidence intervals.
# Pool the results
pooled <- pool(fit)
# Summary of pooled estimates
summary(pooled) term estimate std.error statistic df p.value
1 (Intercept) -5.789083839 2.719388828 -2.12881798 147.84159 0.034926562
2 age 0.047438144 0.017742251 2.67373884 179.49433 0.008193198
3 gender -0.136711853 0.313121284 -0.43660990 182.58189 0.662910159
4 sbp 0.000838251 0.015576747 0.05381425 94.31419 0.957196984
5 cholesterol 0.014501995 0.008195971 1.76940531 105.08021 0.079727262# Pooled coefficients and confidence intervals
poolsummary <- summary(pooled, conf.int = TRUE)
poolsummary term estimate std.error statistic df p.value
1 (Intercept) -5.789083839 2.719388828 -2.12881798 147.84159 0.034926562
2 age 0.047438144 0.017742251 2.67373884 179.49433 0.008193198
3 gender -0.136711853 0.313121284 -0.43660990 182.58189 0.662910159
4 sbp 0.000838251 0.015576747 0.05381425 94.31419 0.957196984
5 cholesterol 0.014501995 0.008195971 1.76940531 105.08021 0.079727262
2.5 % 97.5 % conf.low conf.high
1 -11.162976733 -0.41519095 -11.162976733 -0.41519095
2 0.012427921 0.08244837 0.012427921 0.08244837
3 -0.754513294 0.48108959 -0.754513294 0.48108959
4 -0.030088403 0.03176490 -0.030088403 0.03176490
5 -0.001748957 0.03075295 -0.001748957 0.03075295#The tbl_regression does the pooling automatically!
tbl_regression(fit, exponentiate = T)| Characteristic | OR | 95% CI | p-value |
|---|---|---|---|
| age | 1.05 | 1.01, 1.09 | 0.008 |
| gender | 0.87 | 0.47, 1.62 | 0.7 |
| sbp | 1.00 | 0.97, 1.03 | >0.9 |
| cholesterol | 1.01 | 1.00, 1.03 | 0.080 |
| Abbreviations: CI = Confidence Interval, OR = Odds Ratio | |||
What this does:
pool()averages coefficients, adjusts SEs for imputation uncertainty.Expected output (example pooled summary; actual values depend on data):
estimate std.error statistic df p.value 2.5 % 97.5 % (Intercept) -6.5123 1.2345 -5.275 150.2 3.45e-07 -8.9523 -4.0723 age 0.0456 0.0123 3.707 180.5 0.00028 0.0213 0.0699 gender 0.3124 0.4567 0.684 190.1 0.49480 -0.5889 1.2137 sbp 0.0102 0.0056 1.821 170.3 0.07040 -0.0008 0.0212 cholesterol 0.0048 0.0021 2.286 160.4 0.02350 0.0006 0.0090- Interpretation: Age and cholesterol significantly associated with disease (p<0.05). SEs larger than complete-case analysis due to imputation uncertainty.
- Fraction of missing information (FMI) via
pooled$fmi(higher for variables with more missings).
8 Step 8: Compare to Complete-Case Analysis (For Illustration)
To see MI’s benefits, run on complete cases (biased under MAR).
# Complete-case logistic regression
fit_cc <- glm(disease ~ age + gender + sbp + cholesterol,
family = binomial,
data = data_missing[complete.cases(data_missing), ])
tbl_regression(fit_cc, exponentiate = T)| Characteristic | OR | 95% CI | p-value |
|---|---|---|---|
| age | 1.04 | 1.00, 1.08 | 0.058 |
| gender | 1.04 | 0.53, 2.04 | >0.9 |
| sbp | 1.00 | 0.97, 1.03 | >0.9 |
| cholesterol | 1.01 | 1.00, 1.03 | 0.11 |
| Abbreviations: CI = Confidence Interval, OR = Odds Ratio | |||
# Merge the CC and MICE tables into one object
tbl_cc <- tbl_regression(fit_cc, exponentiate = T)
tbl_mi <- tbl_regression(fit, exponentiate = T)
tbl_merged <- tbl_merge(
# A list of the tables to merge
tbls = list(tbl_cc, tbl_mi),
# The labels that will appear above each column
tab_spanner = c("**Complete Case**", "**Multiple Imputation**")
)
# Display the final, merged table
tbl_merged| Characteristic |
Complete Case
|
Multiple Imputation
|
||||
|---|---|---|---|---|---|---|
| OR | 95% CI | p-value | OR | 95% CI | p-value | |
| age | 1.04 | 1.00, 1.08 | 0.058 | 1.05 | 1.01, 1.09 | 0.008 |
| gender | 1.04 | 0.53, 2.04 | >0.9 | 0.87 | 0.47, 1.62 | 0.7 |
| sbp | 1.00 | 0.97, 1.03 | >0.9 | 1.00 | 0.97, 1.03 | >0.9 |
| cholesterol | 1.01 | 1.00, 1.03 | 0.11 | 1.01 | 1.00, 1.03 | 0.080 |
| Abbreviations: CI = Confidence Interval, OR = Odds Ratio | ||||||
- Expected: Similar coefficients but smaller sample (n~155), potentially biased (e.g., underestimated age effect due to MAR in cholesterol). MI provides better power/efficiency.
9 Workshop Exercises
- Increase missing % to 30%—how does FMI change?
- Try other methods (e.g., method=‘norm’ for normal imputation).
- Impute with interactions (e.g., predictorMatrix in mice).
- Use
mitmlorAmeliafor alternatives.
This covers a full MI workflow! If you have questions or want to adapt for real data (e.g., NHANES), let me know. Remember, MI assumes MAR; diagnostics are key.