Multiple Imputation for Longitudinal RWE

In plain language

Multiple imputation is a method for handling missing values in a dataset without simply throwing away the patients who have gaps. Instead of guessing one single fill-in value, the method creates several (commonly 5 to 40) complete versions of the data, each with slightly different plausible values for the missing cells, then analyzes every version with the same statistical model and combines the results. The combining step — called Rubin's rules — intentionally makes the final confidence interval wider than if the data had been complete, because the width honestly reflects how uncertain we are about what the missing values actually were. The key assumption is that the missingness depends only on information we did observe (called missing at random, or MAR), not on the unobserved value itself.

Multiple imputation (MI)

does not recover information that the data lack; it propagates the uncertainty about missing values under an explicit, untestable assumption — almost always missing at random (MAR) conditional on observed data. The recipe is fixed: (1) specify an imputation model and draw `m` completed datasets that reflect both the predicted value and the residual/posterior uncertainty; (2) fit the identical substantive analysis (Cox, pooled logistic, GEE, a propensity model) in each completed dataset; (3) combine the `m` point estimates and within/between-imputation variances with Rubin's rules so the final standard error includes the extra variance from not knowing the missing values. In RWE the hard part is never the arithmetic — it is making MAR defensible when EHR labs are visit-driven, claims have structurally absent clinical constructs rather than item-level gaps, and longitudinal measures are missing for reasons tied to the very trajectory you want to estimate.

Core conceptual distinction

— two ideas are separable and both matter. (1) MI vs single imputation / complete-case (CCA): single imputation (mean, regression, LOCF) fixes one value and lies about the standard error; CCA discards partial records and is unbiased only under MCAR or, for a regression coefficient, when missingness in covariates is independent of the outcome given the model. MI is preferred when missingness depends on observed data (MAR) and auxiliary variables predict either the missing value or its missingness. (2) Imputation model vs substantive (analysis) model — congeniality: the imputation model must be at least as rich as the analysis model. If the substantive model is a Cox survival model, the imputation model must include the event indicator and the Nelson–Aalen estimate of the cumulative baseline hazard (White & Royston), and any interaction or non-linearity in the analysis must also live in the imputation step. The cleanest way to guarantee this is substantive-model-compatible FCS (SMC-FCS, Bartlett 2014), which imputes from a model derived from the analysis model itself. The estimand is unchanged by MI: you are still estimating the same hazard ratio, risk difference, or marginal mean you would target with complete data — MI only restores valid inference for it.

Pros, cons, and trade-offs

- vs complete-case analysis (CCA): MI uses partial records and auxiliary predictors, so it is more efficient and is unbiased under MAR where CCA is biased (e.g., labs missing more often in monitored, sicker patients). Cost: MI relies on a correctly specified, congenial imputation model and an unverifiable MAR assumption; under MNAR both CCA and MI are biased, sometimes in different directions. Prefer MI when covariate missingness is appreciable and plausibly MAR; prefer CCA when missingness is trivial, or when it is in the outcome and unrelated to covariates (CCA can then be nearly as efficient and more transparent). - vs single imputation / LOCF for repeated measures: MI propagates uncertainty; single imputation understates variance and LOCF additionally assumes a value never changes after the last visit — usually false for labs and PROs and a classic way to bias longitudinal effects toward the null or away from it depending on the dropout pattern. Prefer MI (or a likelihood/mixed-model approach) over LOCF essentially always. - vs likelihood-based / mixed models (MMRM, joint models, full-Bayesian): A correctly specified mixed model is valid under MAR without imputation and avoids imputation-model misspecification; MI is more flexible when missingness spans many heterogeneous variables, when the analysis is not a single likelihood (e.g., a propensity score then an outcome model), or when auxiliary variables outside the analysis model carry information. Prefer MMRM/joint models for a single repeated-measures outcome; prefer MI for multivariable covariate missingness feeding a downstream causal pipeline. - vs inverse-probability weighting for missingness: IPW is robust to imputation-model misspecification but discards partial records and is inefficient when many variables are partly missing; MI is more efficient but needs a correct imputation model. Doubly robust combinations exist but are rarely operationalized in routine RWE.

