Compare Missing Data Imputation Methods

Simulate a clinical trial and quantify bias from three imputation strategies.

Introduction

30 Second Summary

Imagine a new pain medication is being tested, but some patients stop showing up to their appointments halfway through the trial. The conclusions your doctor relies on depend entirely on how researchers handle those gaps in the data.

In this project, you will simulate a two-arm clinical trial with 200 patients, inject realistic dropout patterns where sicker patients leave more often, then apply three different missing data strategies and measure how much each one distorts the true treatment effect.

What You'll Build

A Python simulation and analysis pipeline that demonstrates why imputation method choice is critical in clinical trials.

By the end of this project, you'll have:

  • A fully simulated longitudinal clinical trial dataset with realistic Missing At Random (MAR) dropout patterns.
  • Three imputed datasets (Complete Case, LOCF, Multiple Imputation) with fitted MMRM models producing treatment effect estimates.
  • A comparative bias analysis showing quantitatively how each imputation method distorts or preserves the true treatment effect.
  • Secret Mission: Add a Monte Carlo wrapper that runs 500 replications and plots the sampling distribution of bias under each method.

Are there any prerequisites?

You should be comfortable with Python basics (functions, loops, DataFrames) and have a general understanding of statistics (mean, variance, regression). Familiarity with clinical trial design concepts (treatment arms, placebo, time points) is helpful but not required.

Set Up Your Python Environment

This simulation study depends on several scientific Python packages for data generation, statistical modeling, and imputation. Before writing any code, you need the right version of Python and a clean environment to install those packages into.

In this step, you will verify that Python 3.11+ is installed, create an isolated project folder with a virtual environment, and install all seven required packages.

In this step, get ready to:
  • Verify Python 3.11+ is installed on your machine.
  • Create a project folder and virtual environment.
  • Install the required Python packages.
Verify Python 3.11+ is installed

This project requires Python 3.11 or higher because one of the key packages (scikit-learn 1.8.0) only supports Python 3.11+.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type Terminal (macOS) or Command Prompt (Windows) and press Enter to open it.
  • Check your Python version by running this command:

macOS

python3 --version

Windows

python --version

✔️ I see version 3.11 or higher

You should see output like Python 3.11.x or Python 3.12.x (or higher).

Your Python is ready. Continue to the next substep.

ⓧ I see an older version

If you see Python 3.9 or Python 3.10, your installed version is too old for this project.

On macOS, the system may ship with an older Python by default. Even after installing a newer version from python.org, the python3 command can still point to the older system version.

  • Download Python 3.11 or higher from python.org/downloads.
  • Run the installer and follow the prompts.
  • On macOS, after installing, try running python3.11 --version to confirm the new version is available.
  • If python3 still resolves to the old version on macOS, use python3.11 (the full versioned binary) for the remaining commands in this step.

ⓧ Command not found

Python is not installed on your machine yet.

  • Go to python.org/downloads and download the latest Python 3.11+ installer for your operating system.
  • Run the installer.
  • On Windows, check the box that says Add Python to PATH before clicking Install Now.
  • Close and reopen your terminal, then run the version check again to confirm it works.
Create a project folder and virtual environment

A virtual environment keeps this project's packages isolated from your system Python. This prevents version conflicts with other projects on your machine.

  • Move to your Desktop by running this command:

macOS

cd ~/Desktop

Windows

cd ~\Desktop
  • Create your project folder and move into it by running these commands:
mkdir missing-data-sim
cd missing-data-sim
  • Create and activate a virtual environment by running these commands:

macOS

python3 -m venv venv
source venv/bin/activate

Windows

python -m venv venv
venv\Scripts\activate

What does this do?

The python3 -m venv venv command creates a folder called venv inside your project. This folder contains a self-contained copy of Python and pip that is completely independent from your system installation.

Activating the environment tells your terminal to use this local copy of Python. Any packages you install while the environment is active go into venv/ and nowhere else.

You should see (venv) at the beginning of your terminal prompt. This confirms the virtual environment is active.

Don't see (venv) in your prompt?

  • Make sure you ran the activate command (not just the venv creation command).
  • On Windows, if you see an error about execution policies, run Set-ExecutionPolicy -Scope CurrentUser -ExecutionPolicy RemoteSigned first, then try activating again.
  • On macOS, confirm you used source venv/bin/activate (with source at the start).

Still stuck? Help me activate my virtual environment.

Install the required packages

This project uses seven packages for data manipulation, statistical modeling, machine learning imputation, and visualization.

  • Install all the required packages by running this command:
pip install numpy==2.4.6 pandas==3.0.3 statsmodels==0.14.6 scikit-learn==1.8.0 scipy==1.17.1 matplotlib==3.10.9 seaborn==0.13.2

What are these packages for?

  • numpy and scipy provide numerical computing and statistical functions.
  • pandas handles the tabular data structures you will use to represent patient records.
  • statsmodels fits the Mixed-Effects Model for Repeated Measures (MMRM).
  • scikit-learn provides the IterativeImputer class for Multiple Imputation.
  • matplotlib and seaborn produce the visualizations in the Secret Mission.

The install may take a minute or two. You should see pip downloading and installing each package with a final line reading Successfully installed followed by a list of all the packages.

Seeing an installation error?

  • Confirm your virtual environment is active (you should see (venv) in your prompt). If it is not active, pip may try to install globally and fail due to permissions.
  • If you see an error about no matching distribution for scikit-learn==1.8.0, your Python version is likely below 3.11. Go back to the first substep and upgrade Python.
  • If you see a network timeout, check your internet connection and try running the install command again.

Still stuck? Help me troubleshoot my pip install error.

Your environment is ready. Next up, you will write the simulation script that generates a realistic longitudinal clinical trial dataset with MAR dropout patterns.

Simulate Longitudinal Trial Data

Your environment is ready. Now you need the most critical ingredient for this study: a dataset where you know the exact truth.

You will simulate a two-arm pain trial with 200 patients measured at 4 time points. The true treatment effect is known (-1.5 units at the final visit), so later you can measure precisely how much each imputation method distorts that truth. You will also inject realistic Missing At Random (MAR) dropout, where patients experiencing more pain are more likely to drop out.

In this step, get ready to:
  • Build a data-generating function that produces longitudinal pain scores with a known treatment effect.
  • Implement a MAR dropout mechanism where higher pain increases the probability of missing subsequent visits.
  • Run the simulation and verify dropout rates by visit and treatment arm.
Create the data-generating function

Clinical trials measure patients repeatedly over time. Your simulation models 200 patients (100 per arm) across 4 visits: baseline, week 4, week 8, and week 12.

