9  Missing Data and Multiple Imputation

Open almost any electronic health record and you will find holes. A patient never had their cholesterol checked; a smoking-status field was left blank; a follow-up visit was missed. Missing values are not an occasional nuisance in clinical data — they are the default state. The tempting fix is to throw away every incomplete record and analyse only the “complete” patients. This chapter is largely about why that instinct is dangerous: dropping incomplete rows wastes hard-won data and, worse, can bias your results in a direction you cannot see. We will build up the modern alternative — multiple imputation — gently, from what makes data go missing through to fitting and reporting a model on imputed data.

9.1 Introduction

Suppose you are studying whether a biomarker predicts mortality, and the biomarker was measured in only 70% of patients. The quickest analysis simply drops the 30% with a missing value and runs the model on whoever remains. This is complete-case analysis (also called listwise deletion), and it is still the most common approach in clinical papers.

It has two problems. First, it is wasteful: a patient missing one value out of twenty is discarded entirely, even though their other nineteen values carry information. With several variables each partly missing, complete-case analysis can quietly delete the majority of your sample. Second, and more seriously, it can be biased: if the patients with missing values differ systematically from those without (and they usually do), the complete cases are no longer a fair representation of your population. For a short, clinically grounded overview of these trade-offs aimed at journal readers, see the JAMA Guide to Statistics and Methods piece by Newgard and Lewis (2015).

ImportantMissingness is rarely random in clinical data

A test is ordered because the clinician was worried; a follow-up is missed because the patient was too unwell to attend, or because they recovered and stopped coming. The very reason a value is missing is often linked to the patient’s condition and outcome. That is exactly the situation in which complete-case analysis goes wrong — and exactly why we need the tools in this chapter.

9.2 Why are data missing? MCAR, MAR, MNAR

Before fixing missing data we must think about why it is missing. Statisticians describe three mechanisms, and the right method depends on which one you are facing. The names are unfortunate jargon, so we define each in plain terms with a clinical example.

MCAR — Missing Completely At Random. The missingness has nothing to do with anything, observed or unobserved — as if a coin decided which values to erase. Clinical example: a lab machine breaks for a day, so anyone sampled that day has a missing result regardless of who they are. Under MCAR, the complete cases are a genuinely random subset, so complete-case analysis is unbiased (just wasteful). MCAR is rare in practice.

MAR — Missing At Random. The probability that a value is missing depends only on other variables you have observed, not on the missing value itself once those are accounted for. Clinical example: HbA1c is more often missing in younger patients — but among patients of the same age, whether HbA1c is missing has nothing to do with the HbA1c value. Because age (which you observed) explains the missingness, you can recover the lost information by using age in the imputation. MAR is the assumption under which multiple imputation works, and it is the realistic target for most analyses.

MNAR — Missing Not At Random. The probability of missingness depends on the missing value itself, even after accounting for everything you observed. Clinical example: patients with the highest blood pressure feel unwell and skip the clinic, so the most extreme readings are precisely the ones absent — and no observed variable fully explains this. MNAR is the hardest case; standard imputation can still be biased, and you must reason about the mechanism directly (see the sensitivity-analysis note at the end).

For a health-researcher’s walkthrough of these three mechanisms with worked observational-study examples, Pedersen et al. (2017) in Clinical Epidemiology is a clear and concise companion to this section.

flowchart TD
    Q["Why is the value missing?<br/>What does the chance of<br/>missingness depend on?"]
    Q --> A["MCAR: on nothing<br/>(pure chance, e.g. machine<br/>broke for a day)"]
    Q --> B["MAR: only on OBSERVED<br/>variables (e.g. missing more<br/>in younger patients)"]
    Q --> C["MNAR: on the MISSING<br/>value itself (e.g. sickest<br/>patients skip the visit)"]

Figure 9.1: The three missingness mechanisms. The key question is what the chance of being missing depends on.

NoteYou cannot fully test which mechanism you have

Whether data are MAR or MNAR depends on the unobserved values, so no statistical test can prove it from the data alone. The choice rests on clinical judgement about why values went missing. The pragmatic strategy: make missingness as close to MAR as you can by including rich predictors in the imputation, then test robustness with an MNAR sensitivity analysis.

9.3 From single to multiple imputation

Imputation means filling in a missing value with a plausible estimate. The naive version is single imputation — replace each gap with one number, most commonly the column’s mean (mean imputation). This is tempting but quietly damaging: it pretends you know the filled-in value with certainty. It shrinks the natural variability of the data, distorts correlations, and — crucially — gives standard errors that are too small, because the analysis never registers that you guessed. You end up overconfident, the same trap as complete-case analysis but dressed differently.

Multiple imputation (MI), developed by Rubin, fixes this by being honest about uncertainty. Instead of one filled-in value, it creates several complete datasets, each with the gaps filled in slightly differently to reflect the range of plausible values. You then analyse every dataset and combine the answers. The spread between the datasets becomes part of your final uncertainty — so the standard errors correctly include the cost of having had to guess. The JAMA Guide to Statistics and Methods article by Li et al. (2015) is a particularly accessible one-page introduction to this idea for clinical readers.

flowchart LR
    D["Incomplete<br/>dataset"] --> I["IMPUTE:<br/>create m complete<br/>datasets"]
    I --> A["ANALYSE:<br/>fit your model on<br/>each one separately"]
    A --> P["POOL:<br/>combine estimates with<br/>Rubin's rules"]
    P --> R["One result with<br/>honest standard errors"]

Figure 9.2: Multiple imputation in three steps: impute several complete datasets, analyse each one separately, then pool the results with Rubin’s rules.

9.3.1 MICE: imputation by chained equations

The most popular engine for multiple imputation is MICE — Multiple Imputation by Chained Equations (also called fully conditional specification). The intuition is a round-robin of small regressions: to fill in a variable’s gaps, MICE predicts it from all the other variables; it then moves to the next incomplete variable and does the same, looping through all variables repeatedly until the filled-in values settle. Because each variable is imputed from the others, MICE naturally preserves the relationships between them — a patient imputed to have high blood glucose will also tend to be imputed with the higher weight and HbA1c that usually accompany it. Repeating the whole loop with different random draws is what produces the several distinct completed datasets.

9.3.2 Rubin’s rules: pooling the results

Once you have fitted your model (say a logistic regression) on each of the \(m\) imputed datasets, you have \(m\) slightly different sets of estimates. Rubin’s rules are the recipe for combining them into one:

• The pooled estimate is simply the average of the estimates across the datasets.
• The pooled variance has two parts added together: the average within-dataset variance (your ordinary uncertainty) plus the variance between datasets (the extra uncertainty from not knowing the missing values). This second part is exactly what single imputation forgets.

The result is a single coefficient, standard error, confidence interval and p-value that properly reflect the missingness. Modern software does all of this for you.

9.3.3 How many imputations?