When to use

— appreciable missingness (rough rule: more than a few percent on variables that matter) in baseline covariates, labs, severity scores, PROs, or SDoH fields that you have good reason to treat as MAR given observed data, and you can name auxiliary variables (prior measurements, utilization, site, calendar time, treatment, the outcome) that predict the missing value or its missingness. Use MI to feed propensity or outcome models, and for repeated measures use a longitudinal-aware imputation (multilevel `2l.pmm`/`2l.norm`, or a wide layout that conditions on prior-visit values) rather than treating visits as independent.

When NOT to use — and when it is actively misleading or dangerous

- The construct is structurally absent, not item-missing. Imputing "ECOG" or a lab in claims that never records it is fabricating a variable from its correlates; the imputations carry no real information and falsely shrink the standard error. Do not impute a field that the data source does not collect — restrict, link, or drop the variable. - Missingness is informative (MNAR) and you treat it as MAR. If a PRO is missing because the patient felt too unwell, MAR-based MI biases toward the healthier observed distribution. Here MI under MAR is worse than honest — it launders an untestable assumption into a confident answer. Use delta-adjustment / pattern-mixture sensitivity analyses and report how conclusions move. - Imputing exposure or outcome timing when missingness reflects ascertainment failure. Imputing an index date, an event date, or death when the gap is unobserved follow-up (disenrollment, out-of-network care) manufactures person-time and can create or destroy immortal time. Handle with censoring/linkage, not imputation. - Uncongenial imputation. Imputing covariates without the outcome (or, for survival, without the event indicator and Nelson–Aalen cumulative hazard) biases associations toward the null. Always include the outcome in the imputation model.

Data-source operational depth

- Claims (FFS vs MA): Item-level lab/clinical missingness is rare because claims do not record those constructs at all — so MI is appropriate only for linked clinical variables, not for absent constructs (smoking, BMI, OTC use). A subtler trap: apparent "missing" diagnoses or fills are often unobserved person-time — Medicare Advantage enrollees lack fee-for-service claims, so a missing comorbidity may be an MA gap, not a true negative. Restrict to enrollees with the relevant benefit (A/B/D or commercial medical+pharmacy), and never impute over MA-only spans as if they were observed. - EHR: Missingness is visit-driven and informative — sicker patients are tested more, stable patients vanish. This makes MAR fragile unless the imputation model is rich in the drivers of measurement: site, provider, encounter type, prior measurement frequency, utilization, and the most recent observed value. A test ordered-but-pending vs never-ordered are different missingness mechanisms; encode them. Loss to follow-up is potentially informative; consider inverse probability of observation weights for the visit process rather than naive MI of every gap. - Registry: Stage, biomarker, and PRO missingness commonly varies by hospital, registry version, and calendar period (assays and staging systems change). Include site, calendar year, and registry-version indicators in the imputation model, and run site/calendar sensitivity analyses; pooling sites with different missingness mechanisms without these terms violates MAR. - Linked claims–EHR–registry: The strongest substrate — EHR/registry supply the clinical values, claims supply completeness and continuous enrollment — but linkage is selective (only the linkable subset) and the unlinked records may differ systematically. Treat linkage status as an auxiliary variable and check that the imputation model is stable across linked/unlinked strata.

Worked claims example