Each patient gets a random baseline pain level (random intercept), pain naturally decreases slightly over time, and the treatment arm receives an additional effect that grows linearly to -1.5 at the final visit.

  • Create a new file called simulate_trial.py in your missing-data-sim folder, then paste in the following code:
import numpy as np
import pandas as pd


def simulate_complete_data(n_patients=200, n_timepoints=4, true_treatment_effect=-1.5, seed=42):
    rng = np.random.default_rng(seed)

    patient_ids = np.arange(n_patients)
    treatments = np.array([0] * (n_patients // 2) + [1] * (n_patients // 2))

    # Random intercept per patient (baseline pain level varies)
    random_intercepts = rng.normal(loc=5.0, scale=1.5, size=n_patients)

    # Time effect: natural improvement of -0.3 per visit for placebo
    time_effect_per_visit = -0.3

    # Treatment effect grows linearly across visits (0 at baseline)
    treatment_effects = np.array([0, true_treatment_effect / 3,
                                   2 * true_treatment_effect / 3,
                                   true_treatment_effect])

    records = []
    for i in range(n_patients):
        for t in range(n_timepoints):
            pain = (
                random_intercepts[i]
                + time_effect_per_visit * t
                + treatments[i] * treatment_effects[t]
                + rng.normal(0, 0.8)  # within-subject residual
            )
            records.append({
                "patient_id": patient_ids[i],
                "treatment": treatments[i],
                "visit": t,
                "pain_score": pain
            })

    df = pd.DataFrame(records)
    return df

What does this code do?

  • Each patient starts with a random baseline pain level drawn from a normal distribution centered at 5.0. This simulates natural patient variability.
  • The treatment_effects array controls how the treatment advantage grows over time. At baseline it is 0, then scales linearly to -1.5 at visit 3 (week 12).
  • The residual noise (rng.normal(0, 0.8)) represents within-subject variability. This mimics real trial data where the same patient's pain fluctuates between visits.

Why simulate instead of using real data?

Real clinical trial data is confidential and the true treatment effect is unknown. By simulating, you set the ground truth yourself (-1.5 units). Later, when an imputation method estimates -0.8 instead of -1.5, you can measure exactly how much bias it introduced.

You won't see output yet. The next chunk adds the dropout mechanism, then you will run everything together.

Implement the MAR dropout mechanism

In real trials, patients do not drop out randomly. Patients experiencing worse symptoms are more likely to leave the study. This is the Missing At Random (MAR) assumption: the probability of missingness depends on observed data (current pain), not on future unobserved values.

  • In simulate_trial.py, add the following function below your simulate_complete_data function:
def inject_mar_missingness(df_complete, seed=42):
    rng = np.random.default_rng(seed + 1)
    df = df_complete.copy()

    patient_ids = df["patient_id"].unique()
    dropped_out = set()

    for visit in range(1, df["visit"].max() + 1):
        for pid in patient_ids:
            if pid in dropped_out:
                # Already dropped out, mark as missing
                df.loc[(df["patient_id"] == pid) & (df["visit"] == visit), "pain_score"] = np.nan
                continue

            # Get current pain score to determine dropout probability
            current_pain = df.loc[
                (df["patient_id"] == pid) & (df["visit"] == visit), "pain_score"
            ].values[0]

            # Logistic dropout probability: higher pain -> higher dropout
            # Calibrated to produce ~30% overall dropout by final visit
            logit_p = -3.0 + 0.4 * current_pain
            dropout_prob = 1 / (1 + np.exp(-logit_p))

            if rng.random() < dropout_prob:
                dropped_out.add(pid)
                df.loc[(df["patient_id"] == pid) & (df["visit"] == visit), "pain_score"] = np.nan
                # Set all future visits to NaN (monotone dropout)
                for future_visit in range(visit + 1, df["visit"].max() + 1):
                    df.loc[(df["patient_id"] == pid) & (df["visit"] == future_visit), "pain_score"] = np.nan

    return df

What does this dropout mechanism do?

  • For each post-baseline visit, the function computes a dropout probability using a logistic function: logit_p = -3.0 + 0.4 * current_pain. A patient with pain score 5 has about a 12% chance of dropping out at each visit. A patient with pain score 8 has about a 33% chance.
  • Once a patient drops out, all subsequent visits are set to NaN. This creates a monotone missing pattern, which is the most common dropout structure in longitudinal clinical trials.
  • Since the treatment arm lowers pain, treatment patients are less likely to drop out. This creates differential dropout between arms, which is exactly what makes imputation method choice so consequential.

One more chunk to go. The next section adds the execution block that ties both functions together and produces output.

Generate data and verify dropout patterns

Now you will add the main execution block that calls both functions, prints dropout statistics, and saves the datasets to CSV for the next step.

  • At the bottom of simulate_trial.py, add the following execution block:
if __name__ == "__main__":
    print("=" * 60)
    print("LONGITUDINAL TRIAL SIMULATION WITH MAR DROPOUT")
    print("=" * 60)

    # Generate complete data
    df_complete = simulate_complete_data()
    print(f"\nComplete data: {len(df_complete)} observations, {df_complete['patient_id'].nunique()} patients")

    # Inject missingness
    df_observed = inject_mar_missingness(df_complete)

    # Summary statistics
    print("\n--- Dropout Summary ---")
    for visit in range(4):
        visit_data = df_observed[df_observed["visit"] == visit]
        n_missing = visit_data["pain_score"].isna().sum()
        n_total = len(visit_data)
        print(f"Visit {visit}: {n_missing}/{n_total} missing ({100*n_missing/n_total:.1f}%)")

    print("\n--- Dropout by Treatment Arm ---")
    for arm in [0, 1]:
        arm_data = df_observed[(df_observed["treatment"] == arm) & (df_observed["visit"] == 3)]
        n_missing = arm_data["pain_score"].isna().sum()
        n_total = len(arm_data)
        label = "Placebo" if arm == 0 else "Treatment"
        print(f"{label}: {n_missing}/{n_total} missing at final visit ({100*n_missing/n_total:.1f}%)")

    # Save for use in analysis script
    df_complete.to_csv("trial_complete.csv", index=False)
    df_observed.to_csv("trial_observed.csv", index=False)
    print("\nData saved to trial_complete.csv and trial_observed.csv")

What does this block do?

  • Calls simulate_complete_data() to generate 800 observations (200 patients x 4 visits) with no missing values.
  • Passes the complete data into inject_mar_missingness() to create the realistic dropout pattern.
  • Prints dropout rates by visit (expect increasing missingness over time) and by treatment arm (expect more placebo dropout since those patients have higher pain).
  • Saves both datasets to trial_complete.csv and trial_observed.csv for use in the next step.
  • Save simulate_trial.py.
  • Run the script by entering this command in your terminal:
python simulate_trial.py
  • You should see output showing 800 total observations, dropout rates increasing from 0% at visit 0 to roughly 25-35% at visit 3, and higher dropout in the placebo arm than the treatment arm.

What should the dropout pattern look like?

Visit 0 (baseline) should have 0% missing because all patients are measured before randomization. Each subsequent visit should show increasing missingness as more patients leave the study.

The placebo arm should have more dropout than the treatment arm at the final visit. This is because placebo patients have higher pain on average (no treatment benefit), so the MAR mechanism removes them at a higher rate. This differential dropout is what makes naive imputation methods biased.

Seeing a ModuleNotFoundError?

  • Make sure your virtual environment is activated. You should see (venv) at the start of your terminal prompt. If not, activate it again from inside your missing-data-sim folder.
  • Confirm that numpy and pandas are installed by running pip list and checking they appear in the output.

Why is my simulation script throwing a ModuleNotFoundError?

  • Confirm that two new CSV files appeared by running:
ls
  • You should see trial_complete.csv and trial_observed.csv alongside your simulate_trial.py file.

✔️ Awesome, I've got everything!

Great. Double check you have saved simulate_trial.py and that both CSV files exist in your folder.

ⓧ I'd like to double check the full code

Here is the complete simulate_trial.py file. Compare it against yours:

import numpy as np
import pandas as pd


def simulate_complete_data(n_patients=200, n_timepoints=4, true_treatment_effect=-1.5, seed=42):
    rng = np.random.default_rng(seed)
    patient_ids = np.arange(n_patients)
    treatments = np.array([0] * (n_patients // 2) + [1] * (n_patients // 2))
    random_intercepts = rng.normal(loc=5.0, scale=1.5, size=n_patients)
    time_effect_per_visit = -0.3
    treatment_effects = np.array([0, true_treatment_effect / 3,
                                   2 * true_treatment_effect / 3,
                                   true_treatment_effect])
    records = []
    for i in range(n_patients):
        for t in range(n_timepoints):
            pain = (
                random_intercepts[i]
                + time_effect_per_visit * t
                + treatments[i] * treatment_effects[t]
                + rng.normal(0, 0.8)
            )
            records.append({
                "patient_id": patient_ids[i],
                "treatment": treatments[i],
                "visit": t,
                "pain_score": pain
            })
    df = pd.DataFrame(records)
    return df


def inject_mar_missingness(df_complete, seed=42):
    rng = np.random.default_rng(seed + 1)
    df = df_complete.copy()
    patient_ids = df["patient_id"].unique()
    dropped_out = set()
    for visit in range(1, df["visit"].max() + 1):
        for pid in patient_ids:
            if pid in dropped_out:
                df.loc[(df["patient_id"] == pid) & (df["visit"] == visit), "pain_score"] = np.nan
                continue
            current_pain = df.loc[
                (df["patient_id"] == pid) & (df["visit"] == visit), "pain_score"
            ].values[0]
            logit_p = -3.0 + 0.4 * current_pain
            dropout_prob = 1 / (1 + np.exp(-logit_p))
            if rng.random() < dropout_prob:
                dropped_out.add(pid)
                df.loc[(df["patient_id"] == pid) & (df["visit"] == visit), "pain_score"] = np.nan
                for future_visit in range(visit + 1, df["visit"].max() + 1):
                    df.loc[(df["patient_id"] == pid) & (df["visit"] == future_visit), "pain_score"] = np.nan
    return df


if __name__ == "__main__":
    print("=" * 60)
    print("LONGITUDINAL TRIAL SIMULATION WITH MAR DROPOUT")
    print("=" * 60)
    df_complete = simulate_complete_data()
    print(f"\nComplete data: {len(df_complete)} observations, {df_complete['patient_id'].nunique()} patients")
    df_observed = inject_mar_missingness(df_complete)
    print("\n--- Dropout Summary ---")
    for visit in range(4):
        visit_data = df_observed[df_observed["visit"] == visit]
        n_missing = visit_data["pain_score"].isna().sum()
        n_total = len(visit_data)
        print(f"Visit {visit}: {n_missing}/{n_total} missing ({100*n_missing/n_total:.1f}%)")
    print("\n--- Dropout by Treatment Arm ---")
    for arm in [0, 1]:
        arm_data = df_observed[(df_observed["treatment"] == arm) & (df_observed["visit"] == 3)]
        n_missing = arm_data["pain_score"].isna().sum()
        n_total = len(arm_data)
        label = "Placebo" if arm == 0 else "Treatment"
        print(f"{label}: {n_missing}/{n_total} missing at final visit ({100*n_missing/n_total:.1f}%)")
    df_complete.to_csv("trial_complete.csv", index=False)
    df_observed.to_csv("trial_observed.csv", index=False)
    print("\nData saved to trial_complete.csv and trial_observed.csv")

Your simulated trial data is ready. You have a complete dataset with the known truth, and an observed dataset with realistic MAR dropout. Next up, you will apply three different imputation strategies to handle those missing values and see which one recovers the true treatment effect.

Implement Three Imputation Strategies

You now have a simulated trial dataset with realistic MAR dropout. The next challenge is deciding what to do about those missing values.

Different imputation strategies make fundamentally different assumptions about why data is missing. In this step, you will implement three approaches in a new analysis file and see how each one reconstructs the gaps in your observed data.

In this step, get ready to:
  • Filter to complete cases only and observe the selection bias this introduces.
  • Carry the last observed value forward for each patient who drops out.
  • Run Multiple Imputation using chained equations to generate 20 plausible completed datasets.
Create the analysis file with Complete Case Analysis

The simplest way to handle missing data is to throw it away. Complete Case Analysis keeps only patients who have observations at every single time point. The problem is that dropout is not random. Patients with higher pain left the study, so the remaining "complete" patients are systematically healthier.

  • Create a new file called impute_and_analyze.py in your missing-data-sim folder.
  • Add the following imports at the top of the file:
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.impute import IterativeImputer
from sklearn.linear_model import BayesianRidge

What do these imports do?

  • numpy and pandas handle numerical arrays and DataFrames.
  • statsmodels.formula.api provides the mixed-effects model you will fit later.
  • enable_iterative_imputer unlocks scikit-learn's experimental IterativeImputer class. Without this line, the import on the next line fails.
  • BayesianRidge is the regression estimator that supports sample_posterior=True, which is required for proper Multiple Imputation.
  • You won't see a change yet. The next chunk adds the first analysis function.
  • Below the imports, add the complete_case_analysis function:
def complete_case_analysis(df_observed):
    # Identify patients who have a pain score at every visit
    complete_patients = df_observed.groupby("patient_id")["pain_score"].apply(
        lambda x: x.notna().all()
    )
    # Keep only those patient IDs
    complete_ids = complete_patients[complete_patients].index
    df_cc = df_observed[df_observed["patient_id"].isin(complete_ids)].copy()
    return df_cc

What does this code do?

The function groups the dataset by patient_id and checks whether every pain_score in that group is non-null. Only patients who pass this check are retained.

This introduces selection bias. Since patients with higher pain dropped out, the remaining sample over-represents healthier patients in both arms.

  • Save impute_and_analyze.py.
  • Verify the function loads by running:
python -c "from impute_and_analyze import complete_case_analysis; print('Complete Case function loaded')"

You should see Complete Case function loaded printed in your terminal.

Seeing an ImportError?

  • Make sure your virtual environment is still activated. You should see (venv) in your terminal prompt.
  • Check that the from sklearn.experimental import enable_iterative_imputer line appears BEFORE the from sklearn.impute import IterativeImputer line.
  • help me debug an ImportError in impute_and_analyze.py
Add Last Observation Carried Forward

LOCF fills each missing value with the patient's last recorded pain score. If a patient drops out at week 4, their week-4 score gets copied to week 8 and week 12. This assumes pain stays frozen after dropout, which is biologically implausible in a pain trial where untreated symptoms typically worsen.

  • Below complete_case_analysis, add the locf_imputation function:
def locf_imputation(df_observed):
    df_locf = df_observed.copy()
    # Sort so visits are in order before forward-filling
    df_locf = df_locf.sort_values(["patient_id", "visit"])
    # Carry the last non-null pain score forward within each patient
    df_locf["pain_score"] = df_locf.groupby("patient_id")["pain_score"].ffill()
    return df_locf

What does this code do?

The .ffill() method (forward fill) replaces each NaN with the most recent non-null value within the same patient group. Sorting by visit first ensures the fill goes forward in time, not backward.

The result is a fully populated dataset where every patient appears to have stayed at their last observed pain level. This typically attenuates the treatment effect toward zero because it under-estimates worsening in both arms.

  • Save impute_and_analyze.py.
  • Verify by running:
python -c "from impute_and_analyze import locf_imputation; print('LOCF function loaded')"

You should see LOCF function loaded printed.

Seeing a SyntaxError?

  • Check that the function is indented consistently (4 spaces per level, no mixed tabs).
  • Make sure there is no stray text between complete_case_analysis and locf_imputation.
  • help me fix a SyntaxError in my LOCF function
Add Multiple Imputation with analysis helpers

Multiple Imputation takes a fundamentally different approach. Instead of picking one fill value, it generates multiple plausible datasets by modeling the relationships between all observed variables. Each imputed dataset is slightly different, reflecting the uncertainty about what the missing values could have been.

scikit-learn's IterativeImputer implements this through chained equations. You will run it 20 times with different random seeds and sample_posterior=True to draw from the posterior predictive distribution each time.

  • Below locf_imputation, add the multiple_imputation function:
def multiple_imputation(df_observed, m=20, seed=42):
    # Pivot to wide format: one row per patient, one column per visit
    df_wide = df_observed.pivot(index="patient_id", columns="visit", values="pain_score")
    df_wide.columns = [f"pain_visit_{c}" for c in df_wide.columns]

    # Preserve treatment assignment for reshaping later
    treatment_map = df_observed.drop_duplicates("patient_id").set_index("patient_id")["treatment"]
    df_wide["treatment"] = treatment_map

    # Extract only pain columns for imputation
    pain_cols = [c for c in df_wide.columns if c.startswith("pain_visit_")]
    impute_data = df_wide[pain_cols].values

    imputed_datasets = []
    for i in range(m):
        # Each iteration uses a different random seed for variability
        imputer = IterativeImputer(
            estimator=BayesianRidge(),
            sample_posterior=True,
            max_iter=10,
            random_state=seed + i
        )
        imputed_values = imputer.fit_transform(impute_data)

        # Rebuild wide DataFrame with imputed values
        df_imp_wide = pd.DataFrame(imputed_values, columns=pain_cols, index=df_wide.index)
        df_imp_wide["treatment"] = df_wide["treatment"]

        # Melt back to long format for modeling
        df_imp_long = df_imp_wide.reset_index().melt(
            id_vars=["patient_id", "treatment"],
            value_vars=pain_cols,
            var_name="visit_col",
            value_name="pain_score"
        )
        df_imp_long["visit"] = df_imp_long["visit_col"].str.extract(r"(\d+)").astype(int)
        df_imp_long = df_imp_long.drop(columns=["visit_col"])

        imputed_datasets.append(df_imp_long)

    return imputed_datasets

What does this code do?

  • The data is pivoted to wide format (one row per patient) because IterativeImputer expects a 2D array where each column is a feature.
  • sample_posterior=True means each imputation draws from the Bayesian posterior predictive distribution, not just the point estimate. This is what makes each of the 20 datasets different.
  • BayesianRidge is required as the estimator because it supports return_std in its predict method, which is needed for posterior sampling.
  • After imputation, each dataset is melted back to long format (one row per patient-visit) so it matches the structure expected by the analysis model.

Why 20 imputations?

Modern guidance suggests m=20 or more imputations to stabilize variance estimates. The classic recommendation of m=5 is now considered too few for proper confidence interval coverage under moderate missingness.

  • You won't see a change yet. The MI function needs analysis helpers to produce useful results. You will add those next.

To analyze the imputed datasets, you need two helper functions. The first fits a mixed-effects model. The second pools results across imputations using Rubin's Rules.

  • Below multiple_imputation, add the fit_mmrm function:
def fit_mmrm(df):
    # Prepare categorical columns for the model formula
    df = df.copy()
    df["treatment"] = df["treatment"].astype("category")
    df["visit"] = df["visit"].astype("category")

    # MMRM: pain ~ treatment * visit with random intercept per patient
    model = smf.mixedlm(
        "pain_score ~ C(treatment) * C(visit)",
        data=df,
        groups=df["patient_id"],
        re_formula="~1"
    )
    result = model.fit(reml=True, method=["lbfgs"])

    # Extract treatment effect at the final visit (visit 3)
    params = result.params
    se_params = result.bse

    # Main treatment effect
    trt_effect = params.get("C(treatment)[T.1]", 0)

    # Treatment x Visit=3 interaction adds to the main effect
    interaction_key = "C(treatment)[T.1]:C(visit)[T.3]"
    trt_visit3_interaction = params.get(interaction_key, 0)

    # Total effect at visit 3 = main + interaction
    total_effect = trt_effect + trt_visit3_interaction

    # Approximate SE using variance of sum
    se_trt = se_params.get("C(treatment)[T.1]", 0)
    se_interaction = se_params.get(interaction_key, 0)
    total_se = np.sqrt(se_trt**2 + se_interaction**2)

    return total_effect, total_se

What does this code do?

This function fits a Mixed-Effects Model for Repeated Measures (MMRM). The formula pain_score ~ C(treatment) * C(visit) models pain as a function of treatment group, time point, and their interaction.

The treatment effect at the final visit is the sum of the main treatment coefficient and the treatment-by-visit-3 interaction. The function returns both the estimate and its standard error.

  • You won't see a change yet. One more helper function to go, then you will verify everything loads.
  • Below fit_mmrm, add the pool_mi_results function:
def pool_mi_results(estimates, standard_errors):
    m = len(estimates)
    estimates = np.array(estimates)
    standard_errors = np.array(standard_errors)

    # Pooled estimate: mean across all imputations
    pooled_estimate = np.mean(estimates)

    # Within-imputation variance: average of squared SEs
    within_var = np.mean(standard_errors**2)

    # Between-imputation variance: variance of the estimates themselves
    between_var = np.var(estimates, ddof=1)

    # Total variance combines both sources (Rubin's Rules)
    total_var = within_var + (1 + 1/m) * between_var

    pooled_se = np.sqrt(total_var)

    return pooled_estimate, pooled_se

What are Rubin's Rules?

When you fit a model on each of the 20 imputed datasets, you get 20 slightly different treatment effect estimates. Rubin's Rules combine them into a single pooled result.

The total variance has two components. Within-imputation variance captures the usual sampling uncertainty. Between-imputation variance captures the extra uncertainty from not knowing the true missing values. The (1 + 1/m) correction adjusts for using a finite number of imputations.

  • Save impute_and_analyze.py.
  • Verify all functions load by running:
python -c "from impute_and_analyze import complete_case_analysis, locf_imputation, multiple_imputation, fit_mmrm, pool_mi_results; print('All 5 functions loaded successfully')"

You should see All 5 functions loaded successfully in your terminal.

Import fails or hangs?

  • If you see an error about enable_iterative_imputer, make sure that import line comes before from sklearn.impute import IterativeImputer. Order matters.
  • If you see an IndentationError, check that all function bodies use exactly 4 spaces per indent level.
  • If it hangs indefinitely, press Ctrl+C to cancel and check for infinite loops or misplaced code outside a function.
  • help me debug why my impute_and_analyze.py won't import

✔️ Awesome, I've got everything!

Great. Double check you have saved impute_and_analyze.py with all five functions defined.

ⓧ I'd like to double check the full code

Here is the complete impute_and_analyze.py file at this point:

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.impute import IterativeImputer
from sklearn.linear_model import BayesianRidge


def complete_case_analysis(df_observed):
    # Identify patients who have a pain score at every visit
    complete_patients = df_observed.groupby("patient_id")["pain_score"].apply(
        lambda x: x.notna().all()
    )
    # Keep only those patient IDs
    complete_ids = complete_patients[complete_patients].index
    df_cc = df_observed[df_observed["patient_id"].isin(complete_ids)].copy()
    return df_cc


def locf_imputation(df_observed):
    df_locf = df_observed.copy()
    # Sort so visits are in order before forward-filling
    df_locf = df_locf.sort_values(["patient_id", "visit"])
    # Carry the last non-null pain score forward within each patient
    df_locf["pain_score"] = df_locf.groupby("patient_id")["pain_score"].ffill()
    return df_locf


def multiple_imputation(df_observed, m=20, seed=42):
    # Pivot to wide format: one row per patient, one column per visit
    df_wide = df_observed.pivot(index="patient_id", columns="visit", values="pain_score")
    df_wide.columns = [f"pain_visit_{c}" for c in df_wide.columns]

    # Preserve treatment assignment for reshaping later
    treatment_map = df_observed.drop_duplicates("patient_id").set_index("patient_id")["treatment"]
    df_wide["treatment"] = treatment_map

    # Extract only pain columns for imputation
    pain_cols = [c for c in df_wide.columns if c.startswith("pain_visit_")]
    impute_data = df_wide[pain_cols].values

    imputed_datasets = []
    for i in range(m):
        # Each iteration uses a different random seed for variability
        imputer = IterativeImputer(
            estimator=BayesianRidge(),
            sample_posterior=True,
            max_iter=10,
            random_state=seed + i
        )
        imputed_values = imputer.fit_transform(impute_data)

        # Rebuild wide DataFrame with imputed values
        df_imp_wide = pd.DataFrame(imputed_values, columns=pain_cols, index=df_wide.index)
        df_imp_wide["treatment"] = df_wide["treatment"]

        # Melt back to long format for modeling
        df_imp_long = df_imp_wide.reset_index().melt(
            id_vars=["patient_id", "treatment"],
            value_vars=pain_cols,
            var_name="visit_col",
            value_name="pain_score"
        )
        df_imp_long["visit"] = df_imp_long["visit_col"].str.extract(r"(\d+)").astype(int)
        df_imp_long = df_imp_long.drop(columns=["visit_col"])

        imputed_datasets.append(df_imp_long)

    return imputed_datasets


def fit_mmrm(df):
    # Prepare categorical columns for the model formula
    df = df.copy()
    df["treatment"] = df["treatment"].astype("category")
    df["visit"] = df["visit"].astype("category")

    # MMRM: pain ~ treatment * visit with random intercept per patient
    model = smf.mixedlm(
        "pain_score ~ C(treatment) * C(visit)",
        data=df,
        groups=df["patient_id"],
        re_formula="~1"
    )
    result = model.fit(reml=True, method=["lbfgs"])

    # Extract treatment effect at the final visit (visit 3)
    params = result.params
    se_params = result.bse

    # Main treatment effect
    trt_effect = params.get("C(treatment)[T.1]", 0)

    # Treatment x Visit=3 interaction adds to the main effect
    interaction_key = "C(treatment)[T.1]:C(visit)[T.3]"
    trt_visit3_interaction = params.get(interaction_key, 0)

    # Total effect at visit 3 = main + interaction
    total_effect = trt_effect + trt_visit3_interaction

    # Approximate SE using variance of sum
    se_trt = se_params.get("C(treatment)[T.1]", 0)
    se_interaction = se_params.get(interaction_key, 0)
    total_se = np.sqrt(se_trt**2 + se_interaction**2)

    return total_effect, total_se


def pool_mi_results(estimates, standard_errors):
    m = len(estimates)
    estimates = np.array(estimates)
    standard_errors = np.array(standard_errors)

    # Pooled estimate: mean across all imputations
    pooled_estimate = np.mean(estimates)

    # Within-imputation variance: average of squared SEs
    within_var = np.mean(standard_errors**2)

    # Between-imputation variance: variance of the estimates themselves
    between_var = np.var(estimates, ddof=1)

    # Total variance combines both sources (Rubin's Rules)
    total_var = within_var + (1 + 1/m) * between_var

    pooled_se = np.sqrt(total_var)

    return pooled_estimate, pooled_se

Your three imputation strategies are ready. Next, you will run all three on the observed dataset, fit the MMRM on each result, and build a comparison table showing exactly how much bias each method introduces.

Fit MMRM and Compare Bias

You now have three versions of your trial data: a complete-case subset, an LOCF-imputed dataset, and 20 multiply-imputed datasets. The question is: which one recovers the true treatment effect of -1.5?

In this step, you will fit a Mixed-Effects Model for Repeated Measures (MMRM) to each dataset, then build a comparison table that quantifies how much bias each imputation method introduces. This is the payoff of the entire simulation. Because you know the true effect, you can measure exactly how far each method drifts from reality.

In this step, get ready to:
  • Fit the MMRM on complete-case and LOCF datasets to obtain treatment effect estimates.
  • Pool Multiple Imputation results across 20 datasets using Rubin's Rules.
  • Build a comparison table showing bias, standard errors, and confidence interval coverage.
Fit the MMRM on complete-case and LOCF datasets

The fit_mmrm() function you defined in the previous step handles model fitting. Now you need a __main__ block that loads your CSV data, runs each imputation strategy, and feeds the results into that function.

  • In impute_and_analyze.py, scroll to the bottom of the file (after the pool_mi_results function).
  • Add the following __main__ block to load the data and fit the MMRM on the complete-case and LOCF datasets:
if __name__ == "__main__":
    # Load data
    df_complete = pd.read_csv("trial_complete.csv")
    df_observed = pd.read_csv("trial_observed.csv")

    TRUE_EFFECT = -1.5

    print("=" * 70)
    print("MISSING DATA IMPUTATION COMPARISON STUDY")
    print("=" * 70)

    # --- 1. Complete Case Analysis ---
    print("\n[1/3] Running Complete Case Analysis...")
    df_cc = complete_case_analysis(df_observed)
    n_complete = df_cc["patient_id"].nunique()
    print(f"    Patients retained: {n_complete}/200 ({100*n_complete/200:.0f}%)")
    cc_estimate, cc_se = fit_mmrm(df_cc)
    print(f"    Treatment effect estimate: {cc_estimate:.3f} (SE={cc_se:.3f})")

    # --- 2. LOCF ---
    print("\n[2/3] Running Last Observation Carried Forward...")
    df_locf = locf_imputation(df_observed)
    locf_estimate, locf_se = fit_mmrm(df_locf)
    print(f"    Treatment effect estimate: {locf_estimate:.3f} (SE={locf_se:.3f})")

What does this code do?

This block loads trial_observed.csv (which has NaN dropout values) and runs it through two imputation strategies, fitting the same MMRM model to each.

  • TRUE_EFFECT = -1.5 stores the known ground truth so you can measure bias later.
  • complete_case_analysis() drops patients with any missing observations, then fit_mmrm() estimates the treatment effect on that reduced sample.
  • locf_imputation() carries each patient's last observed score forward, then fit_mmrm() estimates the treatment effect on the filled dataset.

What is the MMRM formula doing?

The model formula pain_score ~ C(treatment) * C(visit) fits main effects for treatment and visit, plus their interaction. The interaction at the final visit (C(treatment)[T.1]:C(visit)[T.3]) captures how the treatment effect evolves over time.

The total treatment effect at the final visit is the sum of the main treatment coefficient and this interaction term. A random intercept (re_formula="~1") accounts for each patient's individual baseline pain level.

Pool Multiple Imputation results with Rubin's Rules

For Multiple Imputation, you cannot just fit one model. You fit the MMRM on each of the 20 imputed datasets, then combine the results using Rubin's Rules to get a single pooled estimate and standard error that properly reflects imputation uncertainty.

  • Directly below the LOCF section you just added, paste the Multiple Imputation fitting and pooling code:
    # --- 3. Multiple Imputation ---
    print("\n[3/3] Running Multiple Imputation (m=20)...")
    imputed_datasets = multiple_imputation(df_observed, m=20)
    mi_estimates = []
    mi_ses = []
    for i, df_imp in enumerate(imputed_datasets):
        est, se = fit_mmrm(df_imp)
        mi_estimates.append(est)
        mi_ses.append(se)
    mi_pooled_estimate, mi_pooled_se = pool_mi_results(mi_estimates, mi_ses)
    print(f"    Pooled treatment effect estimate: {mi_pooled_estimate:.3f} (SE={mi_pooled_se:.3f})")

What does this code do?

This loop fits the identical MMRM model to each of the 20 imputed datasets, collecting 20 treatment effect estimates and 20 standard errors.

  • pool_mi_results() combines them using Rubin's Rules: the pooled estimate is the mean of all 20, and the total variance combines within-imputation variance (average of squared SEs) with between-imputation variance (how much the estimates scatter across imputations).
  • The factor (1 + 1/m) in Rubin's formula inflates the between-imputation variance to account for using a finite number of imputations rather than infinitely many.

Seeing a ConvergenceWarning?

statsmodels may print a warning that the MLE may be on the boundary of the parameter space. This happens when the random intercept variance is small. The results are still usable as long as you see finite estimates and standard errors in the output.

help me understand the ConvergenceWarning from statsmodels MixedLM

Build the comparison table and interpretation

Now for the most revealing part of the study. You will print a formatted table comparing each method's estimate against the known truth, showing the bias, standard error, and whether the 95% confidence interval covers the true value.

  • Directly below the MI pooling code, add the comparison table and interpretation:
    # --- Comparison Table ---
    print("\n" + "=" * 70)
    print("COMPARATIVE RESULTS")
    print("=" * 70)
    print(f"\n{'Method':<25} {'Estimate':<12} {'Bias':<10} {'SE':<10} {'95% CI':<20}")
    print("-" * 70)

    methods = [
        ("Complete Case", cc_estimate, cc_se),
        ("LOCF", locf_estimate, locf_se),
        ("Multiple Imputation", mi_pooled_estimate, mi_pooled_se),
    ]

    for name, est, se in methods:
        bias = est - TRUE_EFFECT
        ci_low = est - 1.96 * se
        ci_high = est + 1.96 * se
        covers = "Yes" if ci_low <= TRUE_EFFECT <= ci_high else "No"
        print(f"{name:<25} {est:<12.3f} {bias:<10.3f} {se:<10.3f} [{ci_low:.2f}, {ci_high:.2f}] Covers: {covers}")

    # Also show the "truth" from complete data
    print(f"\n{'True (full data)':<25} {TRUE_EFFECT:<12.3f} {'0.000':<10} {'---':<10}")

    # --- Interpretation ---
    print("\n" + "=" * 70)
    print("INTERPRETATION")
    print("=" * 70)
    print("""
Why Multiple Imputation outperforms legacy methods under MAR:

1. COMPLETE CASE ANALYSIS drops patients who experienced higher pain (since
   they dropped out more). This creates selection bias: the remaining sample
   over-represents patients who were doing well, potentially in BOTH arms,
   distorting the treatment contrast.

2. LOCF assumes patients' pain stays constant after dropout. In a pain trial,
   treatment-arm patients who drop out likely had rising pain (treatment
   failure), but LOCF freezes them at their earlier (lower) score. This
   typically attenuates the treatment effect toward zero because it under-
   estimates the worsening in both arms, but the effect is differential.

3. MULTIPLE IMPUTATION models the relationships between observed variables
   (treatment, visit, observed pain trajectory) to generate plausible values
   for the missing data. By creating multiple datasets and pooling, it:
   - Properly accounts for the uncertainty due to missing data
   - Uses all available information to predict missing values
   - Produces unbiased estimates when the imputation model is correctly specified
   - Yields valid standard errors and confidence intervals via Rubin's Rules

Under MAR (where missingness depends on observed data), MI is theoretically
unbiased when the imputation model includes the variables that drive missingness.
""")

What does the comparison table reveal?

The "Bias" column shows the difference between each method's estimate and the true effect of -1.5. Closer to zero means less distortion.

  • Complete Case typically shows moderate bias because it discards patients non-randomly (those with higher pain dropped out more).
  • LOCF typically attenuates the treatment effect toward zero because freezing scores at the last observation understates ongoing changes.
  • Multiple Imputation should show the smallest bias because it models the relationships between observed variables and properly accounts for uncertainty.
  • The "Covers" column tells you whether the 95% CI contains the true value. Valid methods should cover about 95% of the time.
  • Save impute_and_analyze.py.
  • Run the analysis script in your terminal:
python impute_and_analyze.py
  • You should see a formatted table with three rows (Complete Case, LOCF, Multiple Imputation) showing each method's estimate, bias, SE, and 95% CI. Multiple Imputation should have the smallest absolute bias.

Script fails with a FileNotFoundError?

Make sure trial_complete.csv and trial_observed.csv exist in your missing-data-sim folder. These are generated by running python simulate_trial.py from the previous step. If you skipped that, run it now before retrying.

my impute_and_analyze.py script can't find the CSV files

✔️ Awesome, I've got everything!

Great. Double check you have saved impute_and_analyze.py and that the comparison table printed successfully with all three methods.

ⓧ I'd like to double check the full code

Here is the complete impute_and_analyze.py file:

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.impute import IterativeImputer
from sklearn.linear_model import BayesianRidge


def complete_case_analysis(df_observed):
    complete_patients = df_observed.groupby("patient_id")["pain_score"].apply(
        lambda x: x.notna().all()
    )
    complete_ids = complete_patients[complete_patients].index
    df_cc = df_observed[df_observed["patient_id"].isin(complete_ids)].copy()
    return df_cc


def locf_imputation(df_observed):
    df_locf = df_observed.copy()
    df_locf = df_locf.sort_values(["patient_id", "visit"])
    df_locf["pain_score"] = df_locf.groupby("patient_id")["pain_score"].ffill()
    return df_locf


def multiple_imputation(df_observed, m=20, seed=42):
    df_wide = df_observed.pivot(index="patient_id", columns="visit", values="pain_score")
    df_wide.columns = [f"pain_visit_{c}" for c in df_wide.columns]
    treatment_map = df_observed.drop_duplicates("patient_id").set_index("patient_id")["treatment"]
    df_wide["treatment"] = treatment_map
    pain_cols = [c for c in df_wide.columns if c.startswith("pain_visit_")]
    impute_data = df_wide[pain_cols].values
    imputed_datasets = []
    for i in range(m):
        imputer = IterativeImputer(
            estimator=BayesianRidge(),
            sample_posterior=True,
            max_iter=10,
            random_state=seed + i
        )
        imputed_values = imputer.fit_transform(impute_data)
        df_imp_wide = pd.DataFrame(imputed_values, columns=pain_cols, index=df_wide.index)
        df_imp_wide["treatment"] = df_wide["treatment"]
        df_imp_long = df_imp_wide.reset_index().melt(
            id_vars=["patient_id", "treatment"],
            value_vars=pain_cols,
            var_name="visit_col",
            value_name="pain_score"
        )
        df_imp_long["visit"] = df_imp_long["visit_col"].str.extract(r"(\d+)").astype(int)
        df_imp_long = df_imp_long.drop(columns=["visit_col"])
        imputed_datasets.append(df_imp_long)
    return imputed_datasets


def fit_mmrm(df):
    df = df.copy()
    df["treatment"] = df["treatment"].astype("category")
    df["visit"] = df["visit"].astype("category")
    model = smf.mixedlm(
        "pain_score ~ C(treatment) * C(visit)",
        data=df,
        groups=df["patient_id"],
        re_formula="~1"
    )
    result = model.fit(reml=True, method=["lbfgs"])
    params = result.params
    se_params = result.bse
    trt_effect = params.get("C(treatment)[T.1]", 0)
    interaction_key = "C(treatment)[T.1]:C(visit)[T.3]"
    trt_visit3_interaction = params.get(interaction_key, 0)
    total_effect = trt_effect + trt_visit3_interaction
    se_trt = se_params.get("C(treatment)[T.1]", 0)
    se_interaction = se_params.get(interaction_key, 0)
    total_se = np.sqrt(se_trt**2 + se_interaction**2)
    return total_effect, total_se


def pool_mi_results(estimates, standard_errors):
    m = len(estimates)
    estimates = np.array(estimates)
    standard_errors = np.array(standard_errors)
    pooled_estimate = np.mean(estimates)
    within_var = np.mean(standard_errors**2)
    between_var = np.var(estimates, ddof=1)
    total_var = within_var + (1 + 1/m) * between_var
    pooled_se = np.sqrt(total_var)
    return pooled_estimate, pooled_se


if __name__ == "__main__":
    df_complete = pd.read_csv("trial_complete.csv")
    df_observed = pd.read_csv("trial_observed.csv")
    TRUE_EFFECT = -1.5
    print("=" * 70)
    print("MISSING DATA IMPUTATION COMPARISON STUDY")
    print("=" * 70)
    print("\n[1/3] Running Complete Case Analysis...")
    df_cc = complete_case_analysis(df_observed)
    n_complete = df_cc["patient_id"].nunique()
    print(f"    Patients retained: {n_complete}/200 ({100*n_complete/200:.0f}%)")
    cc_estimate, cc_se = fit_mmrm(df_cc)
    print(f"    Treatment effect estimate: {cc_estimate:.3f} (SE={cc_se:.3f})")
    print("\n[2/3] Running Last Observation Carried Forward...")
    df_locf = locf_imputation(df_observed)
    locf_estimate, locf_se = fit_mmrm(df_locf)
    print(f"    Treatment effect estimate: {locf_estimate:.3f} (SE={locf_se:.3f})")
    print("\n[3/3] Running Multiple Imputation (m=20)...")
    imputed_datasets = multiple_imputation(df_observed, m=20)
    mi_estimates = []
    mi_ses = []
    for i, df_imp in enumerate(imputed_datasets):
        est, se = fit_mmrm(df_imp)
        mi_estimates.append(est)
        mi_ses.append(se)
    mi_pooled_estimate, mi_pooled_se = pool_mi_results(mi_estimates, mi_ses)
    print(f"    Pooled treatment effect estimate: {mi_pooled_estimate:.3f} (SE={mi_pooled_se:.3f})")
    print("\n" + "=" * 70)
    print("COMPARATIVE RESULTS")
    print("=" * 70)
    print(f"\n{'Method':<25} {'Estimate':<12} {'Bias':<10} {'SE':<10} {'95% CI':<20}")
    print("-" * 70)
    methods = [
        ("Complete Case", cc_estimate, cc_se),
        ("LOCF", locf_estimate, locf_se),
        ("Multiple Imputation", mi_pooled_estimate, mi_pooled_se),
    ]
    for name, est, se in methods:
        bias = est - TRUE_EFFECT
        ci_low = est - 1.96 * se
        ci_high = est + 1.96 * se
        covers = "Yes" if ci_low <= TRUE_EFFECT <= ci_high else "No"
        print(f"{name:<25} {est:<12.3f} {bias:<10.3f} {se:<10.3f} [{ci_low:.2f}, {ci_high:.2f}] Covers: {covers}")
    print(f"\n{'True (full data)':<25} {TRUE_EFFECT:<12.3f} {'0.000':<10} {'---':<10}")
    print("\n" + "=" * 70)
    print("INTERPRETATION")
    print("=" * 70)
    print("""
Why Multiple Imputation outperforms legacy methods under MAR:

1. COMPLETE CASE ANALYSIS drops patients who experienced higher pain (since
   they dropped out more). This creates selection bias: the remaining sample
   over-represents patients who were doing well, potentially in BOTH arms,
   distorting the treatment contrast.

2. LOCF assumes patients' pain stays constant after dropout. In a pain trial,
   treatment-arm patients who drop out likely had rising pain (treatment
   failure), but LOCF freezes them at their earlier (lower) score. This
   typically attenuates the treatment effect toward zero because it under-
   estimates the worsening in both arms, but the effect is differential.

3. MULTIPLE IMPUTATION models the relationships between observed variables
   (treatment, visit, observed pain trajectory) to generate plausible values
   for the missing data. By creating multiple datasets and pooling, it:
   - Properly accounts for the uncertainty due to missing data
   - Uses all available information to predict missing values
   - Produces unbiased estimates when the imputation model is correctly specified
   - Yields valid standard errors and confidence intervals via Rubin's Rules

Under MAR (where missingness depends on observed data), MI is theoretically
unbiased when the imputation model includes the variables that drive missingness.
""")

You have quantitative proof that imputation method choice directly impacts clinical trial conclusions. Multiple Imputation preserves the true treatment effect while legacy methods introduce systematic distortion. Next up, the Secret Mission challenges you to run this entire pipeline 500 times and visualize the sampling distributions.

Secret mission

Monte Carlo Simulation Study

Your single-run comparison is insightful, but how confident can you be that the results hold across different random samples? In this extension, you will wrap your entire pipeline in a Monte Carlo study that runs 500 replications and plots the sampling distribution of bias under each imputation method.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs. Decide whether to keep your files for future reference, or delete them to free up disk space.

Resources you used:

  • The missing-data-sim project folder on your Desktop (contains Python scripts, virtual environment, and generated data files).

Keep everything running

No action needed. Choose this if you want to continue experimenting with different MAR parameters, sample sizes, or imputation models.

  • Your simulation scripts, generated CSVs, and Monte Carlo plot remain in the missing-data-sim folder on your Desktop.
  • The virtual environment venv contains all installed packages. Reactivate it any time to pick up where you left off.

Pause - I'll come back to this later

Since everything is local files, there is nothing to pause. Your project will sit idle on disk until you return.

  • If you have a terminal still running a script, close it to free up memory.
  • Deactivate your virtual environment by running:
deactivate

When you come back, reactivate with the same activation command you used in Step 1.

Delete - I don't want to use this again

Remove all project files, the virtual environment, and generated data. This frees up disk space completely.

  • Deactivate your virtual environment first (if it is still active):
deactivate
  • Delete the entire project folder by running this from your Desktop:

macOS

cd ~/Desktop
rm -rf missing-data-sim

Windows

cd ~\Desktop
rmdir /s /q missing-data-sim
  • Confirm the folder is gone by running ls (macOS) or dir (Windows) on your Desktop. You should no longer see missing-data-sim in the output.

Nice Work!

Nice Work!

You just built a rigorous simulation study that quantifies exactly how much different missing data strategies distort clinical trial conclusions. This is the kind of evidence-based methodology knowledge that biostatistics positions and regulatory submissions demand.

You've learned how to:

  • Simulate a longitudinal clinical trial dataset with correlated repeated measures and inject realistic MAR dropout where patients with higher pain drop out more frequently.
  • Apply three imputation strategies (Complete Case Analysis, LOCF, and Multiple Imputation) and fit MMRM models to compare treatment effect bias against the known truth.
  • Prove quantitatively that Multiple Imputation with Rubin's Rules produces the least biased estimates under MAR, while LOCF attenuates and Complete Case distorts the treatment effect.
  • Built a Monte Carlo simulation study running 500 replications to plot sampling distributions and compute coverage probabilities for each method.

Ready to quiz yourself?