The classic advice was that 5 imputations suffice, and for modest missingness that is often fine. Current practice is more generous: a useful rule of thumb is to set the number of imputations \(m\) at least equal to the percentage of incomplete cases (so 30% missing → roughly \(m = 30\)). Imputations are cheap, so when in doubt use more (20–50 is common); more imputations make your results more stable and reproducible.

TipAlways include the outcome and helpful auxiliary variables

Two rules prevent the commonest MI mistakes. First, include the outcome variable in the imputation model. It feels wrong — “using the outcome to fill in predictors?” — but omitting it biases the very associations you are trying to estimate towards zero. Second, include auxiliary variables: extra variables not in your final analysis model but that help predict the missing values (e.g. a related lab test, an earlier measurement). They make the MAR assumption more plausible and sharpen the imputations. The imputation model should be at least as rich as your analysis model.

9.4 Worked example: complete-case vs multiple imputation

We simulate a clinical cohort, deliberately knock holes in it under a MAR mechanism, then compare three analyses: the full data (the truth we are trying to recover), complete-case analysis, and multiple imputation pooled with Rubin’s rules.

R

library(tidyverse) # data simulation and wrangling

set.seed(42)
n <- 800

# --- Simulate a complete clinical cohort ---
full <- tibble(
  age = rnorm(n, 60, 12),
  bmi = rnorm(n, 28, 5),
  sbp = 100 + 0.4 * age + 0.6 * bmi + rnorm(n, 0, 10), # systolic BP
  # Binary outcome: cardiovascular event, depends on age, bmi, sbp
  event = rbinom(n, 1, plogis(-6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp))
)

# --- Induce MAR missingness in BMI ---
# Older patients are MORE likely to have missing BMI (depends on age,
# an OBSERVED variable -> Missing At Random)
p_missing <- plogis(-2 + 0.05 * (full$age - 60))
miss <- rbinom(n, 1, p_missing) == 1

missing_data <- full
missing_data$bmi[miss] <- NA

cat("Proportion of BMI missing:", round(mean(miss), 3), "\n")
cat("Complete cases retained:", sum(!miss), "of", n, "\n")

What the code shows. We create 800 patients with age, BMI, systolic blood pressure and a binary cardiovascular event whose risk genuinely depends on all three. We then delete BMI for some patients, making deletion more likely in older patients — so missingness depends on age, an observed variable. That is textbook MAR. The printout reports what fraction of BMI is now missing and how many complete cases survive; notice that even moderate per-variable missingness can remove a sizeable chunk of the sample once you require complete records.

library(tidyverse)

set.seed(42)
n <- 800
full <- tibble(
  age = rnorm(n, 60, 12),
  bmi = rnorm(n, 28, 5),
  sbp = 100 + 0.4 * age + 0.6 * bmi + rnorm(n, 0, 10),
  event = rbinom(n, 1, plogis(-6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp))
)
p_missing <- plogis(-2 + 0.05 * (full$age - 60))
missing_data <- full
missing_data$bmi[rbinom(n, 1, p_missing) == 1] <- NA

# --- Reference: model on the FULL (pre-deletion) data ---
fit_full <- glm(event ~ age + bmi + sbp, data = full, family = binomial)

# --- Complete-case analysis (glm drops rows with any NA automatically) ---
fit_cc <- glm(event ~ age + bmi + sbp, data = missing_data, family = binomial)

cat("Full-data BMI coefficient:       ", round(coef(fit_full)["bmi"], 4), "\n")
cat("Complete-case BMI coefficient:   ", round(coef(fit_cc)["bmi"], 4), "\n")
cat("Complete-case n:", length(fit_cc$y), "of", n, "\n")

What the code shows. We fit the same logistic regression twice: once on the full data (our gold-standard answer) and once on the data with holes, where glm silently drops every patient with a missing BMI. Compare the two bmi coefficients and the sample size. The complete-case model is fitted on a smaller, older-skewed subset (because older patients were preferentially deleted), so its estimate is noisier and can drift away from the full-data value. This is the waste-and-bias problem made concrete — and the motivation for imputing instead of deleting.

library(tidyverse)
library(mice) # multiple imputation by chained equations

set.seed(42)
n <- 800
full <- tibble(
  age = rnorm(n, 60, 12),
  bmi = rnorm(n, 28, 5),
  sbp = 100 + 0.4 * age + 0.6 * bmi + rnorm(n, 0, 10),
  event = rbinom(n, 1, plogis(-6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp))
)
p_missing <- plogis(-2 + 0.05 * (full$age - 60))
missing_data <- full
missing_data$bmi[rbinom(n, 1, p_missing) == 1] <- NA

# --- Multiple imputation with mice ---
# m = 20 imputed datasets; the imputation model uses ALL variables,
# INCLUDING the outcome 'event' (important!).
imp <- mice(missing_data, m = 20, method = "pmm", seed = 123, printFlag = FALSE)

# --- Analyse each imputed dataset, then pool with Rubin's rules ---
fits <- with(imp, glm(event ~ age + bmi + sbp, family = binomial))
pooled <- pool(fits)
summary(pooled, conf.int = TRUE)

What the code shows. mice() builds 20 completed datasets using predictive mean matching (pmm, a robust default that fills gaps with observed values from similar patients, so imputed BMIs stay realistic). By default mice uses every column — including the outcome event — as a predictor, which is exactly what we want. with(imp, glm(...)) fits our logistic model separately on all 20 datasets, and pool() applies Rubin’s rules to combine them. The summary() table looks like an ordinary regression output — estimate, standard error, confidence interval, p-value — but every number now reflects both the within-dataset uncertainty and the between-dataset spread. Compare the pooled bmi estimate with the full-data and complete-case values from the previous chunk: the MI estimate should track the full-data truth more closely than complete-case analysis did, while using all 800 patients.

WarningCheck your imputations before trusting them

After running mice, inspect the diagnostics. plot(imp) shows whether the chained-equation loops have settled (the lines for different imputations should mingle, not drift or separate). densityplot(imp) overlays the distribution of imputed values on the observed ones; imputed values should look plausible and broadly similar in shape — wildly different distributions signal a problem with the imputation model.

Python

import numpy as np
import pandas as pd
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.linear_model import LogisticRegression
import statsmodels.api as sm

rng = np.random.default_rng(42)
n = 800

# --- Simulate a complete clinical cohort ---
age = rng.normal(60, 12, n)
bmi = rng.normal(28, 5, n)
sbp = 100 + 0.4 * age + 0.6 * bmi + rng.normal(0, 10, n)
lin = -6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp
event = rng.binomial(1, 1 / (1 + np.exp(-lin)))

full = pd.DataFrame({"age": age, "bmi": bmi, "sbp": sbp, "event": event})

# --- Induce MAR missingness in BMI (depends on observed age) ---
p_missing = 1 / (1 + np.exp(-(-2 + 0.05 * (age - 60))))
miss = rng.binomial(1, p_missing) == 1
data = full.copy()
data.loc[miss, "bmi"] = np.nan
print(f"Proportion BMI missing: {miss.mean():.3f}")