Question: 1-year all-cause mortality after initiating a nephrotoxic oral oncolytic, adjusting for baseline eGFR, in a linked commercial+Medicare FFS claims–EHR cohort. Cohort logic is claims-style: `person_id` with ≥365 days of continuous medical+pharmacy enrollment before the first qualifying `fill_date` (index_date), no MA-only person-time, and the arm/exposure taken from the index NDC. Baseline eGFR (from linked EHR labs in the [index_date−365, index_date] window) is missing for 35% of patients — and missing less often among patients with prior CKD diagnoses and high inpatient utilization, because they are monitored more. Complete-case Cox over-represents these monitored, high-risk patients and biases the eGFR–mortality association. The MI fix: (1) build a Nelson–Aalen estimate of the cumulative hazard from the observed survival times and the death indicator; (2) impute eGFR with predictive mean matching using the death indicator, the Nelson–Aalen term, treatment arm, age, prior CKD/dialysis flags, prior creatinine values, inpatient days, and calendar year as predictors — this enforces congeniality with the downstream Cox model; (3) create `m = 40` completed datasets (rule of thumb: `m` at least the percent missing); (4) fit the same Cox model `Surv(time, death) ~ arm + age + egfr + cci + prior_ip` in each dataset; (5) pool with Rubin's rules and report the pooled HR, its CI, and the fraction of missing information (FMI). (6) Because monitoring-driven missingness could be MNAR, add a delta-adjustment sensitivity analysis that shifts imputed eGFR downward (worse renal function) by a clinically meaningful delta and confirm the conclusion is stable. Report a missingness table, the imputation predictors, `m`, convergence/trace diagnostics, and the sensitivity results — and never collapse the `m` datasets into one averaged dataset before modeling (that destroys the between-imputation variance).

Interpreting the output

Consider the CKD progression example: eGFR is missing in 5 of 10 patients, m = 5 imputed datasets are created, and a Cox model for mortality is fitted in each. The five log-HR estimates are approximately −0.42, −0.38, −0.45, −0.40, and −0.45. Rubin's rules yield a pooled log-HR ≈ −0.42, corresponding to a pooled HR ≈ 0.66 (≈ 34% lower hazard in the treatment arm). The pooled SE is wider than any single imputation SE because it includes a between-imputation variance component reflecting uncertainty about the missing eGFR values themselves.

(1) Formal statistical interpretation. The pooled HR of ≈ 0.66 and its CI summarize the evidence across all m completed datasets under Rubin's rules. The between-imputation component inflates the pooled SE beyond what any single complete-dataset analysis would report — this inflation is correct and necessary: single imputation implicitly treats imputed values as known, understating uncertainty. The fraction of missing information (FMI) quantifies how much of the total variance is attributable to missingness; a high FMI signals that the result is sensitive to the imputation model and mechanism assumption.

(2) Practical interpretation for a decision-maker. The MI estimate of HR ≈ 0.66 uses all available patient-level data — including the five patients with missing eGFR — rather than discarding them as complete-case analysis would. The wider CI compared to a hypothetical complete-data analysis honestly reflects the residual uncertainty from the missingness. If the conclusion is sensitive to a MNAR sensitivity analysis (where imputed eGFR values are shifted to reflect worse kidney function), that finding should be reported alongside the primary MI result.

Worked example

Scenario

A researcher is studying whether a new oral oncolytic reduces one-year mortality compared with standard chemotherapy in a claims-linked EHR cohort of 500 patients. Baseline kidney function, measured as eGFR (estimated glomerular filtration rate, a blood test), is missing for 5 of a small demonstration group of 10 patients. Rather than drop those 5 patients, the analyst uses multiple imputation: five completed datasets are created, each with a different plausible eGFR for the missing patients drawn from a model that uses age, prior kidney-disease diagnosis, and treatment arm. A Cox proportional-hazards regression is then fit in each completed dataset, producing five separate hazard ratio estimates for the treatment effect. Those five estimates are pooled with Rubin's rules to give a single final hazard ratio and a confidence interval that is honestly wider because of the eGFR uncertainty.

Dataset

Raw analytic dataset (10 patients). eGFR is missing for 5 patients. The outcome model will be: Surv(time_days, death) ~ arm + age + eGfr.

Steps

• Run the imputation model (using observed age, arm, death indicator, and observed eGFR values as predictors) five separate times, each time drawing a fresh set of plausible eGFR values for the 5 missing patients. This produces 5 completed datasets.

• Fit the identical Cox regression model — Surv(time_days, death) ~ arm + age + egfr — in each of the 5 completed datasets independently. Each fit produces its own log hazard ratio for the treatment arm and its own standard error.

• Apply Rubin's rules: the pooled point estimate is simply the arithmetic mean of the 5 log hazard ratios. For example, if the 5 log-HR estimates are -0.42, -0.38, -0.45, -0.40, and -0.45, the pooled log-HR = (-0.42 + -0.38 + -0.45 + -0.40 + -0.45) / 5 = -2.10 / 5 = -0.42.

• Compute the within-imputation variance (average of the 5 squared standard errors) and the between-imputation variance (variance of the 5 point estimates across datasets). The total variance for the pooled estimate adds these two components plus a small correction term for finite m, inflating the standard error relative to a complete-data analysis.

• Convert the pooled log-HR to a hazard ratio: exp(-0.42) = 0.66, meaning the treatment arm has a 34 percent lower hazard of death. Report the pooled HR, its 95 percent confidence interval derived from the total variance, and the fraction of missing information (FMI) — which here would be moderate, reflecting that eGFR was missing for 50 percent of this small demonstration group.

Result

Pooled log-HR = (-0.42 + -0.38 + -0.45 + -0.40 + -0.45) / 5 = -0.42. Pooled HR = exp(-0.42) = 0.66 (95% CI wider than complete-data estimate because the between-imputation variance inflates the standard error). The 5 completed datasets were analyzed separately and combined with Rubin's rules; they were never averaged or collapsed into a single dataset before modeling.

Imputed Datasets Table

Each column shows the imputed eGFR drawn for the 5 missing patients in that completed dataset. Values differ across datasets because each draw adds residual uncertainty from the imputation model. The 5 completed datasets are kept separate; they are never averaged before modeling.

Runnable example

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.imputation import mice
from lifelines import NelsonAalenFitter

# 1) Nelson-Aalen cumulative hazard at each subject's follow-up time -> congenial survival imputation predictor.
naf = NelsonAalenFitter().fit(df["time"], event_observed=df["death"])
df["na_cumhaz"] = naf.cumulative_hazard_at_times(df["time"]).values

# 2) Build the MI dataset. Only egfr is imputed; death, na_cumhaz, arm and all covariates are predictors.
imp_cols = ["egfr", "death", "na_cumhaz", "arm", "age", "cci", "prior_ip"]
mi_data = mice.MICEData(df[imp_cols])
# Predictive-mean-matching-style draw for the continuous lab from all other columns.
mi_data.set_imputer("egfr", "egfr ~ death + na_cumhaz + arm + age + cci + prior_ip")

# 3) Fit the SAME Cox model in each completed dataset and pool with Rubin's rules.
#    statsmodels MICE expects a model formula; PHReg uses (time, status) via the `status` kwarg.
def cox_model(formula, data):
    return sm.PHReg.from_formula(
        "time ~ arm + age + egfr + cci + prior_ip",
        status=data["death"], data=data,
    )

mi = mice.MICE("time ~ arm + age + egfr + cci + prior_ip", cox_model, mi_data,
               n_skip=3)
results = mi.fit(n_burnin=10, n_imputations=40)   # m = 40 ~ percent missing
print(results.summary())                          # pooled log-HR, SE, CI, FMI per coefficient

library(mice)
library(survival)

## (A) Baseline MI congenial with Cox: add the event indicator + Nelson-Aalen cumulative hazard (White & Royston).
df_wide$na_cumhaz <- nelsonaalen(df_wide, time, death)

pred <- make.predictorMatrix(df_wide)
pred["egfr", c("person_id")] <- 0            # never predict from the id
imp <- mice(df_wide, m = 40, method = "pmm", predictorMatrix = pred, seed = 2026)
fits <- with(imp, coxph(Surv(time, death) ~ arm + age + egfr + cci + prior_ip))
summary(pool(fits))                          # pooled HR, CI, FMI via Rubin's rules

## (B) Longitudinal MI of a repeated lab using a multilevel imputer (2l.pmm) with person_id as the cluster.
meth <- make.method(df_long)
meth["lab"] <- "2l.pmm"
pm <- make.predictorMatrix(df_long)
pm["lab", "person_id"] <- -2                 # -2 flags the class (cluster) variable for 2l methods
pm["lab", c("prior_lab", "arm", "age", "time_visit")] <- 1
imp_long <- mice(df_long, m = 40, method = meth, predictorMatrix = pm, seed = 2026)
fit_long <- with(imp_long,
                 lme4::lmer(lab ~ arm + time_visit + age + (1 | person_id)))
summary(pool(fit_long))