# --- Complete-case logistic regression ---
cc = data.dropna()
Xcc = sm.add_constant(cc[["age", "bmi", "sbp"]])
fit_cc = sm.Logit(cc["event"], Xcc).fit(disp=0)
print("\nComplete-case BMI coefficient:", round(fit_cc.params["bmi"], 4),
      "  (n =", len(cc), ")")

What the code shows. This mirrors the R setup: a complete cohort, then MAR missingness induced in BMI as a function of age, followed by a complete-case logistic regression via statsmodels. The printed bmi coefficient and reduced sample size show the same waste-and-bias issue as in R — the dropped patients are disproportionately older, so the surviving subset is not representative. Note we import enable_iterative_imputer to switch on scikit-learn’s experimental multiple-imputation machinery, used next.

import numpy as np
import pandas as pd
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.impute import IterativeImputer
import statsmodels.api as sm

rng = np.random.default_rng(42)
n = 800
age = rng.normal(60, 12, n); bmi = rng.normal(28, 5, n)
sbp = 100 + 0.4 * age + 0.6 * bmi + rng.normal(0, 10, n)
event = rng.binomial(1, 1 / (1 + np.exp(-(-6 + 0.04*age + 0.03*bmi + 0.02*sbp))))
data = pd.DataFrame({"age": age, "bmi": bmi, "sbp": sbp, "event": event})
p_missing = 1 / (1 + np.exp(-(-2 + 0.05 * (age - 60))))
data.loc[rng.binomial(1, p_missing) == 1, "bmi"] = np.nan

# --- Multiple imputation: run IterativeImputer m times with sampling on ---
# sample_posterior=True makes each run draw a DIFFERENT plausible fill-in,
# so we get genuinely multiple imputations (not one repeated answer).
# The outcome 'event' is INCLUDED among the imputation predictors.
m = 20
cols = ["age", "bmi", "sbp", "event"]
estimates, variances = [], []

for i in range(m):
    imputer = IterativeImputer(sample_posterior=True, random_state=i)
    completed = pd.DataFrame(imputer.fit_transform(data[cols]), columns=cols)
    X = sm.add_constant(completed[["age", "bmi", "sbp"]])
    fit = sm.Logit(completed["event"].round(), X).fit(disp=0)
    estimates.append(fit.params["bmi"])
    variances.append(fit.bse["bmi"] ** 2)

# --- Pool with Rubin's rules (by hand, to show the formula) ---
estimates = np.array(estimates); variances = np.array(variances)
q_bar = estimates.mean()                 # pooled estimate = average
u_bar = variances.mean()                 # within-imputation variance
b     = estimates.var(ddof=1)            # between-imputation variance
total_var = u_bar + (1 + 1 / m) * b      # Rubin's total variance
se = np.sqrt(total_var)

print(f"Pooled BMI coefficient: {q_bar:.4f}")
print(f"Pooled standard error:  {se:.4f}")
print(f"95% CI: ({q_bar - 1.96*se:.4f}, {q_bar + 1.96*se:.4f})")

What the code shows. scikit-learn’s IterativeImputer is the Python counterpart of MICE: it imputes each variable from the others in a loop. Setting sample_posterior=True and a different random_state each pass makes the 20 runs produce genuinely different completed datasets — the essence of multiple, not single, imputation. We include event among the imputation columns (the same “include the outcome” rule as in R), fit the logistic model on each completed dataset, and then apply Rubin’s rules by hand so the formula is visible: the pooled estimate is the average of the coefficients, and the total variance is the within-imputation variance plus an inflation for the between-imputation spread. The reported standard error is wider than any single imputation’s would be — correctly reflecting the uncertainty from the missing data. The pooled coefficient should land closer to the true full-data value than the complete-case estimate. (For a more turnkey, MICE-faithful workflow in Python, the miceforest and statsmodels MICE packages are popular alternatives.)

9.5 A note on MNAR and sensitivity analysis

Everything above assumes the data are MAR — that observed variables explain the missingness. If you suspect MNAR (the missing values are systematically different in a way no observed variable captures, e.g. the sickest patients dropping out), standard MI can still be biased, and you cannot detect this from the data alone. The responsible response is a sensitivity analysis: deliberately impute the missing values under pessimistic MNAR assumptions and see whether your conclusions hold.

The common technique is a delta adjustment (a pattern-mixture approach): after imputing, shift the imputed values by a fixed amount \(\delta\) to represent “the missing patients were systematically worse” — for example, assume dropouts had blood pressures 10 mmHg higher than the imputation predicted — and refit. Repeat across a plausible range of \(\delta\). If your main finding survives even under a strong, adverse shift, it is robust to MNAR; if it flips under a mild shift, you must be cautious. Reporting such a sensitivity analysis is increasingly expected in clinical trials and high-quality observational studies.

NoteAuxiliary variables that only predict missingness can backfire

A subtle 2024 finding: adding an auxiliary variable that predicts whether a value is missing but is unrelated to the value itself can, under MNAR, actually increase bias. The safe auxiliary variables are those genuinely correlated with the incomplete variable (e.g. an earlier measurement of the same quantity), not merely with the act of dropping out.

9.6 Exercises

TipExercise 1: Classify the mechanism (conceptual)

For each scenario, state whether the missingness is most plausibly MCAR, MAR, or MNAR, and justify in one sentence.

1. A weighing scale was out of service for two weeks; everyone seen in that window has missing weight.
2. Depression-score questionnaires are returned less often by male patients, but among men the chance of returning is unrelated to the score.
3. Patients with the most severe symptoms are too unwell to complete a quality-of-life survey, and severity is not otherwise recorded.

NoteSolution

R

# =============================================================================
# Chapter 6c, Exercise 1: Classify the mechanism (conceptual)
# Decide MCAR / MAR / MNAR for three clinical scenarios, with justification.
# =============================================================================

# This exercise is conceptual: the answers are reasoned below as structured
# comments. No data analysis is required, but we print the answers so the
# script produces visible output when run.

# -----------------------------------------------------------------------------
# (a) A weighing scale was out of service for two weeks; everyone seen in that
#     window has missing weight.
#
#     ANSWER: MCAR (Missing Completely At Random).
#     JUSTIFICATION: The equipment failure is external to the patients -- it has
#     nothing to do with their weight or any of their characteristics, so the
#     missing records are a genuinely random subset of the cohort.
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# (b) Depression-score questionnaires are returned less often by male patients,
#     but among men the chance of returning is unrelated to the score.
#
#     ANSWER: MAR (Missing At Random).
#     JUSTIFICATION: Missingness depends only on sex, an OBSERVED variable, and
#     is unrelated to the depression score itself once sex is accounted for, so
#     the mechanism is recoverable by conditioning on sex in the imputation.
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# (c) Patients with the most severe symptoms are too unwell to complete a
#     quality-of-life survey, and severity is not otherwise recorded.
#
#     ANSWER: MNAR (Missing Not At Random).
#     JUSTIFICATION: The probability of missingness depends on the unmeasured
#     severity, which is exactly what the missing quality-of-life value
#     reflects, and no observed variable captures it -- so the missingness
#     depends on the missing value itself.
# -----------------------------------------------------------------------------

answers <- data.frame(
  scenario = c("(a) scale out of service",
               "(b) fewer returns from men, unrelated to score",
               "(c) sickest skip QoL survey, severity unrecorded"),
  mechanism = c("MCAR", "MAR", "MNAR"),
  stringsAsFactors = FALSE
)

cat("=== Exercise 1: Missingness mechanism classification ===\n\n")
print(answers, row.names = FALSE)
cat("\nKey idea: the mechanism is defined by what the CHANCE of being missing\n")
cat("depends on -- nothing (MCAR), observed variables (MAR), or the missing\n")
cat("value itself (MNAR).\n")

Python

# =============================================================================
# Chapter 6c, Exercise 1: Classify the mechanism (conceptual)
# Decide MCAR / MAR / MNAR for three clinical scenarios, with justification.
# =============================================================================

import pandas as pd

# This exercise is conceptual: the answers are reasoned below as structured
# comments. No data analysis is required, but we print the answers so the
# script produces visible output when run.

# -----------------------------------------------------------------------------
# (a) A weighing scale was out of service for two weeks; everyone seen in that
#     window has missing weight.
#
#     ANSWER: MCAR (Missing Completely At Random).
#     JUSTIFICATION: The equipment failure is external to the patients -- it has
#     nothing to do with their weight or any of their characteristics, so the
#     missing records are a genuinely random subset of the cohort.
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# (b) Depression-score questionnaires are returned less often by male patients,
#     but among men the chance of returning is unrelated to the score.
#
#     ANSWER: MAR (Missing At Random).
#     JUSTIFICATION: Missingness depends only on sex, an OBSERVED variable, and
#     is unrelated to the depression score itself once sex is accounted for, so
#     the mechanism is recoverable by conditioning on sex in the imputation.
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# (c) Patients with the most severe symptoms are too unwell to complete a
#     quality-of-life survey, and severity is not otherwise recorded.
#
#     ANSWER: MNAR (Missing Not At Random).
#     JUSTIFICATION: The probability of missingness depends on the unmeasured
#     severity, which is exactly what the missing quality-of-life value
#     reflects, and no observed variable captures it -- so the missingness
#     depends on the missing value itself.
# -----------------------------------------------------------------------------

answers = pd.DataFrame({
    "scenario": ["(a) scale out of service",
                 "(b) fewer returns from men, unrelated to score",
                 "(c) sickest skip QoL survey, severity unrecorded"],
    "mechanism": ["MCAR", "MAR", "MNAR"],
})

print("=== Exercise 1: Missingness mechanism classification ===\n")
print(answers.to_string(index=False))
print("\nKey idea: the mechanism is defined by what the CHANCE of being missing")
print("depends on -- nothing (MCAR), observed variables (MAR), or the missing")
print("value itself (MNAR).")

TipExercise 2: Quantify the cost of complete-case analysis

Using the simulated cohort from this chapter, additionally induce 20% MAR missingness in sbp (in addition to the missing BMI).

1. How many complete cases remain once both variables have gaps?
2. Fit the complete-case logistic model and compare its age coefficient and standard error to the full-data model.
3. In one sentence, explain why dropping rows became much more costly once a second variable was incomplete.

NoteSolution

R

# =============================================================================
# Chapter 6c, Exercise 2: Quantify the cost of complete-case analysis
# Add 20% MAR missingness in sbp on top of missing BMI, then compare the
# complete-case age coefficient/SE to the full-data model.
# =============================================================================

library(tidyverse)   # data simulation and wrangling

set.seed(42)
n <- 800

# --- Simulate the chapter's complete clinical cohort ---
full <- tibble(
  age = rnorm(n, 60, 12),
  bmi = rnorm(n, 28, 5),
  sbp = 100 + 0.4 * age + 0.6 * bmi + rnorm(n, 0, 10),   # systolic BP
  event = rbinom(n, 1, plogis(-6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp))
)

# --- Induce MAR missingness in BMI (older patients more likely missing) ---
p_missing_bmi <- plogis(-2 + 0.05 * (full$age - 60))
miss_bmi      <- rbinom(n, 1, p_missing_bmi) == 1

# --- ADD ~20% MAR missingness in sbp, also depending on observed age ---
# Intercept chosen so the marginal missing fraction is about 20%.
p_missing_sbp <- plogis(-1.4 + 0.03 * (full$age - 60))
miss_sbp      <- rbinom(n, 1, p_missing_sbp) == 1

missing_data <- full
missing_data$bmi[miss_bmi] <- NA
missing_data$sbp[miss_sbp] <- NA

# --- (a) How many complete cases remain? ---
n_complete <- sum(complete.cases(missing_data))
cat("=== Exercise 2: Cost of complete-case analysis ===\n\n")
cat(sprintf("BMI missing:            %d (%.1f%%)\n",
            sum(miss_bmi), 100 * mean(miss_bmi)))
cat(sprintf("SBP missing:            %d (%.1f%%)\n",
            sum(miss_sbp), 100 * mean(miss_sbp)))
cat(sprintf("Complete cases (a):     %d of %d (%.1f%%)\n\n",
            n_complete, n, 100 * n_complete / n))

# --- (b) Full-data model vs complete-case model: age coefficient & SE ---
fit_full <- glm(event ~ age + bmi + sbp, data = full, family = binomial)
fit_cc   <- glm(event ~ age + bmi + sbp, data = missing_data, family = binomial)

sm_full <- summary(fit_full)$coefficients
sm_cc   <- summary(fit_cc)$coefficients

cat("--- (b) age coefficient (data-generating truth = 0.04) ---\n")
cat(sprintf("Full-data:      coef = %+.4f   SE = %.4f   (n = %d)\n",
            sm_full["age", "Estimate"], sm_full["age", "Std. Error"],
            length(fit_full$y)))
cat(sprintf("Complete-case:  coef = %+.4f   SE = %.4f   (n = %d)\n\n",
            sm_cc["age", "Estimate"], sm_cc["age", "Std. Error"],
            length(fit_cc$y)))
cat(sprintf("SE inflation (complete-case / full-data): %.2fx\n\n",
            sm_cc["age", "Std. Error"] / sm_full["age", "Std. Error"]))

# --- (c) Why dropping rows became much more costly ---
cat("--- (c) Comment ---\n")
cat("Requiring BOTH bmi and sbp to be present removes any patient missing\n")
cat("either one, so the two missing-data fractions compound -- the surviving\n")
cat("subset shrinks far more than either variable alone, and because both\n")
cat("gaps are age-driven the remainder is increasingly younger-skewed,\n")
cat("giving a smaller, less representative sample with larger standard errors.\n")

Python

# =============================================================================
# Chapter 6c, Exercise 2: Quantify the cost of complete-case analysis
# Add 20% MAR missingness in sbp on top of missing BMI, then compare the
# complete-case age coefficient/SE to the full-data model.
# =============================================================================

import numpy as np
import pandas as pd
import statsmodels.api as sm

rng = np.random.default_rng(42)
n = 800

# --- Simulate the chapter's complete clinical cohort ---
age = rng.normal(60, 12, n)
bmi = rng.normal(28, 5, n)
sbp = 100 + 0.4 * age + 0.6 * bmi + rng.normal(0, 10, n)      # systolic BP
lin = -6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp
event = rng.binomial(1, 1 / (1 + np.exp(-lin)))
full = pd.DataFrame({"age": age, "bmi": bmi, "sbp": sbp, "event": event})

# --- Induce MAR missingness in BMI (older patients more likely missing) ---
p_missing_bmi = 1 / (1 + np.exp(-(-2 + 0.05 * (age - 60))))
miss_bmi = rng.binomial(1, p_missing_bmi) == 1

# --- ADD ~20% MAR missingness in sbp, also depending on observed age ---
# Intercept chosen so the marginal missing fraction is about 20%.
p_missing_sbp = 1 / (1 + np.exp(-(-1.4 + 0.03 * (age - 60))))
miss_sbp = rng.binomial(1, p_missing_sbp) == 1

data = full.copy()
data.loc[miss_bmi, "bmi"] = np.nan
data.loc[miss_sbp, "sbp"] = np.nan

# --- (a) How many complete cases remain? ---
n_complete = int(data.dropna().shape[0])
print("=== Exercise 2: Cost of complete-case analysis ===\n")
print(f"BMI missing:            {miss_bmi.sum()} ({100 * miss_bmi.mean():.1f}%)")
print(f"SBP missing:            {miss_sbp.sum()} ({100 * miss_sbp.mean():.1f}%)")
print(f"Complete cases (a):     {n_complete} of {n} "
      f"({100 * n_complete / n:.1f}%)\n")


def fit_logit(df):
    X = sm.add_constant(df[["age", "bmi", "sbp"]])
    return sm.Logit(df["event"], X).fit(disp=0)


# --- (b) Full-data model vs complete-case model: age coefficient & SE ---
fit_full = fit_logit(full)
fit_cc = fit_logit(data.dropna())

print("--- (b) age coefficient (data-generating truth = 0.04) ---")
print(f"Full-data:      coef = {fit_full.params['age']:+.4f}   "
      f"SE = {fit_full.bse['age']:.4f}   (n = {int(fit_full.nobs)})")
print(f"Complete-case:  coef = {fit_cc.params['age']:+.4f}   "
      f"SE = {fit_cc.bse['age']:.4f}   (n = {int(fit_cc.nobs)})\n")
print(f"SE inflation (complete-case / full-data): "
      f"{fit_cc.bse['age'] / fit_full.bse['age']:.2f}x\n")

# --- (c) Why dropping rows became much more costly ---
print("--- (c) Comment ---")
print("Requiring BOTH bmi and sbp to be present removes any patient missing")
print("either one, so the two missing-data fractions compound -- the surviving")
print("subset shrinks far more than either variable alone, and because both")
print("gaps are age-driven the remainder is increasingly younger-skewed,")
print("giving a smaller, less representative sample with larger standard errors.")

TipExercise 3: Multiple imputation end to end

Using the same dataset with missing BMI and SBP:

1. Run multiple imputation with m = 30 (R mice, or Python IterativeImputer).
2. Fit the logistic model on each imputed dataset and pool with Rubin’s rules.
3. Compare the pooled estimates and standard errors to both the complete-case and full-data models. Which approach recovers the truth best?
4. Inspect the imputation diagnostics. Do the imputed values look plausible?

NoteSolution

R

# =============================================================================
# Chapter 6c, Exercise 3: Multiple imputation end to end
# m = 30 imputations with mice on the cohort with missing BMI and SBP; pool
# with Rubin's rules and compare pooled / complete-case / full-data estimates.
# =============================================================================

library(tidyverse)   # data simulation and wrangling
library(mice)        # multiple imputation by chained equations

set.seed(42)
n <- 800

# --- Simulate the chapter's complete clinical cohort ---
full <- tibble(
  age = rnorm(n, 60, 12),
  bmi = rnorm(n, 28, 5),
  sbp = 100 + 0.4 * age + 0.6 * bmi + rnorm(n, 0, 10),
  event = rbinom(n, 1, plogis(-6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp))
)

# --- Induce MAR missingness in BMI and SBP (both depend on observed age) ---
p_missing_bmi <- plogis(-2 + 0.05 * (full$age - 60))
p_missing_sbp <- plogis(-1.4 + 0.03 * (full$age - 60))
missing_data  <- full
missing_data$bmi[rbinom(n, 1, p_missing_bmi) == 1] <- NA
missing_data$sbp[rbinom(n, 1, p_missing_sbp) == 1] <- NA

cat("=== Exercise 3: Multiple imputation end to end ===\n\n")
cat(sprintf("Complete cases: %d of %d\n\n",
            sum(complete.cases(missing_data)), n))

# --- Reference models: full data and complete-case ---
fit_full <- glm(event ~ age + bmi + sbp, data = full,         family = binomial)
fit_cc   <- glm(event ~ age + bmi + sbp, data = missing_data, family = binomial)

# --- (a) Multiple imputation with m = 30 ---
# Percentage of incomplete cases is high, so we follow the "m >= % missing"
# rule of thumb. The imputation model uses ALL columns, INCLUDING the outcome.
imp <- mice(missing_data, m = 30, method = "pmm", seed = 123, printFlag = FALSE)

# --- (b) Fit the logistic model on each imputed set and pool (Rubin's rules) ---
fits   <- with(imp, glm(event ~ age + bmi + sbp, family = binomial))
pooled <- summary(pool(fits), conf.int = TRUE)
rownames(pooled) <- pooled$term

# --- (c) Compare pooled vs complete-case vs full-data (focus: age coef) ---
get_age <- function(fit) {
  s <- summary(fit)$coefficients
  c(coef = s["age", "Estimate"], se = s["age", "Std. Error"])
}
full_age <- get_age(fit_full)
cc_age   <- get_age(fit_cc)
pool_age <- c(coef = pooled["age", "estimate"], se = pooled["age", "std.error"])

cat("--- (c) age coefficient (data-generating truth = 0.04) ---\n")
cmp <- data.frame(
  method = c("Full data (truth)", "Complete-case", "MI pooled (m=30)"),
  coef   = c(full_age["coef"], cc_age["coef"], pool_age["coef"]),
  se     = c(full_age["se"],   cc_age["se"],   pool_age["se"])
)
print(format(cmp, digits = 4), row.names = FALSE)

cat("\nThe MI pooled estimate uses all 800 patients and should sit between the\n")
cat("complete-case value and the full-data truth, recovering the truth best\n")
cat("while giving honest standard errors (wider than a naive single imputation).\n\n")

# Also show the full pooled table (all coefficients).
cat("--- Full pooled summary (Rubin's rules) ---\n")
print(pooled[, c("term", "estimate", "std.error", "conf.low", "conf.high", "p.value")],
      row.names = FALSE, digits = 4)

# --- (d) Imputation diagnostics ---
cat("\n--- (d) Imputation diagnostics ---\n")
cat("Observed vs imputed summaries (mean [sd]) across the 30 imputations:\n")
diag_var <- function(var) {
  obs <- missing_data[[var]][!is.na(missing_data[[var]])]
  imp_vals <- unlist(lapply(seq_len(imp$m),
                            function(k) complete(imp, k)[[var]][is.na(missing_data[[var]])]))
  cat(sprintf("  %-4s observed: %6.2f [%.2f]   imputed: %6.2f [%.2f]\n",
              var, mean(obs), sd(obs), mean(imp_vals), sd(imp_vals)))
}
diag_var("bmi")
diag_var("sbp")
cat("Imputed means/spreads that track the observed ones indicate plausible\n")
cat("imputations. For interactive checks use plot(imp) (convergence) and\n")
cat("densityplot(imp) (imputed vs observed distributions).\n")

Python

# =============================================================================
# Chapter 6c, Exercise 3: Multiple imputation end to end
# m = 30 imputations with IterativeImputer on the cohort with missing BMI and
# SBP; pool with Rubin's rules and compare pooled / complete-case / full-data.
# =============================================================================

import numpy as np
import pandas as pd
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.impute import IterativeImputer
import statsmodels.api as sm

rng = np.random.default_rng(42)
n = 800

# --- Simulate the chapter's complete clinical cohort ---
age = rng.normal(60, 12, n)
bmi = rng.normal(28, 5, n)
sbp = 100 + 0.4 * age + 0.6 * bmi + rng.normal(0, 10, n)
event = rng.binomial(1, 1 / (1 + np.exp(-(-6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp))))
full = pd.DataFrame({"age": age, "bmi": bmi, "sbp": sbp, "event": event})

# --- Induce MAR missingness in BMI and SBP (both depend on observed age) ---
p_missing_bmi = 1 / (1 + np.exp(-(-2 + 0.05 * (age - 60))))
p_missing_sbp = 1 / (1 + np.exp(-(-1.4 + 0.03 * (age - 60))))
data = full.copy()
data.loc[rng.binomial(1, p_missing_bmi) == 1, "bmi"] = np.nan
data.loc[rng.binomial(1, p_missing_sbp) == 1, "sbp"] = np.nan

n_complete = int(data.dropna().shape[0])
print("=== Exercise 3: Multiple imputation end to end ===\n")
print(f"Complete cases: {n_complete} of {n}\n")


def fit_logit(df):
    X = sm.add_constant(df[["age", "bmi", "sbp"]])
    return sm.Logit(df["event"], X).fit(disp=0)


# --- Reference models: full data and complete-case ---
fit_full = fit_logit(full)
fit_cc = fit_logit(data.dropna())

# --- (a) Multiple imputation with m = 30 ---
# IterativeImputer is the sklearn counterpart of MICE. sample_posterior=True
# with a different random_state each pass yields genuinely different completed
# datasets (multiple, not single, imputation). The outcome 'event' is INCLUDED
# among the imputation predictors.
m = 30
cols = ["age", "bmi", "sbp", "event"]
coefs, ses = {"age": [], "bmi": [], "sbp": []}, {"age": [], "bmi": [], "sbp": []}
imputed_bmi, imputed_sbp = [], []
mask_bmi = data["bmi"].isna().to_numpy()
mask_sbp = data["sbp"].isna().to_numpy()

for i in range(m):
    imputer = IterativeImputer(sample_posterior=True, random_state=i)
    completed = pd.DataFrame(imputer.fit_transform(data[cols]), columns=cols)
    completed["event"] = data["event"].to_numpy()          # outcome is observed
    fit = fit_logit(completed)
    for v in ("age", "bmi", "sbp"):
        coefs[v].append(fit.params[v])
        ses[v].append(fit.bse[v])
    imputed_bmi.append(completed["bmi"].to_numpy()[mask_bmi])
    imputed_sbp.append(completed["sbp"].to_numpy()[mask_sbp])


# --- (b) Pool with Rubin's rules (implemented by hand) ---
def rubin(coef_list, se_list, m):
    q = np.array(coef_list)
    u = np.array(se_list) ** 2
    q_bar = q.mean()                       # pooled estimate = average of coefs
    u_bar = u.mean()                       # within-imputation variance
    b = q.var(ddof=1)                      # between-imputation variance
    total_var = u_bar + (1 + 1 / m) * b    # Rubin's total variance
    return q_bar, np.sqrt(total_var)


pooled = {v: rubin(coefs[v], ses[v], m) for v in ("age", "bmi", "sbp")}

# --- (c) Compare pooled vs complete-case vs full-data (focus: age coef) ---
print("--- (c) age coefficient (data-generating truth = 0.04) ---")
cmp = pd.DataFrame({
    "method": ["Full data (truth)", "Complete-case", "MI pooled (m=30)"],
    "coef": [fit_full.params["age"], fit_cc.params["age"], pooled["age"][0]],
    "se":   [fit_full.bse["age"],    fit_cc.bse["age"],    pooled["age"][1]],
})
print(cmp.to_string(index=False, float_format=lambda x: f"{x:.4f}"))
print("\nThe MI pooled estimate uses all 800 patients and should sit between the")
print("complete-case value and the full-data truth, recovering the truth best")
print("while giving honest standard errors (wider than a naive single imputation).\n")

print("--- Full pooled summary (Rubin's rules) ---")
full_tab = pd.DataFrame({
    "term": ["age", "bmi", "sbp"],
    "estimate": [pooled[v][0] for v in ("age", "bmi", "sbp")],
    "std.error": [pooled[v][1] for v in ("age", "bmi", "sbp")],
    "truth": [0.04, 0.03, 0.02],
})
print(full_tab.to_string(index=False, float_format=lambda x: f"{x:.4f}"))

# --- (d) Imputation diagnostics ---
print("\n--- (d) Imputation diagnostics ---")
print("Observed vs imputed summaries (mean [sd]) across the 30 imputations:")
for name, mask, imp_vals in (("bmi", mask_bmi, imputed_bmi),
                             ("sbp", mask_sbp, imputed_sbp)):
    obs = data[name].dropna().to_numpy()
    imp_all = np.concatenate(imp_vals)
    print(f"  {name:<4} observed: {obs.mean():6.2f} [{obs.std(ddof=1):.2f}]   "
          f"imputed: {imp_all.mean():6.2f} [{imp_all.std(ddof=1):.2f}]")
print("Imputed means/spreads that track the observed ones indicate plausible")
print("imputations (the sklearn analogue of mice's densityplot / convergence checks).")

TipExercise 4: MNAR sensitivity analysis (conceptual + code)

Suppose BMI is suspected MNAR: patients with very high BMI were less likely to be weighed.

1. Explain why standard MI (which assumes MAR) might underestimate the true association between BMI and the outcome here.
2. Describe a delta-adjustment sensitivity analysis: which way would you shift the imputed BMIs, and by how much?
3. Sketch (in words or code) how you would report the result across a range of delta values.

NoteSolution

R

# =============================================================================
# Chapter 6c, Exercise 4: MNAR sensitivity analysis (conceptual + code sketch)
# BMI suspected MNAR (high-BMI patients less likely to be weighed): reason
# about the bias and sketch a delta-adjustment sensitivity analysis.
# =============================================================================

library(tidyverse)
library(mice)

# -----------------------------------------------------------------------------
# (a) Why standard MI (which assumes MAR) may UNDERESTIMATE the association.
#
#     Standard MI fills the gaps using the observed data under a MAR model, so
#     imputed BMIs are drawn towards the observed (lower) range. If the truly
#     missing patients had systematically HIGHER BMI, the imputations are too
#     low and the upper tail of BMI -- the part most strongly linked to the
#     outcome -- is under-represented, so the fitted BMI-outcome association is
#     biased towards zero (attenuated).
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# (b) Delta-adjustment: direction and magnitude.
#
#     A delta (pattern-mixture) adjustment adds a fixed offset delta to the
#     imputed BMIs to represent "the unweighed patients were heavier than MAR
#     predicts". Here delta should be POSITIVE (shift imputed BMIs UP), because
#     the suspected mechanism removes high values. The magnitude spans a
#     clinically plausible range, e.g. 0 to +5 BMI units (0, +1, +2, +3, +5),
#     ideally anchored by external knowledge of how much heavier the missing
#     group is thought to be.
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# (c) Reporting across a range of delta values (illustrative code).
#     For each delta: impute under MAR, add delta to the imputed BMIs only,
#     refit the logistic model on each completed set, pool with Rubin's rules,
#     and tabulate the pooled BMI coefficient (+CI) as a function of delta.
#     A finding that stays clearly non-null across the range is robust to MNAR.
# -----------------------------------------------------------------------------

set.seed(42)
n <- 800
full <- tibble(
  age = rnorm(n, 60, 12),
  bmi = rnorm(n, 28, 5),
  sbp = 100 + 0.4 * age + 0.6 * bmi + rnorm(n, 0, 10),
  event = rbinom(n, 1, plogis(-6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp))
)
# MNAR-style deletion: higher BMI => more likely missing (for illustration).
p_missing_bmi <- plogis(-2 + 0.15 * (full$bmi - 28))
miss_bmi      <- rbinom(n, 1, p_missing_bmi) == 1
missing_data  <- full
missing_data$bmi[miss_bmi] <- NA

cat("=== Exercise 4: MNAR delta-adjustment sensitivity analysis ===\n\n")
cat(sprintf("BMI missing (MNAR mechanism): %d (%.1f%%)\n\n",
            sum(miss_bmi), 100 * mean(miss_bmi)))

# Impute ONCE (m = 20) under MAR, then re-use the imputations with each shift.
imp   <- mice(missing_data, m = 20, method = "pmm", seed = 123, printFlag = FALSE)
where_bmi <- is.na(missing_data$bmi)

deltas <- c(0, 1, 2, 3, 5)   # positive shifts on the BMI scale
results <- lapply(deltas, function(d) {
  shifted <- lapply(seq_len(imp$m), function(k) {
    dat <- complete(imp, k)
    dat$bmi[where_bmi] <- dat$bmi[where_bmi] + d   # delta adjustment
    dat
  })
  ests <- sapply(shifted, function(dat) coef(glm(event ~ age + bmi + sbp,
                                                 data = dat, family = binomial))["bmi"])
  vars <- sapply(shifted, function(dat) {
    s <- summary(glm(event ~ age + bmi + sbp, data = dat, family = binomial))
    s$coefficients["bmi", "Std. Error"]^2
  })
  qbar <- mean(ests)
  ubar <- mean(vars)
  b    <- var(ests)
  se   <- sqrt(ubar + (1 + 1 / imp$m) * b)         # Rubin's total variance
  c(delta = d, bmi_coef = qbar, se = se,
    lo = qbar - 1.96 * se, hi = qbar + 1.96 * se)
})
tab <- as.data.frame(do.call(rbind, results))

cat("Pooled BMI coefficient across delta (truth = 0.03):\n")
print(format(tab, digits = 4), row.names = FALSE)
cat("\nReading the table: as delta increases (assuming heavier unweighed\n")
cat("patients), the pooled BMI coefficient moves away from the attenuated\n")
cat("MAR value (delta = 0). If the coefficient and CI stay clearly positive\n")
cat("across the plausible delta range, the BMI-outcome association is robust\n")
cat("to this MNAR concern; if it collapses under a mild shift, interpret with\n")
cat("caution.\n")

Python

# =============================================================================
# Chapter 6c, Exercise 4: MNAR sensitivity analysis (conceptual + code sketch)
# BMI suspected MNAR (high-BMI patients less likely to be weighed): reason
# about the bias and sketch a delta-adjustment sensitivity analysis.
# =============================================================================

import numpy as np
import pandas as pd
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.impute import IterativeImputer
import statsmodels.api as sm

# -----------------------------------------------------------------------------
# (a) Why standard MI (which assumes MAR) may UNDERESTIMATE the association.
#
#     Standard MI fills the gaps using the observed data under a MAR model, so
#     imputed BMIs are drawn towards the observed (lower) range. If the truly
#     missing patients had systematically HIGHER BMI, the imputations are too
#     low and the upper tail of BMI -- the part most strongly linked to the
#     outcome -- is under-represented, so the fitted BMI-outcome association is
#     biased towards zero (attenuated).
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# (b) Delta-adjustment: direction and magnitude.
#
#     A delta (pattern-mixture) adjustment adds a fixed offset delta to the
#     imputed BMIs to represent "the unweighed patients were heavier than MAR
#     predicts". Here delta should be POSITIVE (shift imputed BMIs UP), because
#     the suspected mechanism removes high values. The magnitude spans a
#     clinically plausible range, e.g. 0 to +5 BMI units (0, +1, +2, +3, +5),
#     ideally anchored by external knowledge of how much heavier the missing
#     group is thought to be.
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# (c) Reporting across a range of delta values (illustrative code).
#     For each delta: impute under MAR, add delta to the imputed BMIs only,
#     refit the logistic model on each completed set, pool with Rubin's rules,
#     and tabulate the pooled BMI coefficient (+CI) as a function of delta.
#     A finding that stays clearly non-null across the range is robust to MNAR.
# -----------------------------------------------------------------------------

rng = np.random.default_rng(42)
n = 800
age = rng.normal(60, 12, n)
bmi = rng.normal(28, 5, n)
sbp = 100 + 0.4 * age + 0.6 * bmi + rng.normal(0, 10, n)
event = rng.binomial(1, 1 / (1 + np.exp(-(-6 + 0.04 * age + 0.03 * bmi + 0.02 * sbp))))
full = pd.DataFrame({"age": age, "bmi": bmi, "sbp": sbp, "event": event})

# MNAR-style deletion: higher BMI => more likely missing (for illustration).
p_missing_bmi = 1 / (1 + np.exp(-(-2 + 0.15 * (bmi - 28))))
mask_bmi = rng.binomial(1, p_missing_bmi) == 1
data = full.copy()
data.loc[mask_bmi, "bmi"] = np.nan

print("=== Exercise 4: MNAR delta-adjustment sensitivity analysis ===\n")
print(f"BMI missing (MNAR mechanism): {mask_bmi.sum()} ({100 * mask_bmi.mean():.1f}%)\n")


def fit_logit(df):
    X = sm.add_constant(df[["age", "bmi", "sbp"]])
    return sm.Logit(df["event"], X).fit(disp=0)


def rubin(coef_list, se_list, m):
    q = np.array(coef_list)
    u = np.array(se_list) ** 2
    q_bar = q.mean()
    total_var = u.mean() + (1 + 1 / m) * q.var(ddof=1)
    return q_bar, np.sqrt(total_var)


# Impute m sets ONCE under MAR, then re-use them with each delta shift.
m = 20
cols = ["age", "bmi", "sbp", "event"]
mask = data["bmi"].isna().to_numpy()
completed_sets = []
for i in range(m):
    imputer = IterativeImputer(sample_posterior=True, random_state=i)
    comp = pd.DataFrame(imputer.fit_transform(data[cols]), columns=cols)
    comp["event"] = data["event"].to_numpy()
    completed_sets.append(comp)

deltas = [0, 1, 2, 3, 5]      # positive shifts on the BMI scale
rows = []
for d in deltas:
    coefs, ses = [], []
    for comp in completed_sets:
        shifted = comp.copy()
        shifted.loc[mask, "bmi"] = shifted.loc[mask, "bmi"] + d   # delta adjustment
        fit = fit_logit(shifted)
        coefs.append(fit.params["bmi"])
        ses.append(fit.bse["bmi"])
    est, se = rubin(coefs, ses, m)
    rows.append({"delta": d, "bmi_coef": est, "se": se,
                 "lo": est - 1.96 * se, "hi": est + 1.96 * se})

tab = pd.DataFrame(rows)
print("Pooled BMI coefficient across delta (truth = 0.03):")
print(tab.to_string(index=False, float_format=lambda x: f"{x:.4f}"))
print("\nReading the table: as delta increases (assuming heavier unweighed")
print("patients), the pooled BMI coefficient moves away from the attenuated")
print("MAR value (delta = 0). If the coefficient and CI stay clearly positive")
print("across the plausible delta range, the BMI-outcome association is robust")
print("to this MNAR concern; if it collapses under a mild shift, interpret with")
print("caution.")

9.7 Summary

Missing data are ubiquitous in clinical research, and the instinctive fix — dropping incomplete records (complete-case analysis) — is both wasteful and potentially biased, because patients with missing values usually differ from those without. The three missingness mechanisms (MCAR, MAR, MNAR) describe why values are absent and determine what is safe to do; MAR, where missingness depends only on observed variables, is the realistic target and the assumption under which multiple imputation is valid. Single imputation (e.g. mean filling) is dangerous because it ignores the uncertainty of guessing. Multiple imputation, typically via MICE, instead creates several completed datasets, fits the analysis on each, and combines the results with Rubin’s rules, so the final standard errors honestly include the cost of the missing data. Always include the outcome and helpful auxiliary variables in the imputation model, use enough imputations (roughly the percentage missing), and, where MNAR is plausible, add an MNAR sensitivity analysis.

TipKey Takeaways

1. Complete-case analysis wastes data and can bias results — avoid it as a default.
2. MCAR (pure chance), MAR (depends on observed variables), MNAR (depends on the missing value itself); the mechanism dictates the method.
3. Single imputation (mean filling) is unsafe: it understates uncertainty and shrinks standard errors.
4. Multiple imputation = impute several times → analyse each → pool with Rubin’s rules.
5. Include the outcome and auxiliary variables in the imputation model; the imputation model should be at least as rich as the analysis model.
6. Use enough imputations (rule of thumb: at least the percentage of incomplete cases).
7. When MNAR is plausible, run a delta-adjustment sensitivity analysis — robust conclusions survive an adverse shift.

9.8 References and Further Reading

• For missing data theory and multiple imputation foundations, see Rubin (1987) and Little and Rubin (2019).
• For practical guidance on implementation, see Buuren (2018), Buuren and Groothuis-Oudshoorn (2011), Sterne et al. (2009), Newgard and Lewis (2015), Li et al. (2015), Pedersen et al. (2017), and Curnow et al. (2024).

Buuren, Stef van. 2018. Flexible Imputation of Missing Data. 2nd ed. Chapman; Hall/CRC.

Buuren, Stef van, and Karin Groothuis-Oudshoorn. 2011. “mice: Multivariate Imputation by Chained Equations in R.” Journal of Statistical Software 45 (3): 1–67. https://doi.org/10.18637/jss.v045.i03.

Curnow, Elinor, James R Carpenter, and Kate Tilling. 2024. “Multiple Imputation Using Auxiliary Imputation Variables That Only Predict Missingness Can Increase Bias Due to Data Missing Not at Random.” BMC Medical Research Methodology 24. https://doi.org/10.1186/s12874-024-02353-9.

Li, Peng, Elizabeth A Stuart, and David B Allison. 2015. “Multiple Imputation: A Flexible Tool for Handling Missing Data.” JAMA 314 (18): 1966–67. https://doi.org/10.1001/jama.2015.15281.

Little, Roderick J A, and Donald B Rubin. 2019. Statistical Analysis with Missing Data. 3rd ed. John Wiley & Sons.

Newgard, Craig D, and Roger J Lewis. 2015. “Missing Data: How to Best Account for What Is Not Known.” JAMA 314 (9): 940–41. https://doi.org/10.1001/jama.2015.10516.

Pedersen, Alma B, Ellen M Mikkelsen, Deirdre Cronin-Fenton, et al. 2017. “Missing Data and Multiple Imputation in Clinical Epidemiological Research.” Clinical Epidemiology 9: 157–66. https://doi.org/10.2147/CLEP.S129785.

Rubin, Donald B. 1987. Multiple Imputation for Nonresponse in Surveys. John Wiley & Sons.

Sterne, Jonathan A C, Ian R White, John B Carlin, et al. 2009. “Multiple Imputation for Missing Data in Epidemiological and Clinical Research: Potential and Pitfalls.” BMJ 338: b2393. https://doi.org/10.1136/bmj.b2393.