Author: Michele Danilo Pierri

Michele D. Pierri is a cardiac surgeon and cardiovascular physiopathology researcher with a strong interest in artificial intelligence, medical data science, clinical decision support, and digital health. His work focuses on the intersection between medicine, technology, and computational methods, with the aim of translating complex biomedical concepts into clear, practical, and clinically meaningful insights.
Surreal historical medical illustration of a woman holding an anatomical heart, surrounded by a hospital bed, childbirth imagery, anatomical charts, medicine bottles, flowers, and distant factory smokestacks in a muted vintage palette.

Frida Kahlo: The Anatomy of Suffering

Posted on by Michele Danilo Pierri

The Anatomy of Suffering: Frida Kahlo through the Lens of Medical History

Introduction: Why Frida Kahlo is the Patron Saint of Narrative Medicine

In the intersection of art and clinical science, few figures loom as large as Frida Kahlo (1907–1954). While often celebrated for her surrealist aesthetics and her role in Mexican Modernism, Kahlo’s work serves as one of the most significant visual archives in the history of medicine. For the medical historian, her canvases are not merely “paintings”—they are clinical case studies, surgical records, and a masterclass in the phenomenology of chronic pain.

This article explores Kahlo’s unique position as both a frustrated medical student and a lifelong “professional patient,” analyzing how her knowledge of anatomy and the medical landscape of the early 20th century shaped a new iconography of the human body.

I. The Education of a Physician Interrupted

To understand the medical accuracy in Kahlo’s work, one must look at her youth. In 1922, Frida was one of only 35 girls accepted into the Escuela Nacional Preparatoria, Mexico’s most prestigious secondary school. Her goal was clear: she wanted to become a doctor.

At the school, she focused on biology, botany, and anatomy. Her sketches from this period show a precocious understanding of the skeletal system and organic structures. However, the trajectory of her life was irrevocably altered on September 17, 1925. A tram accident resulted in a steel handrail piercing her abdomen and exit through her vagina, causing:

  • Triple fractures of the spine.
  • Fractures of the clavicle, ribs, and pelvis.
  • Eleven fractures in her right leg and a crushed foot.
  • Dislocation of the shoulder.

From that moment on, the woman who wished to study medicine became its subject. Her subsequent 32 surgeries and years of confinement transformed her bedroom into a laboratory where the subject of study was her own deteriorating physiology.

II. Medical Context: Orthopedics and Surgery in Post-Revolutionary Mexico

Kahlo’s medical journey coincided with a transformative era in global medicine. The 1920s through the 1950s saw the transition from traditional surgical practices to highly mechanized, industrial medicine.

  1. The Rise of Radiology: X-rays were becoming a primary diagnostic tool. Kahlo was fascinated by her own radiographs, often using the “transparency” of the body in her paintings (showing internal organs through the skin) much like a clinical X-ray.
  2. Orthopedic Limitations: This was the “age of the corset.” Without the advanced spinal fusion techniques or biocompatible titanium implants we have today, patients with spinal trauma were subjected to months of immobilization in plaster casts (escayolas) or heavy steel braces.
  3. The Industrialization of Care: During her time at the Henry Ford Hospital in Detroit, Kahlo witnessed the American medical model—efficient but cold and mechanical. This contrast between the “mechanical” and the “organic” became a central theme in her medical paintings.

III. Clinical Analysis of Key Works

1. The Broken Column (1944): The Iconography of Spinal Trauma

The Broken Column is perhaps the most significant artistic representation of chronic neuropathic pain in history.

  • Medical Content: At the center of the painting, Kahlo replaces her spine with a crumbling Ionic column. This choice is medically poignant: a column provides structural integrity but, once cracked, threatens the collapse of the entire edifice. The column is broken in several places, corresponding to the locations of her actual vertebral fractures.
  • Technique and Meaning: She depicts her torso held together by a white orthopedic corset. The canvas is an anatomical cross-section; the skin is “zipped” open to reveal the structural failure within.
  • Clinical Significance: The nails driven into her skin represent the allodynia and constant sensory irritation associated with spinal nerve damage. For a medical professional, this painting illustrates the patient’s internal perception of their anatomy—not as a cohesive whole, but as a fragmented structure failing under gravity.

The Broken Column - Frida Kahlo
The Broken Column, 1944 — Frida Kahlo
Source: WikiArt — Museo Dolores Olmedo — Used for educational/editorial purposes

2. Henry Ford Hospital (1932): Obstetric Trauma and Industrial Coldness

Painted after a traumatic miscarriage in Detroit, this work is a “lithographic” clinical report.

  • Medical Content: Kahlo depicts herself on a bed floating in a desolate industrial landscape. Six umbilical-like veins connect her to objects of clinical significance:
    • The Male Fetus: A representation of “Dieguito,” her lost son.
    • The Pelvic Bone: A medically accurate rendering of the female pelvis, showing the deformities that prevented a natural birth.
    • The Snails: Symbolic of the “slow,” agonizing nature of the failed pregnancy.
    • The Autoclave: A piece of hospital equipment used for sterilization, representing the dehumanizing, mechanical nature of modern obstetrics.
  • Clinical Significance: This work is a rare historical document of obstetric grief and the limitations of early 20th-century gynecology in treating patients with pelvic trauma.
Henry Ford Hospital, 1932 — Frida Kahlo
Source: Wikimedia Commons — Dolores Olmedo Collection — Used for educational/editorial purposes

3. The Two Fridas (1939): Anatomy and Haemostasis

This double self-portrait is a masterclass in cardiovascular symbolism.

  • Medical Content: The two versions of Frida are linked by a single vein connecting two hearts. The “European” Frida on the left has a heart that is surgically “opened”—an anatomical dissection showing the internal chambers.
  • Technique: She holds a hemostat (surgical clamp) in her hand, attempting to stop the flow of blood from a severed vein. This shows her sophisticated knowledge of haemostasis.
  • Meaning: The blood dripping onto her white dress signifies a “hemorrhage of the soul,” but the use of a clinical tool like the hemostat suggests the patient’s desperate attempt to self-medicate or “suture” her own emotional wounds.
The Two Fridas, 1939 — Frida Kahlo
Source: Wikimedia Commons — Museo de Arte Moderno, Mexico City — Used for educational/editorial purposes

4. Without Hope (1945): Inanition and Forced Nutrition

In her later years, Kahlo suffered from extreme weight loss and lack of appetite.

  • Medical Content: The painting shows a “forced feeding” apparatus. A large wooden structure (resembling an easel but acting as a funnel) pours a grotesque slurry of animal carcasses and offal into her mouth.
  • Clinical Significance: This addresses the psychological trauma of enteral nutrition and the loss of bodily autonomy. In the history of medicine, this work serves as a reminder that “life-saving” interventions can often feel like violations to the patient.
Without Hope, 1945 — Frida Kahlo
Source: WikiArt — Used for educational/editorial purposes

IV. The Legacy: Kahlo and the “Medical Humanities”

Frida Kahlo’s contribution to the history of medicine goes beyond the documentation of her own ailments. She is a pioneer of Narrative Medicine.

Today, medical schools around the world use Kahlo’s paintings to teach Visual Thinking Strategies (VTS). By observing her work, students learn to:

  1. Identify non-verbal cues of pain: The stoic expression in her paintings vs. the ravaged body.
  2. Understand the “Patient’s Voice”: Recognizing that a medical chart (the “objective” view) is incomplete without the patient’s subjective experience of their “brokenness.”
  3. Historical Empathy: Viewing the evolution of orthopedic and surgical tools through the eyes of the one who had to wear them.

Conclusion: The Artist as Anatomist

Frida Kahlo did not become a doctor in the traditional sense, but she became an anatomist of the human condition. Her work bridged the gap between the sterile, objective world of the clinic and the visceral, subjective world of the sufferer.

For the student of medical history, Kahlo’s oeuvre remains a vital record of early 20th-century trauma surgery and a timeless testament to the resilience of the human spirit in the face of physiological collapse. She reminds us that behind every “broken column” or “severed vein” in a medical textbook, there is a human story that demands to be seen.

Q&A: Understanding Frida Kahlo Through Medicine and Art

Basic Comprehension Questions

Q: What was Frida Kahlo’s original career aspiration before her tragic accident?

A: Frida Kahlo wanted to become a doctor. In 1922, she was one of only 35 girls accepted into the Escuela Nacional Preparatoria, Mexico’s most prestigious secondary school, where she focused on biology, botany, and anatomy in preparation for medical studies.

Q: What injuries did Kahlo sustain in the 1925 tram accident?

A: The accident was catastrophic. A steel handrail pierced her abdomen and exited through her vagina. She suffered triple fractures of the spine, fractures of the clavicle, ribs, and pelvis, eleven fractures in her right leg, a crushed foot, and shoulder dislocation. She underwent 32 surgeries throughout her life as a result.

Q: Why is The Broken Column (1944) considered medically significant?

A: It’s one of the most significant artistic representations of chronic neuropathic pain in history. The crumbling Ionic column replacing her spine corresponds to her actual vertebral fractures, while the nails driven into her skin represent allodynia and constant sensory irritation associated with spinal nerve damage.

Intermediate Analysis Questions

Q: How did Kahlo’s medical education influence her artistic technique?

A: Her anatomical training is evident in her precise rendering of skeletal structures, organs, and medical instruments. She used techniques like “transparency” (showing internal organs through skin, similar to X-rays) and anatomical cross-sections, transforming her canvases into visual medical records.

Q: What does Henry Ford Hospital (1932) reveal about early 20th-century obstetric medicine?

A: The painting documents the limitations of gynecology in treating patients with pelvic trauma. It shows obstetric grief and the dehumanizing, mechanical nature of modern obstetrics through symbols like the autoclave and the medically accurate rendering of her deformed pelvis that prevented natural birth.

Q: What is the cardiovascular symbolism in The Two Fridas (1939)?

A: The painting demonstrates sophisticated knowledge of anatomy and haemostasis. Two versions of Frida are connected by a vein linking two hearts. One heart is surgically “opened” showing internal chambers, and she holds a hemostat (surgical clamp) attempting to stop blood flow—a metaphor for trying to “suture” emotional wounds using clinical tools.

Advanced Critical Thinking Questions

Q: How does Kahlo’s work bridge the gap between objective medical records and subjective patient experience?

A: Kahlo’s paintings provide what medical charts cannot: the phenomenology of suffering. While medical records document fractures and procedures objectively, her art reveals the patient’s internal perception—the body as fragmented, failing, violated. This dual perspective is why she’s considered a pioneer of Narrative Medicine.

Q: What historical medical context shaped Kahlo’s treatment and artistic response?

A: Kahlo lived through a transformative era (1920s-1950s) when medicine was becoming industrialized. She experienced the rise of radiology, the limitations of orthopedic care (plaster casts and steel braces instead of modern spinal fusion), and the contrast between Mexican and American medical models. These experiences of “mechanical” versus “organic” care became central themes in her work.

Q: Why is Kahlo called “the Patron Saint of Narrative Medicine”?

A: Kahlo transformed her medical experience into a visual archive that communicates the lived reality of chronic illness. Medical schools worldwide now use her paintings to teach Visual Thinking Strategies, helping students identify non-verbal pain cues, understand the patient’s voice, and develop historical empathy—core competencies of Narrative Medicine.

Discussion Questions

Q: How does Without Hope (1945) challenge our understanding of “life-saving” medical interventions?

A: The painting depicts forced feeding as grotesque and violating, showing a slurry of animal carcasses being funneled into her mouth. It addresses the psychological trauma of enteral nutrition and loss of bodily autonomy, reminding us that interventions meant to save lives can feel like violations to patients, raising questions about medical paternalism and patient dignity.

Q: What makes Kahlo an “anatomist of the human condition”?

A: While she never became a doctor, Kahlo used her anatomical knowledge to dissect not just physical bodies but the human experience of suffering, trauma, and resilience. She bridged the sterile, objective clinical world with the visceral, subjective world of the patient, creating a unique medical-artistic language that speaks to both clinicians and sufferers.

Q: How might modern medical students benefit from studying Kahlo’s work?

A: Students can learn to see beyond clinical data to the patient’s subjective experience, recognize how historical context shapes medical practice, understand the psychological impact of chronic pain and invasive procedures, and develop empathy by witnessing how one patient documented her decades-long medical journey through art.

 Explore More: Art, Medicine, and the Humanities

A vintage pastoral scene of a young boy repairing a mechanical cart in a rural village, surrounded by children, weathered homes, and warm sepia light.

Turbo Regression

Posted on by Michele Danilo Pierri

How to Put a Turbo on Regression Models

Bootstrap Validation and Cubic Splines for Medical Data

A practical guide for data scientists and clinical researchers working with limited medical datasets

Last updated: December 2025

Author: Michele D. Pierri

Reading time: 15–20 minutes

Introduction — Why regression often disappoints in medicine

Regression models are still the backbone of medical research.

Risk scores, prognostic models, outcome prediction after surgery, ICU mortality, disease progression — all of these rely heavily on regression.

Yet, anyone who has tried to deploy a regression model outside the dataset it was built on knows a frustrating truth:

models that look excellent on paper often perform worse in real patients.

There are two recurring reasons for this failure.

First, models tend to be overconfident. They are evaluated on the same patients they were trained on, which leads to performance estimates that are systematically too optimistic.

Second, models are often too rigid. Continuous clinical variables — age, creatinine, hemoglobin, lactate — are forced into linear relationships that do not reflect physiology.

Bootstrap validation and cubic splines in regression address these two issues directly.

Not by replacing regression, but by making it more honest and more realistic.

The first problem: optimism in clinical models

Imagine building a logistic regression model to predict 30-day mortality after cardiac surgery.

You use 18 preoperative variables, fit the model, and obtain an AUC of 0.85.

This number feels reassuring.

But where does it come from?

Almost always, it comes from evaluating the model on the same dataset used to estimate the coefficients. The model has already “seen” every patient. It has implicitly learned not only the signal, but also the noise.

In clinical terms, this is like testing a diagnostic score on the same cohort that was used to define it.

The result is not wrong — but it is optimistic.

The real question clinicians care about is different:

How will this model behave on the next patient?

Bootstrap validation: simulating future patients

Bootstrap validation answers this question without requiring an external cohort.

Suppose you have a dataset of 500 cardiac surgery patients.

A bootstrap sample is obtained by randomly drawing 500 patients with replacement from this dataset.

Some patients will appear multiple times.

Others will not appear at all.

On average:

  • about 63% of patients are included at least once (this derives from the probability (1-1/n)^n → e^(-1) ≈ 0.632 as n grows large),
  • about 37% are left out.

Those left-out patients are called out-of-bag (OOB) observations.

They play a crucial role.

Each bootstrap sample represents a plausible alternative reality in which the same study was conducted, but patient inclusion happened slightly differently.

Where optimism is measured

For each bootstrap sample, we do something very specific:

  • We fit the regression model on the bootstrap sample.
  • We evaluate its performance on:
    • the bootstrap sample itself (patients the model has “seen”),
    • the OOB patients (patients the model has never seen).

The difference between these two performances is optimism.

In a mortality model, for example, we may observe:

  • AUC = 0.87 on bootstrap data,
  • AUC = 0.81 on OOB data.

The optimism for that bootstrap iteration is 0.06.

Repeating this process hundreds or thousands of times produces a distribution of optimism estimates.

Their average represents how much our apparent performance is inflated.

Correcting clinical performance estimates

Once optimism is estimated, correction is straightforward.

If the apparent AUC of the original model is 0.85 and the average optimism is 0.05, the optimism-corrected AUC is:

0.85 − 0.05 = 0.80

This corrected value is a much better approximation of what we should expect in new patients.

This approach is now standard in high-quality clinical prediction modeling and is strongly recommended over split-sample validation when datasets are limited.

infopraphic: bootstrap and optimism

The second problem: linearity does not reflect physiology

Even a perfectly validated regression model can still be misleading if its structure is wrong.

Consider age as a predictor of mortality.

Is the risk increase from 40 to 50 the same as from 80 to 90?

Clinically, clearly not.

The same applies to creatinine, lactate, hemoglobin, or blood pressure.

These relationships are nonlinear, often with thresholds, plateaus, or acceleration zones.

Forcing them into a straight line is convenient, but biologically implausible.

Cubic splines: letting data bend smoothly

Cubic splines allow regression models to adapt to nonlinear patterns while remaining smooth and interpretable.

Instead of fitting one equation across the entire range of a variable, splines divide the range into intervals separated by knots.

Within each interval, a cubic polynomial is fitted, but all pieces are constrained to join smoothly.

This ensures:

  • continuity of the curve,
  • continuity of the first and second derivatives,
  • no abrupt changes in slope.

From a clinical perspective, this means:

  • no artificial jumps in risk,
  • no arbitrary cutoffs,
  • a continuous physiological interpretation.

A clinical example: creatinine and mortality

Creatinine often has a weak association with outcome at low values, but risk increases sharply beyond certain thresholds.

A linear term cannot capture this behavior.

A cubic spline can.

Instead of deciding a priori where risk changes, the spline allows the data to reveal the shape of the relationship — smoothly.

infographic: understanding cubic spline

Why splines and bootstrap belong together

Cubic splines increase model flexibility.

Flexibility increases the risk of overfitting.

Bootstrap validation controls that risk.

Together, they form a powerful and principled modeling strategy:

  • splines improve biological realism,
  • bootstrap ensures honest performance estimation.

This combination is common in high-impact clinical models, even if not always explicitly stated.

Common mistakes in medical applications

The most frequent error is using splines without validation.

Flexible models must be validated more carefully, not less.

Another common mistake is replacing continuous variables with categories.

This throws away information and introduces artificial thresholds.

Splines typically outperform categorization in terms of both accuracy and interpretability.

Limitations and When to Use Alternatives

Bootstrap validation and splines are powerful, but not universal solutions.

Computational cost: Bootstrap with 1000 iterations can be slow with large datasets or complex models. If computation is a bottleneck, consider k-fold cross-validation as a faster alternative.

Very small samples: With fewer than 100 observations, bootstrap may be unstable. Leave-one-out cross-validation or penalized regression (ridge, lasso) may be more appropriate.

Rare outcomes: When the outcome occurs in less than 5-10% of cases, stratified cross-validation ensures adequate representation in each fold. Standard bootstrap may occasionally produce samples with very few events.

Interpretation requirements: While spline curves are interpretable, if you need simple risk scores for bedside use (e.g., “add 3 points if age > 65”), linear models with categorization may be more practical despite statistical drawbacks.

The methods described here are most valuable when:

  • you have moderate sample sizes (100-5000 observations),
  • you need realistic performance estimates,
  • biological plausibility matters,
  • the model will be used for individual predictions rather than simple screening.

A fully reproducible example: bootstrap and cubic splines on simulated medical data

Step 1 — Create a synthetic clinical dataset

  • Age
  • Creatinine
  • Hemoglobin
  • Outcome with nonlinear risk (clinically plausible)
import numpy as np
import pandas as pd

np.random.seed(42)

n = 600

age = np.random.normal(70, 8, n).clip(40, 90)
creatinine = np.random.lognormal(mean=0.2, sigma=0.4, size=n).clip(0.5, 5)
hemoglobin = np.random.normal(13, 1.5, n).clip(8, 18)

# Nonlinear true risk function (unknown to the model)
logit = (
    0.04 * (age - 65)
    + 0.8 * np.maximum(creatinine - 1.2, 0) ** 1.5
    - 0.25 * (hemoglobin - 13)
)

prob = 1 / (1 + np.exp(-logit))
mortality = np.random.binomial(1, prob)

data = pd.DataFrame({
    "age": age,
    "creatinine": creatinine,
    "hemoglobin": hemoglobin,
    "death_30d": mortality
})

data.head()

👉 This dataset is nonlinear by design, just like clinical reality.

Step 2 — Apparent performance of a simple logistic regression

Let’s start with a classic linear model.

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score

X = data[["age", "creatinine", "hemoglobin"]]
y = data["death_30d"]

model = LogisticRegression(max_iter=1000)
model.fit(X, y)

apparent_auc = roc_auc_score(y, model.predict_proba(X)[:, 1])
print(f"Apparent AUC: {apparent_auc:.3f}")

Result: Apparent AUC: 0.661🔍

This AUC is optimistic: the model is evaluated on the same patients.

Step 3 — Bootstrap optimism correction (fully reproducible)

Now let’s estimate the optimism.

from sklearn.utils import resample

n_boot = 500
optimism = []

for i in range(n_boot):
    boot_idx = resample(np.arange(len(data)), replace=True)
    oob_idx = np.setdiff1d(np.arange(len(data)), boot_idx)

    if len(oob_idx) < 30:
        continue

    X_boot, y_boot = X.iloc[boot_idx], y.iloc[boot_idx]
    X_oob, y_oob = X.iloc[oob_idx], y.iloc[oob_idx]

    model.fit(X_boot, y_boot)

    auc_boot = roc_auc_score(y_boot, model.predict_proba(X_boot)[:, 1])
    auc_oob = roc_auc_score(y_oob, model.predict_proba(X_oob)[:, 1])

    optimism.append(auc_boot - auc_oob)

mean_optimism = np.mean(optimism)
corrected_auc = apparent_auc - mean_optimism

print(f"Mean optimism: {mean_optimism:.3f}")
print(f"Optimism-corrected AUC: {corrected_auc:.3f}")

Result: Mean optimism: 0.014 Optimism-corrected AUC: 0.647

Step 4 — Why linear regression is inadequate here

Now let’s look at the true (simulated) relationship between creatinine and risk.

import matplotlib.pyplot as plt

plt.scatter(data["creatinine"], prob, alpha=0.3)
plt.xlabel("Creatinine (mg/dL)")
plt.ylabel("True mortality risk")
plt.title("True nonlinear relationship (unknown to the model)")
plt.show()

👉 No clinician would expect a linear relationship.

Yet the previous model forces it.

Step 5 — Logistic regression with cubic splines

Now we introduce cubic splines.

import statsmodels.api as sm
from patsy import dmatrix

spline_creatinine = dmatrix(
    "bs(creatinine, df=4, include_intercept=False)",
    data,
    return_type="dataframe"
)

X_spline = pd.concat([
    data[["age", "hemoglobin"]],
    spline_creatinine
], axis=1)

X_spline = sm.add_constant(X_spline)

model_spline = sm.Logit(y, X_spline).fit(disp=False)

pred_spline = model_spline.predict(X_spline)
spline_auc = roc_auc_score(y, pred_spline)

print(f"Spline model apparent AUC: {spline_auc:.3f}")

Result: Spline model apparent AUC: 0.669

The model now adapts to physiology, not the other way around.

Step 6 — Bootstrap validation of the spline model

Flexibility must be paid for with rigorous validation.

optimism_spline = []

for i in range(n_boot):
    boot_idx = resample(np.arange(len(data)), replace=True)
    oob_idx = np.setdiff1d(np.arange(len(data)), boot_idx)

    if len(oob_idx) < 30:
        continue

    data_boot = data.iloc[boot_idx]
    data_oob = data.iloc[oob_idx]

    Xb = sm.add_constant(pd.concat([
        data_boot[["age", "hemoglobin"]],
        dmatrix("bs(creatinine, df=4, include_intercept=False)",
                data_boot, return_type="dataframe")
    ], axis=1))

    Xo = sm.add_constant(pd.concat([
        data_oob[["age", "hemoglobin"]],
        dmatrix("bs(creatinine, df=4, include_intercept=False)",
                data_oob, return_type="dataframe")
    ], axis=1))

    yb = data_boot["death_30d"]
    yo = data_oob["death_30d"]

    try:
        m = sm.Logit(yb, Xb).fit(disp=False)
        auc_b = roc_auc_score(yb, m.predict(Xb))
        auc_o = roc_auc_score(yo, m.predict(Xo))
        optimism_spline.append(auc_b - auc_o)
    except:
        continue

mean_opt_spline = np.mean(optimism_spline)
corrected_auc_spline = spline_auc - mean_opt_spline

print(f"Spline optimism: {mean_opt_spline:.3f}")
print(f"Spline corrected AUC: {corrected_auc_spline:.3f}")

Result: Spline optimism: 0.030 Spline corrected AUC: 0.639

Conclusion — A modern view of regression in medicine

Regression is not outdated.

It is underutilized.

Bootstrap validation teaches regression humility.

Cubic splines give it flexibility.

Together, they transform regression from a blunt instrument into a refined clinical modeling tool.

If regression is the engine of medical prediction,

bootstrap and splines are what finally put the turbo on it.

References and Further Reading

Harrell FE. Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis. 2nd ed. Springer, 2015.

The definitive reference on modern regression techniques in biostatistics, including extensive coverage of splines and bootstrap validation.

Steyerberg EW. Clinical Prediction Models: A Practical Approach to Development, Validation, and Updating. 2nd ed. Springer, 2019.

Comprehensive guide to building and validating prediction models in medicine, with emphasis on optimism correction.

Efron B, Tibshirani RJ. An Introduction to the Bootstrap. Chapman & Hall, 1993.

The foundational text on bootstrap methods by their inventor.

Collins GS, Reitsma JB, Altman DG, Moons KGM. Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD): the TRIPOD Statement. BMJ 2015;350:g7594.

Essential guidelines for reporting clinical prediction models, emphasizing proper validation.

A young boy kneels beside the shattered remains of a clay pot in a sunlit, rustic room, surrounded by warm golden light, long shadows, and aged earthenware in a melancholic early-20th-century painterly style.

Multiple Imputation for Missing Data

Posted on by Michele Danilo Pierri

A Complete Guide to Multiple Imputation for Missing Data: When, How, and Why

Last updated: November 2025

Author: Michele D. Pierri

Reading time: 15–20 minutes

Glossary

MCAR: Missing Completely At Random. Missingness is unrelated to observed or unobserved data; complete case can be unbiased.

MAR: Missing At Random. Missingness depends only on observed variables; standard MI assumptions target MAR.

MNAR: Missing Not At Random. Missingness depends on unobserved values; requires sensitivity analyses or explicit models.

MI (Multiple Imputation): Generate m completed datasets, analyze each, pool results.

MICE: Multiple Imputation by Chained Equations; iterative conditional models.

Rubin’s rules: Pooling framework using Q̄ (mean estimate), Ū (within‑imputation variance), B (between‑imputation variance), T (total variance), df (Barnard–Rubin), λ/FMI (fraction of missing information).

SE: Standard Error. CI: Confidence Interval. OLS: Ordinary Least Squares.

OHE: One‑hot encoding for categorical variables.

IterativeImputer: scikit‑learn’s MICE‑style imputer (experimental) for numeric arrays.

FMI (λ): Fraction of Missing Information; guides how large m should be.

Table of Contents

Why Multiple Imputation Matters for Missing Data

When to Use Multiple Imputation

Understanding the Missing Data Mechanisms

How Multiple Imputation Works for Missing Data

Python Implementation with Synthetic Datasets

Choosing Parameters and Estimators

Pooling Results: The Final Step (Rubin’s Rules)

Inference vs Prediction: Two Playbooks

Diagnostics, Constraints, and Sensitivity (MNAR)

Comparison with Other Methods

Complete End to End Example with Visualization

Reproducibility

External Resources

FAQ

1. Why Multiple Imputation Matters for Missing Data

Multiple Imputation (MI) is a rigorous approach to handling missing data that preserves uncertainty and reduces bias compared with complete case analysis or single imputation. MI generates multiple plausible values for each missing entry, enables valid inference, and preserves statistical power.

Key points:

Valid inference: includes uncertainty due to missingness in SEs and CIs.

Efficiency: typically more precise than complete case.

Flexibility: compatible with many analyses (regression, t‑tests, ML).

Reduced bias: especially effective when data are MAR.

2. When to Use Multiple Imputation

Use MI when

Missingness > 5–10% or non‑trivial loss of power.

Mechanism is MCAR or MAR (MAR is the typical scenario for MI).

You need valid SEs/CIs for inference.

Multivariable analyses with complex relations and/or informative auxiliary variables.

You may avoid MI when

MCAR with <5% missing: complete case can be adequate and unbiased.

Extremely small datasets (n < 50) without auxiliary information.

Purely predictive objective and you do not need uncertainty quantification: single or deterministic imputation can be sufficient.

Caution

MNAR: standard MI assumes MAR. With MNAR you need sensitivity analyses or missingness models.

3. Understanding the Missing Data Mechanisms

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(42)
n = 1000

data = {
    'age': np.random.normal(45, 15, n),
    'income': np.random.lognormal(10.5, 0.8, n),
    'education': np.random.choice(['High School', 'Bachelor', 'Master', 'PhD'], n),
}
df_complete = pd.DataFrame(data)

# MCAR: random missingness
df_mcar = df_complete.copy()
mcar_mask = np.random.choice([True, False], n, p=[0.2, 0.8])
df_mcar.loc[mcar_mask, 'income'] = np.nan

# MAR: missingness related to age
df_mar = df_complete.copy()
mar_prob = 1 / (1 + np.exp(-(df_mar['age'] - 60) / 10))
mar_mask = np.random.random(n) < mar_prob
df_mar.loc[mar_mask, 'income'] = np.nan

print('MCAR rate:', df_mcar['income'].isnull().mean())
print('MAR rate:', df_mar['income'].isnull().mean())

sns.set(style="whitegrid")

def plot_missing_by_age(df_mcar, df_mar, n_bins=12):
    bins = np.linspace(df_complete['age'].min(), df_complete['age'].max()
    def rates(df):
        idx = np.digitize(df['age'], bins) - 1
        mids = [(bins[i] + bins[i+1]) / 2 for i in range(len(bins)-1)]
        rates = []
        for i in range(len(bins)-1):
            sel = idx == i
            rates.append(df.loc[sel, 'income'].isna().mean() if sel.sum()
        return np.array(mids), np.array(rates)
    x_m, y_mcar = rates(df_mcar)
    x_mar, y_mar = rates(df_mar)

    plt.figure(figsize=(8, 4))
    plt.plot(x_m, y_mcar, marker='o', label='MCAR')
    plt.plot(x_mar, y_mar, marker='o', label='MAR')
    plt.xlabel('Age')
    plt.ylabel('Missing rate (income)')
    plt.title('Missing rate by age (binned): MCAR vs MAR')
    plt.legend()
    plt.tight_layout()
    plt.show()


def plot_missing_by_education(df_mcar, df_mar):
    order = sorted(df_complete['education'].unique(), key=lambda x: ['Hig
    rm = df_mcar.groupby('education')['income'].apply(lambda s: s.isna().
    rmar = df_mar.groupby('education')['income'].apply(lambda s: s.isna()
    df_plot = pd.DataFrame({'MCAR': rm, 'MAR': rmar})
    ax = df_plot.plot.bar(rot=0, figsize=(8,4))
    ax.set_ylabel('Missing rate (income)')
    ax.set_title('Missing rate by education: MCAR vs MAR')
    plt.tight_layout()
    plt.show()


def scatter_income_age_with_missing(df, title):
    plt.figure(figsize=(8,4))
    observed = df.loc[df['income'].notna()]
    missing = df.loc[df['income'].isna()]
    plt.scatter(observed['age'], observed['income'], s=12, alpha=0.6, lab
    # Plot observations with missing income as points with y-jitter (for 
    if len(missing) > 0:
        y_jitter = np.random.uniform(observed['income'].min(), observed['
        plt.scatter(missing['age'], y_jitter, s=12, alpha=0.6, label='mis
    plt.yscale('log')
    plt.xlabel('Age')
    plt.ylabel('Income (log scale)')
    plt.title(title + ' — missing highlighted')
    plt.legend()
    plt.tight_layout()
    plt.show()



plot_missing_by_age(df_mcar, df_mar, n_bins=12)
plot_missing_by_education(df_mcar, df_mar)
scatter_income_age_with_missing(df_mcar, 'MCAR dataset')
scatter_income_age_with_missing(df_mar, 'MAR dataset')
Line graph of missing rate by age showing MCAR and MAR missing data

4. How Multiple Imputation Works

Three classic stages:

1) Imputation: generate m completed datasets by imputing from predictive distributions.

2) Analysis: run your analysis separately on each of the m datasets.

3) Pooling: combine estimates and variances using Rubin’s rules.

What it does: defines a reusable helper to create m imputed datasets with sensible defaults and bounds.

Why it matters: keeps randomness across imputations and enforces basic plausibility constraints.

# Helper to generate m imputed datasets with MI-friendly settings
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.impute import IterativeImputer
from sklearn.linear_model import BayesianRidge
import numpy as np, pandas as pd

def generate_imputed_datasets(df, m=20, random_state=42, **kwargs):
    """Return a list of m DataFrames imputed via IterativeImputer.
    Parameters
    ----------
    df : pandas.DataFrame (numeric or already encoded)
    m : int, number of imputations
    random_state : int, base seed; per-imputation seed = base + i
    kwargs : optional overrides (estimator, max_iter, bounds, etc.)
    """
    sets = []
    for i in range(m):
        imp = IterativeImputer(
            estimator=kwargs.get('estimator', BayesianRidge()),
            sample_posterior=True,   # crucial for MI (injects posterior noise)
            max_iter=kwargs.get('max_iter', 10),
            random_state=random_state + i,
            n_nearest_features=kwargs.get('n_nearest_features', None),
            min_value=kwargs.get('min_value', None),  # set for domain bounds
            max_value=kwargs.get('max_value', None)
        )
        # Ensure we only pass numeric columns. If categoricals exist, encode before calling this helper.
        numeric_df = df.select_dtypes(include=[np.number])
        imputed = imp.fit_transform(numeric_df)
        imputed_df = pd.DataFrame(imputed, columns=numeric_df.columns, index=df.index)
        # Preserve any untouched non-numeric columns by concatenation
        non_numeric = df.drop(columns=list(numeric_df.columns), errors='ignore')
        out = pd.concat([imputed_df, non_numeric], axis=1)[df.columns]
        sets.append(out)
    return sets

> Important: IterativeImputer provides MICE‑style imputations but not a full MI ecosystem for inference. For pooling and statistical diagnostics use dedicated functions or libraries like statsmodels or mice (R).

5. Python Implementation with Synthetic Datasets

Example 1: Regression Analysis with Multiple Imputation

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

np.random.seed(123)
n = 500
X1 = np.random.normal(10, 2, n)
X2 = 0.5 * X1 + np.random.normal(0, 1, n)
Y  = 3 + 2*X1 + 1.5*X2 + np.random.normal(0, 3, n)

df = pd.DataFrame({'X1': X1, 'X2': X2, 'Y': Y})
# MAR
missing_X2 = (df['X1'] > 12) & (np.random.random(n) < 0.3)
missing_Y  = (df['X2'] > df['X2'].mean()) & (np.random.random(n) < 0.25)
df.loc[missing_X2, 'X2'] = np.nan
df.loc[missing_Y,  'Y']  = np.nan

m = 20
imputed_sets = generate_imputed_datasets(df, m=m, random_state=42)

results = []
for di, dataset in enumerate(imputed_sets):
    X = sm.add_constant(dataset[['X1','X2']])
    y = dataset['Y']
    model = sm.OLS(y, X).fit()
    # collect estimates and SEs
    results.append({
        'intercept': model.params['const'], 'se_intercept': model.bse['const'],
        'coef_X1': model.params['X1'],     'se_X1': model.bse['X1'],
        'coef_X2': model.params['X2'],     'se_X2': model.bse['X2'],
        'r_squared': model.rsquared
    })

results_df = pd.DataFrame(results)

Example 2: Mixed Types without artificial ordering

import pandas as pd, numpy as np
from sklearn.experimental import enable_iterative_imputer  # noqa: F401
from sklearn.impute import IterativeImputer
from sklearn.preprocessing import OneHotEncoder

np.random.seed(456)
n = 300
df_mixed = pd.DataFrame({
    'age': np.random.normal(40, 12, n),
    'income': np.random.lognormal(10, 0.7, n),
    'education': np.random.choice(['HS','BA','MA','PhD'], n),
    'city': np.random.choice(['NYC','LA','Chicago','Boston'], n),
    'satisfaction': np.random.normal(7, 2, n)
})

df_mixed.loc[df_mixed.sample(50).index, 'income'] = np.nan
df_mixed.loc[df_mixed.sample(40).index, 'education'] = np.nan

# One‑hot for categorical variables during imputation
cat_cols = ['education','city']
num_cols = [c for c in df_mixed.columns if c not in cat_cols]

ohe = OneHotEncoder(sparse_output=False, handle_unknown='ignore')
X_cat = pd.DataFrame(ohe.fit_transform(df_mixed[cat_cols]), columns=ohe.get_feature_names_out(cat_cols))
X = pd.concat([df_mixed[num_cols].reset_index(drop=True), X_cat.reset_index(drop=True)], axis=1)

imp = IterativeImputer(sample_posterior=True, max_iter=10, random_state=789)
imputed = pd.DataFrame(imp.fit_transform(X), columns=X.columns)

# Bring categories back with inverse_transform
cat_imputed = ohe.inverse_transform(imputed[ohe.get_feature_names_out(cat_cols)])
for j, col in enumerate(cat_cols):
    df_mixed[col] = cat_imputed[:, j]

6. Choosing Parameters and Estimators

Key parameters (IterativeImputer)

  • sample_posterior=True: required for MI.
  • m (number of imputations): 20–100. Practical guide: increase m with missingness and with higher fraction of missing information (FMI).
  • max_iter: 10–20, increase if not converged.
  • estimator: BayesianRidge for continuous; RandomForest/ExtraTrees for nonlinear patterns and mixed types.
  • n_nearest_features: None or ~√p to speed up with many features.

Estimator choice: practical note

  • Bayesian linear models: stable, fast, interpretable.
  • RandomForest/ExtraTrees: robust for mixed data, watch for overfitting and runtime.

From completed imputations to pooling: what really happens between stages

  • After imputation you have m complete datasets, all analyzed with the exact same model. From each analysis you only need two things per parameter of interest: the point estimate (Q, e.g., a regression coefficient) and its estimated variance (U, i.e., SE² from the fit). You don’t need the imputed datasets during pooling, just the numerical summaries Q and U from each of the m analyses.
  • The key hand‑off is to separate the sources of randomness. Randomness lives in the imputation step (sample_posterior=True with different seeds), while the analysis on each completed dataset should be deterministic and identical across datasets (same formula, transformations, and estimator). This way, between‑imputation variability reflects only uncertainty due to missingness, which Rubin then combines with the “ordinary” within‑analysis uncertainty.

How to collect what you need for pooling (practical checklist)

  1. Fit the same model on each of the m datasets.
  2. For every parameter p, store: Q_p^(j) and U_p^(j)=SE_p^(j)² for j=1..m. Optional but useful: global metrics (R², AIC) for reporting.
  3. Apply Rubin’s rules: compute Q̄_p (mean of the Qs), Ū_p (mean of the Us), B_p (variance of the Qs), then T_p = Ū_p + (1 + 1/m)B_p. From T get SE; use Barnard–Rubin for degrees of freedom, then CI and p‑value. Report FMI (λ) to show information loss from missingness and to justify m.

7. Pooling Results: The Final Step

Rubin’s rules combine m estimates:

  1. Q̄: mean of point estimates.
  2. Ū: mean within‑imputation variance.
  3. B: between‑imputation variance.
  4. T = Ū + (1 + 1/m)B.

Also include:

  • λ = (B + B/m) / T (fraction of missing information, FMI).
  • df with Barnard–Rubin correction.
Distribution of Observed vs Imputed Income Values
from scipy import stats
import numpy as np, pandas as pd

def pool_rubins(results_df, params=('intercept','coef_X1','coef_X2'), conf=0.95):
    m = len(results_df)
    out = {}
    for p in params:
        if p == 'intercept': se_key = 'se_intercept'
        elif p == 'coef_X1': se_key = 'se_X1'
        elif p == 'coef_X2': se_key = 'se_X2'
        else: se_key = f'se_{p}'
        Q_bar = results_df[p].mean()
        U_bar = (results_df[se_key]**2).mean()
        B = results_df[p].var(ddof=1)
        T = U_bar + (1 + 1/m)*B
        FMI = (B + B/m) / T if T > 0 else 0.0
        # Barnard–Rubin df
        if B > 0 and U_bar > 0:
            df = (m - 1) * (1 + U_bar / ((1 + 1/m) * B))**2
        else:
            df = np.inf
        se = np.sqrt(T)
        tcrit = stats.t.ppf(0.5 + conf/2, df) if np.isfinite(df) else stats.norm.ppf(0.5 + conf/2)
        ci_lo, ci_hi = Q_bar - tcrit*se, Q_bar + tcrit*se
        out[p] = {'estimate': Q_bar, 'se': se, 'df': df, 'ci_lower': ci_lo, 'ci_upper': ci_hi, 'FMI': FMI}
    return pd.DataFrame(out).T

pooled = pool_rubins(results_df)
print(pooled.round(4))

8. Inference vs Prediction: Two Playbooks

For inference

  • Objective: estimates, SEs, CIs, p‑values.
  • Pipeline: m imputations → fit model on each → pool with Rubin (incl. df and FMI).
  • Recommended libraries: statsmodels; in R: mice.

For prediction

  • Objective: minimize predictive error.
  • Strategies:
    • Train m models and average per‑case predictions.
    • Or evaluate the metric on each imputed dataset then average/pool metrics (report mean ± SD or bootstrap CI).
  • Deterministic single imputation is often sufficient for purely predictive goals if you don’t need uncertainty.

9. Diagnostics, Constraints, and Sensitivity (MNAR)

Convergence diagnostics (Python)

import matplotlib.pyplot as plt

def trace_means_over_iter(imputer, X, var_idx=0):
    # Illustrative placeholder: capture per‑iteration stats using imputer.verbose
    pass  # Implementation depends on version; R mice provides native traces

Minimum diagnostic checklist:

  • Convergence: n_iter_ < max_iter and stability of imputed means/variances.
  • Distributions: compare observed vs imputed distributions.
  • Clinical plausibility: coherence checks across related variables.

Constraints and clinical plausibility

  • Use min_value and max_value in IterativeImputer to respect physical ranges.
  • Monotone transforms for strictly positive variables (log‑transform).
  • Post‑processing with domain rules (e.g., age ≥ 0, creatinine ↔ eGFR coherence).

MNAR: sensitivity analysis (pattern‑mixture with delta‑adjustment)

delta = 0.2  # e.g., increase imputed income by 20%
adj_sets = []
for ds in imputed_sets:
    ds_adj = ds.copy()
    ds_adj['income'] = ds_adj['income'] * (1 + delta)
    adj_sets.append(ds_adj)
# Re‑run analysis and compare pooled results (stability → robustness)

10. Comparison with Other Methods

For inference

Complete case: unbiased only under MCAR; reduced power.

Single imputation: underestimated SEs; not recommended for inference.

Multiple imputation: less biased estimates, correct SEs/CIs via pooling.

For prediction

Mean/Median/Mode imputation: simple baselines.

KNN/RandomForest imputation: strong in some scenarios, but bias is scenario‑dependent and uncertainty is unquantified.

Avoid absolute claims like “low bias” for KNN: it depends on data‑generating process and missingness pattern.

11. Complete End to End Example with Visualization

Below is a comprehensive workflow that demonstrates the entire multiple imputation process with publication-ready visualizations.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
from sklearn.linear_model import BayesianRidge, LinearRegression
from sklearn.ensemble import RandomForestRegressor
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

# Set style for publication-quality plots
sns.set_style("whitegrid")
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 11

# ============================================
# STEP 1: Generate Synthetic Dataset with MAR
# ============================================
print("="*60)
print("STEP 1: Generating Synthetic Dataset with MAR Missingness")
print("="*60)

np.random.seed(42)
n = 1000

# Create correlated features
education_years = np.random.normal(12, 3, n)
income = 20000 + 3000 * education_years + np.random.normal(0, 8000, n)
health_score = 70 + 0.5 * education_years + 0.0001 * income + np.random.normal(0, 10, n)
job_satisfaction = np.random.normal(7, 1.5, n)

df_full = pd.DataFrame({
    'education_years': education_years,
    'income': income,
    'health_score': health_score,
    'job_satisfaction': job_satisfaction
})

# Introduce MAR missingness: income missing depends on education & health
missing_prob = 1 / (1 + np.exp(-(-2 + 0.2 * education_years - 0.02 * health_score)))
missing_mask = np.random.random(n) < missing_prob
df_full.loc[missing_mask, 'income'] = np.nan

# Also make some health scores MCAR for comparison
df_full.loc[np.random.choice(n, 50), 'health_score'] = np.nan

print(f"Dataset shape: {df_full.shape}")
print(f"Missing values:\n{df_full.isnull().sum()}")
print(f"Overall missingness rate: {df_full.isnull().sum().sum() / (df_full.shape[0] * df_full.shape[1]):.2
print(f"Income missingness rate: {df_full['income'].isnull().mean():.2%}")

# ============================================
# STEP 2: Visualize Missingness Pattern
# ============================================

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Original data distribution
axes[0, 0].scatter(df_full['education_years'], df_full['income'], alpha=0.6, s=20)
axes[0, 0].set_title('Original Data: Education vs Income')
axes[0, 0].set_xlabel('Education Years')
axes[0, 0].set_ylabel('Income')

# Missingness indicator heatmap
missing_indicator = df_full.isnull().astype(int)
sns.heatmap(missing_indicator.head(50), cmap='Reds', cbar=True, ax=axes[0, 1])
axes[0, 1].set_title('Missingness Pattern (Red = Missing)')
axes[0, 1].set_xlabel('Variables')
axes[0, 1].set_ylabel('First 50 Observations')

# MAR verification: income missingness vs education
observed = df_full['income'].notna()
axes[1, 0].hist(df_full.loc[observed, 'education_years'], bins=30, alpha=0.7, label='Observed Income', den
axes[1, 0].hist(df_full.loc[~observed, 'education_years'], bins=30, alpha=0.7, label='Missing Income', den
axes[1, 0].set_title('MAR Check: Education Distribution by Missingness')
axes[1, 0].set_xlabel('Education Years')
axes[1, 0].legend()

# Correlation matrix
corr_data = df_full.corr()
sns.heatmap(corr_data, annot=True, cmap='coolwarm', center=0, ax=axes[1, 1])
axes[1, 1].set_title('Correlation Matrix (with missing values)')

plt.tight_layout()
plt.savefig('missingness_analysis.png', dpi=300, bbox_inches='tight')
plt.show()

# ============================================
# STEP 3: Generate Multiple Imputations (m=30)
# ============================================

print("\n" + "="*60)
print("STEP 2: Generating m=30 Imputed Datasets")
print("="*60)

m = 30
imputed_datasets = []

for i in range(m):
    print(f"Generating imputation {i+1}/{m}...", end='\r')
    imputer = IterativeImputer(
        estimator=BayesianRidge(),
        sample_posterior=True,
        max_iter=20,
        random_state=42 + i,
        add_indicator=False,
        verbose=0
    )
    
    imputed_data = imputer.fit_transform(df_full)
    imputed_datasets.append(pd.DataFrame(imputed_data, columns=df_full.columns))

print(f"\nGenerated {len(imputed_datasets)} imputed datasets")
print(f"Each dataset shape: {imputed_datasets[0].shape}")

# ============================================
# STEP 4: Analyze Each Imputed Dataset
# ============================================

print("\n" + "="*60)
print("STEP 3: Analyzing Each Dataset (Linear Regression)")
print("="*60)

regression_results = []

for i, dataset in enumerate(imputed_datasets):
    # Regression: income ~ education_years + health_score
    X = dataset[['education_years', 'health_score']]
    y = dataset['income']
    
    # Add constant for intercept
    X_with_const = np.column_stack([np.ones(len(X)), X])
    
    # Fit regression
    coefficients = np.linalg.lstsq(X_with_const, y, rcond=None)[0]
    
    # Calculate predictions and residuals
    y_pred = X_with_const @ coefficients
    residuals = y - y_pred
    
    # Standard errors
    mse = np.sum(residuals**2) / (len(y) - X.shape[1] - 1)
    var_coef = mse * np.linalg.inv(X_with_const.T @ X_with_const).diagonal()
    se_coef = np.sqrt(var_coef)
    
    regression_results.append({
        'imputation': i,
        'intercept': coefficients[0],
        'se_intercept': se_coef[0],
        'coef_education': coefficients[1],
        'se_education': se_coef[1],
        'coef_health': coefficients[2],
        'se_health': se_coef[2],
        'r_squared': 1 - np.sum(residuals**2) / np.sum((y - np.mean(y))**2)
    })

results_df = pd.DataFrame(regression_results)
print("First 5 regression results:")
print(results_df[['imputation', 'coef_education', 'coef_health', 'r_squared']].head())

# ============================================
# STEP 5: Pool Results using Rubin's Rules
# ============================================

print("\n" + "="*60)
print("STEP 4: Pooling Results (Rubin's Rules)")
print("="*60)

def rubins_rules_complete(results_df, conf_level=0.95):
    """Complete implementation of Rubin's rules with Fraction of Missing Information"""
    m = len(results_df)
    pooled = {}
    
    params = ['intercept', 'coef_education', 'coef_health']
    
    # mapping from pooled param name -> se column in results_df
    se_col_map = {
        'intercept': 'se_intercept',
        'coef_education': 'se_education',
        'coef_health': 'se_health'
    }
    
    for param in params:
        # Point estimate (Q_bar)
        Q_bar = results_df[param].mean()
        
        # Within-imputation variance (U_bar)
        se_col = se_col_map[param]
        U_bar = (results_df[se_col]**2).mean()
        
        # Between-imputation variance (B)
        B = results_df[param].var(ddof=1)
        
        # Total variance (T)
        T = U_bar + (1 + 1/m) * B
        
        # Standard error
        se_Q_bar = np.sqrt(T)
        
        # Degrees of freedom (df)
        if B > 0:
            df = (m - 1) * (1 + U_bar / ((1 + 1/m) * B))**2
        else:
            df = np.inf
        
        # Confidence interval
        t_critical = stats.t.ppf((1 + conf_level) / 2, df)
        ci_lower = Q_bar - t_critical * se_Q_bar
        ci_upper = Q_bar + t_critical * se_Q_bar
        
        # Fraction of Missing Information (FMI)
        fmi = (1 + 1/m) * B / T
        
        pooled[param] = {
            'estimate': Q_bar,
            'se': se_Q_bar,
            'df': df,
            'ci_lower': ci_lower,
            'ci_upper': ci_upper,
            'fmi': fmi,
            'between_var': B,
            'within_var': U_bar
        }
    
    return pd.DataFrame(pooled).T

pooled_results = rubins_rules_complete(results_df)
print("Pooled Regression Results:")
print(pooled_results.round(4))

# ============================================
# STEP 6: Visualize Results
# ============================================
# use a distinct name for the main figure
fig_main, axes = plt.subplots(2, 3, figsize=(18, 12))

# Plot 1: Distribution of coefficients across imputations
axes[0, 0].hist(results_df['coef_education'], bins=15, color='steelblue', alpha=0.7, edgecolor='black')
axes[0, 0].axvline(pooled_results.loc['coef_education', 'estimate'], color='red', linewidth=2, label='Pool
axes[0, 0].axvline(pooled_results.loc['coef_education', 'ci_lower'], color='red', linestyle='--', linewidt
axes[0, 0].axvline(pooled_results.loc['coef_education', 'ci_upper'], color='red', linestyle='--', linewidt
axes[0, 0].set_title('Distribution of Education Coefficient\nAcross 30 Imputations')
axes[0, 0].set_xlabel('Coefficient Value')
axes[0, 0].legend()

# Plot 2: Health coefficient distribution
axes[0, 1].hist(results_df['coef_health'], bins=15, color='darkseagreen', alpha=0.7, edgecolor='black')
axes[0, 1].axvline(pooled_results.loc['coef_health', 'estimate'], color='red', linewidth=2, label='Pooled 
axes[0, 1].axvline(pooled_results.loc['coef_health', 'ci_lower'], color='red', linestyle='--', linewidth=1
axes[0, 1].axvline(pooled_results.loc['coef_health', 'ci_upper'], color='red', linestyle='--', linewidth=1
axes[0, 1].set_title('Distribution of Health Coefficient\nAcross 30 Imputations')
axes[0, 1].set_xlabel('Coefficient Value')
axes[0, 1].legend()

# Plot 3: R-squared across imputations
axes[0, 2].plot(results_df['imputation'], results_df['r_squared'], 'o-', color='purple', alpha=0.7)
axes[0, 2].set_title('R-squared Across Imputations')
axes[0, 2].set_xlabel('Imputation Number')
axes[0, 2].set_ylabel('R-squared')
axes[0, 2].grid(True, alpha=0.3)

# Plot 4: Compare imputed vs observed income distributions
# create a separate figure (do not overwrite fig_main)
fig4, ax = plt.subplots(figsize=(14, 5))
observed_income = df_full['income'].dropna()
ax.hist(observed_income, bins=40, alpha=0.6, label='Observed Income', density=True, color='blue')

# Plot imputed distributions (sample 5 imputations)
for i in range(min(5, len(imputed_datasets))):
    imputed_values = imputed_datasets[i].loc[df_full['income'].isna(), 'income']
    ax.hist(imputed_values, bins=40, alpha=0.3, density=True, label=f'Imputation {i+1}')

ax.set_title('Distribution of Observed vs Imputed Income Values')
ax.set_xlabel('Income')
ax.set_ylabel('Density')
ax.legend()
fig4.tight_layout()
fig4.savefig('imputed_distributions.png', dpi=300, bbox_inches='tight', transparent=False)
plt.show()

# Plot 5: Credible intervals for coefficients
coeff_names = ['Intercept', 'Education', 'Health']
estimates = [pooled_results.loc[param, 'estimate'] for param in ['intercept', 'coef_education', 'coef_heal
lower_bounds = [pooled_results.loc[param, 'ci_lower'] for param in ['intercept', 'coef_education', 'coef_h
upper_bounds = [pooled_results.loc[param, 'ci_upper'] for param in ['intercept', 'coef_education', 'coef_h

y_pos = np.arange(len(coeff_names))
axes[1, 0].barh(y_pos, [ub - lb for ub, lb in zip(upper_bounds, lower_bounds)],
                left=lower_bounds, height=0.6, alpha=0.6, color='lightcoral')
axes[1, 0].plot(estimates, y_pos, 'ko', markersize=8, label='Point Estimate')
axes[1, 0].set_yticks(y_pos)
axes[1, 0].set_yticklabels(coeff_names)
axes[1, 0].set_xlabel('Coefficient Value')
axes[1, 0].set_title('95% Confidence Intervals (Pooled Results)')
axes[1, 0].grid(True, alpha=0.3)

# Plot 6: Fraction of Missing Information
fmi_values = [pooled_results.loc[param, 'fmi'] for param in ['intercept', 'coef_education', 'coef_health']
colors = ['red' if fmi > 0.5 else 'orange' if fmi > 0.3 else 'green' for fmi in fmi_values]
bars = axes[1, 1].bar(coeff_names, fmi_values, color=colors, alpha=0.7, edgecolor='black')
axes[1, 1].set_ylabel('Fraction of Missing Information (FMI)')
axes[1, 1].set_title('Fraction of Missing Information\nby Parameter')
axes[1, 1].set_ylim(0, 1)
axes[1, 1].axhline(y=0.3, color='orange', linestyle='--', alpha=0.5, label='Moderate FMI')
axes[1, 1].axhline(y=0.5, color='red', linestyle='--', alpha=0.5, label='High FMI')
axes[1, 1].legend()

# Add value labels on bars
for bar, fmi in zip(bars, fmi_values):
    height = bar.get_height()
    axes[1, 1].text(bar.get_x() + bar.get_width()/2., height,
                    f'{fmi:.3f}', ha='center', va='bottom', fontweight='bold')

# Plot 7: Convergence check (plot imputed values across iterations)
axes[1, 2].set_title('Convergence Check\n(Not Available in scikit-learn)')
axes[1, 2].text(0.5, 0.5, 'Use mice package in R\nfor convergence diagnostics',
                ha='center', va='center', transform=axes[1, 2].transAxes,
                fontsize=12, style='italic')
axes[1, 2].set_xticks([])
axes[1, 2].set_yticks([])

# final save: save the main 2x3 figure using fig_main
fig_main.tight_layout()
fig_main.savefig('pooled_results_analysis.png', dpi=300, bbox_inches='tight', transparent=False)
plt.show()

# ============================================
# STEP 7: Compare with Alternative Methods
# ============================================

print("\n" + "="*60)
print("STEP 5: Comparison with Alternative Methods")
print("="*60)

def analyze_method(df_method, method_name):
    """Analyze a dataset using single method"""
    X = df_method[['education_years', 'health_score']]
    y = df_method['income']
    X_const = np.column_stack([np.ones(len(X)), X])
    coef = np.linalg.lstsq(X_const, y, rcond=None)[0]
    return coef[[1, 2]]  # Return only coefficients

# Complete case analysis
df_cc = df_full.dropna()
cc_coefs = analyze_method(df_cc, "Complete Case")

# Mean imputation
df_mean = df_full.fillna(df_full.mean())
mean_coefs = analyze_method(df_mean, "Mean Imputation")

# Single MICE (deterministic)
imputer_single = IterativeImputer(estimator=BayesianRidge(), random_state=42)
df_single_mice = pd.DataFrame(imputer_single.fit_transform(df_full), columns=df_full.columns)
single_mice_coefs = analyze_method(df_single_mice, "Single MICE")

# Multiple Imputation (pooled)
mi_coefs = [pooled_results.loc['coef_education', 'estimate'],
            pooled_results.loc['coef_health', 'estimate']]

comparison = pd.DataFrame({
    'Method': ['Complete Case', 'Mean Imputation', 'Single MICE', 'Multiple Imputation'],
    'Education Coef': [cc_coefs[0], mean_coefs[0], single_mice_coefs[0], mi_coefs[0]],
    'Health Coef': [cc_coefs[1], mean_coefs[1], single_mice_coefs[1], mi_coefs[1]],
    'Sample Size': [len(df_cc), len(df_full), len(df_full), len(df_full)]
})

print("\nComparison of Coefficients Across Methods:")
print(comparison.round(4))

# ============================================
# STEP 8: Summary Statistics
# ============================================

print("\n" + "="*60)
print("STEP 6: Summary Statistics")
print("="*60)

summary_stats = pd.DataFrame({
    'Original': df_full.mean(),
    'Complete Case': df_cc.mean(),
    'Mean Imputed': df_mean.mean()
})

# Multiple imputation mean (average across all imputations)
mi_means = []
for dataset in imputed_datasets:
    mi_means.append(dataset.mean())
mi_summary = pd.DataFrame(mi_means).mean()
summary_stats['Multiple Imputation'] = mi_summary

# Standard errors for MI (using Rubin's rules)
mi_se = []
for col in df_full.columns:
    if col in pooled_results.index:
        mi_se.append(pooled_results.loc[col, 'se'])
    else:
        mi_se.append(np.nan)
summary_stats.loc['Standard Error'] = mi_se

print("\nVariable Means Across Methods:")
print(summary_stats.round(2))

# ============================================
# STEP 9: Final Summary
# ============================================

print("\n" + "="*60)
print("STEP 7: Final Summary and Recommendations")
print("="*60)

print(f"""
ANALYSIS SUMMARY
================
- Dataset: n={len(df_full)}, p={df_full.shape[1]}
- Missingness: {df_full.isnull().sum().sum() / (len(df_full) * df_full.shape[1]):.1%}
- Imputations: m={m}
- Method: Bayesian Ridge MICE with posterior sampling

KEY FINDINGS
============
- Education coefficient: {mi_coefs[0]:.3f} (95% CI: {pooled_results.loc['coef_education', 'ci_lower']:.3f}
- Health coefficient: {mi_coefs[1]:.3f} (95% CI: {pooled_results.loc['coef_health', 'ci_lower']:.3f} to {p
- Fraction of Missing Information: {pooled_results['fmi'].mean():.3f} (average)

DIAGNOSTICS
===========
- Convergence: ✓ (max_iter={imputer.max_iter}, n_iter={imputer.n_iter_ if hasattr(imputer, 'n_iter_') els
- FMI < 0.3: {'✓' if all(pooled_results['fmi'] < 0.3) else '✗'} (indicates low missing information)
- Between-variance < Within-variance: {'✓' if all(pooled_results['between_var'] < pooled_results['within_

RECOMMENDATIONS
===============
1. Use Multiple Imputation when analyzing relationships with income
2. Report pooled estimates with proper confidence intervals
3. Include FMI values to indicate uncertainty from missingness
4. For prediction only, consider single imputation to save computation
""")

# Save all results to CSV
pooled_results.to_csv("pooled_regression_results.csv")
comparison.to_csv("method_comparison.csv", index=False)
summary_stats.to_csv("summary_statistics.csv")

print("\nFiles saved:")
print("- pooled_regression_results.csv")
print("- method_comparison.csv") 
print("- summary_statistics.csv")
print("- missingness_analysis.png")
print("- pooled_results_analysis.png")
print("- imputed_distributions.png")

11. Reproducibility

  • Python: 3.11+
  • scikit‑learn: 1.5+
  • statsmodels: 0.14+
  • numpy/pandas: recent stable versions
  • Random seeds fixed when possible; report package versions used.

12. External Resources

Books and papers

  • Rubin DB (1987). Multiple Imputation for Nonresponse in Surveys. Wiley.
  • van Buuren S (2018). Flexible Imputation of Missing Data. Chapman & Hall/CRC.
  • White IR, Royston P, Wood AM (2011). Multiple imputation using chained equations: Issues and guidance.

Python documentation

Reporting guidelines (clinical)

13. FAQ

Q1: How many imputations (m)?

Modern practice: 20–100. Better to guide m by FMI: increase m until CIs and estimates stabilize. Relative efficiency increases with m and decreases with FMI.

Q2: Can I use MI for prediction?

Yes. For pure prediction you can average predictions over m models or use deterministic single imputation if you don’t need uncertainty. Keep predictive vs inferential goals separate.

Q3: Which variables to include in the imputation model?

All variables used in the final analysis + auxiliaries that predict missingness/values, including the outcome when appropriate.

Q4: Categorical variables?

Avoid “round back”. During imputation use one‑hot and inverse_transform afterwards, or methods that natively handle categoricals (mice in R, forests).

Q5: What if data are MNAR?

Run sensitivity analyses (pattern‑mixture with delta‑adjustment or selection models). Interpret standard MI with caution.

Q6: How do I check convergence?

Check n_iter_ and stability of imputed statistics. R mice provides native trace plots; in Python implement ad‑hoc checks.

Q7: How to report MI in publications?

Report m, software, method, assumptions (MAR), diagnostics, applied constraints, and pooled results (estimate, SE, CI, FMI, df). Add notes on any sensitivity analyses.

Appendix — MI Reporting Checklist

  • Context and objective
    • Research question and goal (inference vs prediction)
    • Population/dataset, period, inclusion/exclusion criteria
  • Missing data
    • Percent missing by variable and overall
    • Hypothesized mechanism (MCAR/MAR/MNAR) with rationale and evidence
    • Auxiliary variables included in imputation
  • Imputation strategy
    • Software and version (e.g., scikit‑learn 1.5, statsmodels 0.14)
    • Method (e.g., MICE/IterativeImputer, estimator) and key settings
      • m (number of imputations) and selection criterion
      • max_iter, sample_posterior, n_nearest_features
      • Applied constraints (min_value, max_value) and transforms (e.g., log)
    • Handling categoricals (one‑hot, inverse_transform, or native methods)
  • Diagnostics
    • Evidence of convergence or stability (n_iter_, trace/summary of imputed stats)
    • Compare observed vs imputed distributions
    • Clinical/subject‑matter plausibility and coherence checks
  • Analysis on imputed datasets
    • Models fitted on each of the m datasets
    • How variances/SEs were computed per fit (e.g., statsmodels OLS)
  • Pooling
    • Rubin’s rules applied, report df (Barnard–Rubin) and FMI
    • Final estimates with SE, CI and, if relevant, p‑values
  • Sensitivity analysis
    • Strategies for MNAR or robustness (e.g., pattern‑mixture with delta‑adjustment)
    • Impact on primary results
  • Reproducibility
    • Package versions, seeds, relevant hardware
    • Link to code/notebook or repository
  • Final reporting
    • Clear statement of assumptions (MAR vs MNAR)
    • Known limitations and potential residual biases
    • Practical implications and usage recommendations
Edvard Munch Melancholy

Edvard Munch and Illness: The Ailing Body, the Soul Aflame

Posted on by Michele Danilo Pierri

“Sickness, madness and death were the dark angels that had accompanied me since childhood.”

Edvard Munch, Diaries, 1890

In the cold north of Europe, between the late nineteenth and early twentieth centuries, a Norwegian painter transformed his pain into vision. Not abstract pain, but one rooted in flesh, lungs, and nerves: a pain named tuberculosis, grief, existential anguish.

Edvard Munch did not paint illness as an external theme: he lived it, breathed it, embodied it. His art became an autopsy of the soul, a clinical specimen of the invisible.

For physicians (accustomed to reading the body through signs, symptoms, and examinations) Munch’s work offers a unique field of observation. Not for anatomical precision (which he deliberately neglects), but for the phenomenological representation of suffering.

Here, illness is not merely pathology but lived experience, laden with emotional, symbolic, even metaphysical meaning.

Key takeaways

Munch’s illness imagery foregrounds experience over diagnosis and remains clinically relevant.

“The Sick Child”, “The Scream”, and “Madonna” trace grief, anxiety, and eros‑thanatos with striking clarity.

“Self‑Portrait with the Spanish Flu” offers a rare, direct view of post‑viral exhaustion.

Serial reworking of motifs functions like a visual follow‑up, paralleling clinical trajectories.

1) The Sick Child: Tuberculosis and Childhood Trauma

Munch - the sick child

In 1868, when Edvard was only five, his mother Laura died of tuberculosis. Eight years later, his fifteen-year-old sister Sophie succumbed to the same disease.

The young Munch witnessed her slow decline: fever, hemoptysis, waxy pallor, breathing that grew increasingly labored. The bedroom became a theater of death.

The White Plague: Tubercolosis and Childhood Trauma

In nineteenth-century Europe, tuberculosis was not an isolated tragedy but a mass epidemic. Known as “the white plague,” it accounted for nearly one in four deaths across the continent, with mortality rates reaching 400 per 100,000 in urban areas like Oslo.

Children and young adults were particularly vulnerable: pulmonary tuberculosis in adolescents carried a mortality approaching 50% before the advent of antibiotics.

The disease was not only medical but deeply social—associated with poverty, overcrowding, and the rapid industrialization that defined the era.

For families like the Munchs, tuberculosis was an almost inescapable presence. The household became a site of contagion and vigil, where the slow consumption of the body unfolded over months or years.

There was no effective treatment: rest, fresh air, and cod liver oil were the main prescriptions, while the coughing, night sweats, and hemoptysis continued relentlessly. The emotional toll was compounded by the sense of inevitability: once the diagnosis was made, survival was uncertain at best.

Visual Analysis

This episode marks the genesis of The Sick Child (1886, with a more mature version from 1896 held at the Munchmuseet in Oslo). The work is painted in oil on canvas, with dense, textured brushstrokes in the later version that evoke the physical weight of pain.

The bed sits at the center of the composition, tilted diagonally, creating a spatial tension that denies the viewer visual rest.

Sophie lies supine, her pale face turned toward a window through which cold, almost lunar light enters. Her hands are crossed on her chest, a pose reminiscent of the deceased in nineteenth-century funeral portraits.

Beside the bed, a second figure, likely her younger sister Inger, bends forward, her face hidden in her arms. Her black dress merges with the shadow of the room, suggesting that mourning is not merely an event but an atmosphere enveloping everything.

The floor is bare, the walls stripped of ornament: there is no religious consolation, no visible medical intervention. Only the looming presence of the end.

Clinical Relevance

For the contemporary physician accustomed to viewing tuberculosis as a curable disease this scene may appear remote. Yet it remains profoundly relevant in its existential dimension.

Munch does not document Koch’s bacilli; he documents the wait for death, the void that opens around the sick person’s bed, the solitude of the dying body.

Museum link: Edvard Munch, Public domain, via Wikimedia Commons – Munchmuseet

2) Self-Portrait with Spanish Flu: Post-Viral Exhaustion Visualized

Munch Self Portrait with the Spanish Flu

For a clinically unambiguous self‑image, Edvard Munch’s Self‑Portrait with the Spanish Flu (1919) holds a truly pivotal place in art history. In this haunting painting, Munch is depicted wearing a dressing gown, looking gaunt and visibly weakened, with an unmade bed looming behind him as a stark reminder of his recent illness.

The greenish pallor of his skin, the deeply hollowed eyes, and the slightly parted mouth all vividly suggest the profound exhaustion and lingering effects of a severe viral infection.

While his body remains almost motionless, there is an unsettling sense that the world around him continues to vibrate and pulse with restless energy.

Technical execution

The work is executed in oil on canvas with quick, vibrant brushstrokes and a restricted palette of acid yellows, earthy greens, and browns. The forms are not naturalistic: Munch deforms in order to express the subjective experience of illness.

The background is a claustral interior, not a winter landscape. The blanket replaces the nude—here the body is covered and fragile, not displayed. The air in the room feels dense and stagnant

Historical Resonance: Spanish Flu and Long COVID

In the aftermath of COVID-19, this 1919 portrait resonates with striking immediacy. The Spanish flu infected one-third of humanity and left survivors with prolonged fatigue, cognitive fog, and a sense of bodily disconnection—symptoms now recognized as post-acute sequelae of viral infections (PASC).

Munch’s visual testimony predates our contemporary language for “long COVID”by a century, offering clinicians a phenomenological map of what 10–30% of COVID survivors still experience: not recovery, but a protracted negotiation with the afterlife of illness.

Clinical Lens

In this self-portrait, Munch depicts himself not as a subject of diagnosis but as a person undergoing a profound erosion of bodily integrity.

The gaunt posture, hollow eyes, and static environment convey the disorientation typical of severe post-viral states.

Clinically, the painting invites reflection on the subjective dimension of recovery: how illness reshapes one’s sense of time, agency, and embodiment. It also reminds clinicians that post-viral fatigue is not only biological but deeply narrative—an experience that destabilizes identity before it restores it.

Museum links: Edvard Munch, Public domain, via Wikimedia Commons • Munchmuseet

3) The Scream: Anxiety as Physiological Experience

Munch The Scream

While tuberculosis marked his body, anguish dominated his psyche. Here Munch makes an extraordinary conceptual leap: he transforms an interior state into a universal image.

The Scream (1893), of which four versions exist (two paintings, two pastels), is not about mental illness in a nosographic sense. It portrays neither a psychotic patient nor a diagnosable panic attack.

Yet no other image has expressed with such force the sensation of the self dissolving.

Composizional Brilliance

The composition is brilliant in its simplicity: an androgynous figure with an elongated skull and hands pressed against its ears stands at the center of a bridge crossing a fjord.

The sky ignites with undulating streaks of red, orange, and yellow, while the landscape, trees, water, mountains, liquefies into sinuous curves, as if seen through distorting glass.

Color does not describe reality but reality’s effect on the nervous system.

The 1893 pastel version (now at the Munchmuseet) is particularly intense: pigments are applied with quick, almost feverish gestures. The lines do not delineate forms but vibrations.

The figure’s face has no defined human features: eyes, mouth, and nose are reduced to holes, as if the body were losing its material consistency.

Clilnical Lens

Rather than illustrating a psychiatric label, The Scream visualizes the physiology of alarm: derealization, autonomic surge, and overwhelming allostatic load.

The melting landscape echoes the perceptual distortions seen in acute anxiety or panic states, where reality loses its boundaries.

For clinicians, the work challenges the reduction of anxiety to symptom checklists, urging a more phenomenological understanding of what it feels like when the self approaches fragmentation.

It is a reminder that distress can reach clinical thresholds without ever fitting neatly into diagnostic categories.

Museum link: Edvard Munch, Public domain, via Wikimedia Commons Munchmuseet

4) Madonna: Eros, Thanatos, and Medical Anxieties

Munch "Madonna"

In Edvard Munch’s artistic work, the themes of eros and thanatos are intricately intertwined, creating a complex exploration of life and death forces.

In his piece “Madonna”, a nude female figure is depicted in an ecstatic, almost rapturous pose, reclining on a deep red bed. The use of the traditional Marian title “Madonna” is deliberately subverted here, blending together elements of ecstatic experience, fertility, and an undercurrent of danger.

In the lithograph versions of this work, the addition of the iconic sperm border makes the themes of desire and creation much more explicit and visually pronounced.

Medical Context: The Syphilis Epidemic

In the 1890s, syphilis was epidemic across Europe, with urban prevalence reaching 5–15%. Known as “the great imitator,” it manifested in stages (chancre, rash, neurological collapse) and carried profound moral stigma.

With no cure until Salvarsan (1910), the disease fueled anxieties about hereditary “degeneration” and sexual transgression.

Munch’s “Madonna” captures this era’s ambivalence: the female body as both sacred and dangerous, desire as life force and contagion. The red halo evokes fever as much as divinity.

For clinicians today, it remains a potent reminder that bodies are never purely medical—they are always morally interpreted, particularly in domains of sexuality and reproduction.

Visual Symbolism

The work showcases chromatic symbolism at its finest. The red is not realistic but psychological—it invades the space, suffocates, envelops.

The woman’s skin is milky white, almost translucent, evoking a fragile, threatened purity.

In the lithographic version, Munch added an inscription: “You are the woman I love—and who destroys me.”

Clinical Lens

Munch’s Madonna collapses binaries—healthy/sick, sacred/profane, desire/danger—exposing the body as a site of vulnerability and ambiguity.

The sensual posture, combined with undertones of threat and mortality, reflects the anxieties of his era around sexuality, syphilis, and hereditary “madness.”

Clinically, the painting opens a space to consider desire as a dimension of care, not merely a risk factor. It foregrounds the emotional and existential tensions that patients carry, reminding clinicians that the body is never purely medical but always symbolic, relational, and morally interpreted.

Museum links: Edvard Munch, Public domain, via Wikimedia Commons – Munchmuseet

5) Serial Imagery as Clinical Follow-Up

“ I do not paint what I see, but what I have seen.”

Munch returned obsessively to the same motifs across decades. “The Sick Child” exists in at least six painted versions (1885–1927), each subtly different in tone, texture, and emotional intensity.

This was not repetition but longitudinal observation—a visual follow-up of unresolved grief.

Clinically, this mirrors how we track chronic conditions over time: not as single events but as trajectories requiring serial engagement.

The 1885 version is raw, almost violent; the 1896 version more resolved but still aching; the 1927 version muted, resigned.

What changed was not the event (Sophie’s death in 1877) but Munch’s relationship to the memory—the slow, non-linear work of integration.

This has implications for practice. Many conditions—chronic pain, bereavement, PTSD—do not follow linear paths toward “closure.” They require repeated visits, adjusted narratives, sustained witness.

Munch’s serial practice models this: care as persistent presence, not cure. In clinic, not every pain is resolved, but every pain deserves recognition.

Sometimes the highest form of care is simply to remain and witness.

Conclusion: Munch, medicine, and narrative care

As medicine advances technologically and molecularly, Munch restores the human, narrative, embodied dimension.

Narrative medicine has shown how illness narratives improve listening, diagnosis, adherence, and alliance. Munch offers no solutions—only visual questions that re‑center the sickbed as a site of truth.

“My painting is a confession. And every confession is an act of healing—not for the speaker, but for the listener.”

 Explore More: Art, Medicine, and the Humanities

Image use and attributions

Image via Wikimedia Commons Public Domain: “Melancholy”, “The Sick Child”, “Self-Portrait with the Spanish Flu”, “The Scream”, “Madonna”

Munchmuseet

FAQ: Edvard Munch and illness in art

What illnesses appear in Edvard Munch’s art?

Tuberculosis, anxiety, depression, and the 1918–19 influenza appear as biographical and symbolic forces.

Why is The Sick Child important?

It reframes pediatric illness and grief as lived experience, not medical spectacle.

What does The Scream represent in mental health terms?

A powerful visualization of anxiety, derealization, and autonomic overload without pathologizing the subject.

Did Munch depict his own illness?

Yes, most explicitly in Self‑Portrait with the Spanish Flu (1919).

How is Munch relevant to clinicians today?

His images support reflective practice, empathy, and narrative competence in care.

Villagers and children in period clothing gaze upward as dozens of muted red, gold, blue, and cream balloons drift above a cobbled old European street in a warm, painterly historical scene.

Softmax

Posted on by Michele Danilo Pierri

Introduction

The softmax function is essential in mathematics and machine learning.

It transforms a vector of real numbers into a probability distribution.

Put simply, it converts a set of numbers into probabilities.

To understand how it works, let’s look at a list of risk values for some patients:

PatientsSurgical Risk
Patient A2
Patient B1
Patient C0

When we apply the softmax function to this series of numbers [2,1,0] we get:

PatientsSurgical Risk
Patient A66% (0.66)
Patient B24% (0.24)
Patient C9% (0.9)

Note that the resulting probability values from the softmax function always sum to 1 (or 100%).

How does it work?

The Softmax transformation consists of two key operations: exponentiation and normalization.

During exponentiation, we calculate e^x for each number. This amplifies larger numbers while reducing smaller ones in the series.

In the normalization step, we sum all the numbers and divide each by that total. This produces values between 0 and 1 that always sum to 1.

These two steps together create our final probability distribution.

The graphs illustrate how Softmax transforms numerical values into probabilities, with the resulting probabilities always summing to one.

From initial scores to final probability bar graphs

Probability function with softmax

Mathematical Formula for Softmax

For a vector z = [z₁, z₂, z₃…zₙ], the Softmax formula is:

\text{softmax}(z_i) = \frac{e^{z_i}}{\sum_{j=1}^{n} e^{z_j}}

Where: – z_i is the i-th element of the input vector z, – e^{z_i} is the exponential of z_i, – \sum_{j=1}^{n} e^{z_j} is the sum of the exponentials of all elements in the vector z, – n is the total number of elements in the vector z. This formula ensures that each output value lies between 0 and 1, and the sum of all outputs equals 1.

​

Graphical Examples of Softmax

For two classes, the Softmax function simplifies to a sigmoid function, where the first class has probability p1 and the second class has probability 1-p1 (since probabilities must sum to 1). As one class’s probability increases, the other’s must decrease proportionally.

Softmax function for two classes

For three classes, we can visualize the function using a three-dimensional graph. When we treat the first two classes as variables and fix the third as a constant, the graph displays the probabilities of the first two classes. The third class’s probability is then calculated as 1 minus the sum of the first two class probabilities.

Softmax function for three classes

Numerical Saturation and Normalization by Maximum

A key challenge when applying the Softmax function occurs with extremely large or small z values.

These extreme values can cause overflow or underflow—situations where numbers become too large or too small for a computer to represent accurately.

Consider a vector z with large numbers:

z=[1000, 1001, 1002]

​Calculating e^{1000}, e^{1001}, and e^{1002} for Softmax would produce enormous numbers that cause overflow.

To solve this, we can normalize the vector. The new vector z’ then has much smaller values:

z’ = [1000 – 1002, 1001 – 1002, 1002 – 1002] = [-2, -1, 0]

​This normalization gives us manageable exponential values:

e^{-2} \approx 0.135, \quad e^{-1} \approx 0.368, \quad e^{0} = 1

​Finally, we calculate the sum of exponentials and apply the softmax function:

\text{Sum} = 0.135 + 0.368 + 1 = 1.503


\text{Softmax} = \left[\frac{0.135}{1.503}, \frac{0.368}{1.503}, \frac{1}{1.503}\right] \approx [0.09, 0.24, 0.67]

​

This same approach works for very small z values.

For both extremely large and small values, we can use a modified formula:

\text{softmax}(z_i) = \frac{e^{z_i - \max(z)}}{\sum_{j=1}^{n} e^{z_j - \max(z)}}

​

Softmax with Python

Let’s explore how to implement the Softmax function in Python, covering both single vector applications and matrix operations.

# Example of applying softmax to a NumPy array

import numpy as np

def softmax(z):
    # Subtract maximum value to prevent numerical overflow
    z = z - np.max(z)
    exp_z = np.exp(z)
    return exp_z / np.sum(exp_z)

# Example usage
z = np.array([2.0, 1.0, 0.1])
print("Input:", z)
print("Softmax Output:", softmax(z))

# Matrix application example

def softmax_batch(z):
    # Subtract the maximum along axis 1 (per row)
    z = z - np.max(z, axis=1, keepdims=True)
    exp_z = np.exp(z)
    return exp_z / np.sum(exp_z, axis=1, keepdims=True)

# Usage example
z_batch = np.array([[2.0, 1.0, 0.1], [1.0, 2.0, 3.0]])
print("Input Batch:\n", z_batch)
print("Softmax Output Batch:\n", softmax_batch(z_batch))

​

Applications

The Softmax function has three main applications:

In Multiclass Classification (Machine Learning): It converts raw scores into a probability distribution across possible classes

In Neural Networks: It serves as the activation function in the output layer

In Reinforcement Learning: It transforms action scores into selection probabilities

Conclusion

The Softmax function plays a vital role in machine learning, especially for multiclass classification tasks. Its mathematical properties and computational efficiency have made it indispensable in neural networks and predictive models. Proper implementation and careful handling of numerical challenges are key to achieving optimal results.

A physician in a pale coat stands at a forked mountain path at sunset, holding a glowing compass between two wooden signs marked “TREAT NONE” and “TREAT ALL,” symbolizing a clinical decision between opposing treatment strategies.

Decision Curve Analysis

Posted on by Michele Danilo Pierri

Introduction

Decision Curve Analysis (DCA) is a powerful tool for evaluating the clinical utility of predictive models and diagnostic tests. Unlike traditional metrics like AUC or calibration, DCA focuses on what truly matters in practice: whether using a model leads to better decisions and outcomes.

DCA was introduced by Vickers and Elkin in 2006 to address the limitations of conventional performance metrics. It calculates the net benefit of using a model across a range of threshold probabilities—helping clinicians decide whether a model improves decision-making compared to default strategies like treating all or no patients.

In this guide, we’ll walk through the core concepts of DCA, how to interpret decision curves, calculate net benefit, and apply these insights to real-world clinical scenarios.

Case Example: Prostate Cancer Biopsy

Imagine a patient who presents with elevated PSA (prostate-specific antigen) levels during routine screening. A predictive model is then applied to estimate the patient’s individualized risk of hosting high-grade prostate cancer. Rather than proceeding with a biopsy for every single patient who has elevated PSA—which would result in many unnecessary and invasive procedures—Decision Curve Analysis helps clinicians determine whether incorporating the predictive model into their decision-making process actually leads to better clinical outcomes. Specifically, it evaluates whether the model results in fewer unnecessary biopsies being performed on patients who do not have aggressive disease, while simultaneously achieving more accurate and timely identification of those patients who do have aggressive cancers that require intervention.

Key concepts

  • Net benefit: balances true positives and false positives using the threshold probability.
  • Threshold probability (pt): the minimum predicted risk at which intervention is justified.
  • Exchange rate: pt/(1−pt), the implied trade-off between one false negative and false positives.
  • Strategies: Treat all, Treat none, and Model-based.

Interpreting a decision curve

Decision curve

Net benefit across clinically relevant thresholds. Compare the model to Treat all and Treat none.

Step-by-Step Interpretation:

  1. Higher Curve = Greater Benefit: The model with the highest curve offers the most clinical value.
  2. Preferences Matter: Some patients prioritize avoiding disease; others fear unnecessary procedures.
  3. Threshold Probability Is the Decision Point: It defines when intervention becomes justified.
  4. Net Benefit Is Like Net Profit: It’s the clinical “gain” after accounting for harms.
  5. Can Be Expressed as Interventions Avoided: Useful for communicating impact in practical terms.
Model comparison with decision curve

Decision curves can help to compare models:

  • Higher curve means greater clinical utility at that threshold.
  • Focus on the clinically plausible threshold range for the condition and intervention.
  • Compare your model against Treat all and Treat none to contextualize gains.
  • Translate net benefit into people terms when communicating with clinicians.

Net benefit formula

\text{Net Benefit} \,=\, \frac{\text{TP}}{n} \, - \, \frac{\text{FP}}{n} \cdot \frac{p_t}{1-p_t}

Where TP and FP are counts on a sample of size n, and ptpt​ is the threshold probability. The factor pt1−pt1−pt​pt​​ is the exchange rate between harms of false positives and benefits of true positives.

Equivalent impact metrics per 100 patients:

True-positive equivalents per 100 = 100×NB100×NB​

Interventions avoided per 100 = 100 \times \text{NB} \times \tfrac{1-p_t}{p_t}

Clinical Impact Plot

Clinical impact plot

This graph shows two lines, both answering a practical question:

“If I use this prediction model on 100 patients, what actually happens?”

  • Green line (True Positives): The number of sick patients who are correctly identified and would receive treatment. → These are the people who benefit from using the model.
  • Orange line (False Positives): The number of healthy patients who are mistakenly flagged as high-risk and would receive unnecessary treatment. → These are the people who might be harmed (side effects, anxiety, cost) due to overuse.

Both numbers change as you adjust the threshold (the risk level at which you decide to treat). For example:

  • At a low threshold (e.g., 5%), you treat almost everyone → you catch more sick people (high green), but also treat many healthy ones (high orange).
  • At a high threshold (e.g., 40%), you treat only the highest-risk patients → fewer healthy people are treated (low orange), but you also miss some sick patients (low green).

Why this matters:

Doctors don’t think in “net benefit” or “AUC”—they think in people.

This plot translates abstract model performance into real human outcomes:

“At a 15% risk threshold, I’ll help about 22 patients—but 18 others will get treatment they don’t need.”

This helps clinicians choose a threshold that feels right for their patients, their setting, and the seriousness of the treatment.

Calibration Plot

Calibration plot

This graph checks whether the model’s predicted risks match reality.

X-axis: What the model says the risk is (e.g., “20% chance of disease”).

Y-axis: What actually happened in patients with that predicted risk (e.g., did 20% of them really get sick?).

The dotted diagonal line represents perfect trust: if the model says 30%, then 30% of those patients get sick.

The blue markers show how your model actually performed.

You’ll also see vertical dashed lines at common decision points (10%, 20%, 30%)—these are the thresholds doctors might use to decide who to treat.

Why this matters:

Decision Curve Analysis assumes your model’s probabilities are honest.

But if your model is overconfident (e.g., says “10% risk” but 25% actually get sick), then:

Using a 10% treatment threshold would mean treating far too few people.

Your DCA results could look good—but in reality, you’re missing many at-risk patients.

Similarly, if the model is underconfident (says “40%” but only 15% get sick), you’d overtreat healthy people.

Recommended analysis workflow

Split or use external validation. Freeze indices for reproducibility.

Fit model and obtain predicted probabilities on the test set.

Calibrate if needed (Platt scaling or isotonic regression with proper cross-validation).

Compute net benefit across a clinically justified threshold range.

Quantify uncertainty with bootstrap confidence intervals.

Report clinical impact and, for model comparisons, differences in net benefit at pre-specified thresholds with CIs.

Python implementation

The code below generates four figures and adds bootstrap CIs for net benefit at selected thresholds. Save the figures and embed them here.

Step 1 — Data setup

Use a small, self‑contained example to demonstrate the full workflow. In practice, replace the simulated dataset with your real features X and labels y.

Step 2 — Train/test split

Split the data into training and test sets. Train on the former and evaluate on the latter to obtain unbiased estimates of performance and decision utility.

Step 3 — Base model

Fit a simple baseline model (logistic regression here) to produce predicted probabilities for the outcome of interest. Any probabilistic classifier can be used as long as it outputs calibrated probabilities.

Step 4 — Probability calibration (if needed)

If probabilities are miscalibrated around clinically relevant thresholds, calibrate them using Platt scaling or isotonic regression with proper cross‑validation. DCA relies on well‑calibrated probabilities near the decision thresholds.

Step 5 — Decision curve calculation

For a clinically justified range of threshold probabilities, compute net benefit for three strategies: Treat none, Treat all, and the Model. This shows the clinical utility of using the model versus simple baselines.

Step 6 — Uncertainty via bootstrap

Use bootstrap resampling of the test set to obtain confidence intervals for net benefit curves and key threshold points. Report CIs for both single‑model curves and differences between models.

Step 7 — Figures and saving

Generate four outputs: decision curve, model comparison, clinical impact plot, and calibration with threshold markers. Save figures to disk and embed them in this page for reporting and discussion.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.calibration import calibration_curve
from sklearn.isotonic import IsotonicRegression

# -------------------------------
# 1) Simulated dataset (replace with your data)
# -------------------------------
RNG_SEED = 42
np.random.seed(RNG_SEED)
X, y = make_classification(
    n_samples=2000,
    n_features=10,
    n_informative=6,
    n_redundant=4,
    weights=[0.7, 0.3],    # ~30% event rate
    flip_y=0.05,
    random_state=RNG_SEED,
)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=RNG_SEED
)

# -------------------------------
# 2) Base model
# -------------------------------
model = LogisticRegression(max_iter=1000)
model.fit(X_train, y_train)
y_proba_raw = model.predict_proba(X_test)[:, 1]

# -------------------------------
# 3) Optional post-hoc calibration (isotonic)
#    Fit on train via CV in real studies. Here, a simple holdout for demo.
# -------------------------------
# Map raw scores from (X_train) via isotonic then apply to (X_test) would require CV.
# For demonstration, we calibrate on test predictions vs. test labels cautiously.
ir = IsotonicRegression(out_of_bounds='clip')
ir.fit(y_proba_raw, y_test)
y_proba = ir.transform(y_proba_raw)

# -------------------------------
# Helpers
# -------------------------------

def net_benefit(y_true, y_prob, thresholds):
    y_true = np.asarray(y_true)
    y_prob = np.asarray(y_prob)
    n = len(y_true)
    nb = []
    for t in thresholds:
        pred = (y_prob >= t).astype(int)
        tp = np.sum((y_true == 1) & (pred == 1))
        fp = np.sum((y_true == 0) & (pred == 1))
        nb_val = (tp / n) - (fp / n) * (t / (1 - t))
        nb.append(nb_val)
    return np.array(nb)


def net_benefit_at_threshold(y_true, y_prob, t):
    y_true = np.asarray(y_true)
    y_prob = np.asarray(y_prob)
    n = len(y_true)
    pred = (y_prob >= t).astype(int)
    tp = np.sum((y_true == 1) & (pred == 1))
    fp = np.sum((y_true == 0) & (pred == 1))
    return (tp / n) - (fp / n) * (t / (1 - t))


def clinical_impact(y_true, y_prob, thresholds, per=100):
    y_true = np.asarray(y_true)
    y_prob = np.asarray(y_prob)
    n = len(y_true)
    tp_list, fp_list = [], []
    for t in thresholds:
        pred = (y_prob >= t).astype(int)
        tp = np.sum((y_true == 1) & (pred == 1))
        fp = np.sum((y_true == 0) & (pred == 1))
        tp_list.append(tp * per / n)
        fp_list.append(fp * per / n)
    return np.array(tp_list), np.array(fp_list)


def bootstrap_ci_nb(y_true, y_prob, t, B=1000, rng=None):
    rng = np.random.default_rng(rng)
    n = len(y_true)
    vals = np.empty(B)
    for b in range(B):
        idx = rng.integers(0, n, n)
        vals[b] = net_benefit_at_threshold(y_true[idx], y_prob[idx], t)
    return np.quantile(vals, [0.025, 0.5, 0.975])

# -------------------------------
# 4) Threshold range and baselines
# -------------------------------
thresholds = np.linspace(0.05, 0.50, 100)  # adjust to clinical range
prevalence = np.mean(y_test)
nb_treat_all = prevalence - (1 - prevalence) * (thresholds / (1 - thresholds))
nb_treat_none = np.zeros_like(thresholds)

# -------------------------------
# 5) Figures
# -------------------------------
# Plot 1: DCA
nb_model = net_benefit(y_test, y_proba, thresholds)
plt.figure(figsize=(8, 5.5))
plt.plot(thresholds, nb_model, label='Prediction model', color='steelblue', lw=2.5)
plt.plot(thresholds, nb_treat_all, label='Treat all', color='crimson', ls='--', lw=2)
plt.plot(thresholds, nb_treat_none, label='Treat none', color='gray', ls=':', lw=2)
plt.xlabel('Threshold probability')
plt.ylabel('Net Benefit')
plt.title('Decision Curve Analysis')
plt.legend(frameon=False)
plt.grid(True, ls='--', alpha=0.6)
plt.xlim([thresholds.min(), thresholds.max()])
plt.tight_layout()
plt.savefig('figure_dca.png', dpi=300, bbox_inches='tight')

# Plot 2: Clinical Impact (per 100 patients)
thresholds_ci = np.linspace(0.05, 0.50, 60)
tp_vals, fp_vals = clinical_impact(y_test, y_proba, thresholds_ci, per=100)
plt.figure(figsize=(8, 5.5))
plt.plot(thresholds_ci, tp_vals, label='True positives per 100', color='green', lw=2.5)
plt.plot(thresholds_ci, fp_vals, label='False positives per 100', color='orange', lw=2.5)
plt.xlabel('Threshold probability')
plt.ylabel('Number of patients (per 100)')
plt.title('Clinical Impact Plot')
plt.legend(frameon=False)
plt.grid(True, ls='--', alpha=0.6)
plt.xlim([thresholds_ci.min(), thresholds_ci.max()])
plt.tight_layout()
plt.savefig('figure_clinical_impact.png', dpi=300, bbox_inches='tight')

# Plot 3: Model comparison (simulate weaker model)
y_proba_weak = np.clip(y_proba + np.random.normal(0, 0.15, size=y_proba.shape), 0, 1)
nb_strong = net_benefit(y_test, y_proba, thresholds)
nb_weak = net_benefit(y_test, y_proba_weak, thresholds)
plt.figure(figsize=(8, 5.5))
plt.plot(thresholds, nb_strong, label='Enhanced model', color='steelblue', lw=2.5)
plt.plot(thresholds, nb_weak, label='Basic model', color='purple', ls='-.', lw=2.5)
plt.plot(thresholds, nb_treat_all, label='Treat all', color='crimson', ls='--', lw=2)
plt.plot(thresholds, nb_treat_none, label='Treat none', color='gray', ls=':', lw=2)
plt.xlabel('Threshold probability')
plt.ylabel('Net Benefit')
plt.title('Model Comparison via Net Benefit')
plt.legend(frameon=False)
plt.grid(True, ls='--', alpha=0.6)
plt.xlim([thresholds.min(), thresholds.max()])
plt.tight_layout()
plt.savefig('figure_model_comparison.png', dpi=300, bbox_inches='tight')

# Plot 4: Calibration with threshold markers
frac_pos, mean_pred = calibration_curve(y_test, y_proba, n_bins=10, strategy='quantile')
plt.figure(figsize=(7, 7))
plt.plot(mean_pred, frac_pos, 's-', color='darkblue', label='Model', lw=2, ms=6)
plt.plot([0, 1], [0, 1], 'k:', label='Perfect calibration', lw=1.5)
for th in [0.1, 0.2, 0.3]:
    plt.axvline(x=th, color='gray', ls='--', alpha=0.7)
    plt.text(th + 0.01, 0.02, f'{int(th*100)}%', rotation=90, color='gray', fontsize=10)
plt.xlabel('Predicted probability')
plt.ylabel('Observed frequency')
plt.title('Calibration Plot (zoom to clinical range as needed)')
plt.legend(frameon=False)
plt.grid(True, ls='--', alpha=0.6)
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.tight_layout()
plt.savefig('figure_calibration.png', dpi=300, bbox_inches='tight')

# -------------------------------
# 6) Example: bootstrap CIs at pre-specified thresholds
# -------------------------------
for t in [0.10, 0.20, 0.30]:
    lo, med, hi = bootstrap_ci_nb(y_test, y_proba, t, B=500, rng=123)
    print(f"Threshold {t:.2f}: NB median {med:.4f} (95% CI {lo:.4f} to {hi:.4f})")

Reporting recommendations

  • State the clinical rationale for the chosen threshold range. For prostate biopsy, 5–30% is commonly discussed in the literature.
  • Provide decision curves with Treat all and Treat none baselines.
  • Include calibration assessment, ideally focusing on the threshold range of interest.
  • Quantify uncertainty with bootstrap CIs for net benefit and for differences between models at pre-specified thresholds.
  • Validate externally when possible; results are population- and prevalence-dependent.

Limitations

  • DCA is not a cost-effectiveness analysis, though it reflects preferences through p_t.
  • Miscalibration near decision thresholds can mislead net benefit.
  • Model utility depends on implementation burden and downstream harms, not just curve separation.

References

  • Vickers AJ, Elkin EB. Decision curve analysis: a novel method for evaluating prediction models. Medical Decision Making, 2006.
  • Vickers AJ, Van Calster B, Steyerberg EW. Net benefit approaches to the evaluation of prediction models. Epidemiology, 2016.
  • Van Calster B, McLernon DJ, van Smeden M, Wynants L, Steyerberg EW. Calibration: the Achilles heel of predictive analytics. BMJ, 2019.
A lone early-20th-century traveler stands on a rocky ridge, holding a compass toward the light while surveying a vast golden mountain valley crossed by winding paths, in a warm, antique painterly style.

Orientation in Dicom

Posted on by Michele Danilo Pierri

Understanding DICOM Coordinate Systems and Image Orientation: Why Your 3D Volume Looks Upside Down

1. Introduction — Why Orientation Matters

Have you ever opened a medical image and found the anatomy upside down or mirrored?

It’s not your viewer’s fault — it’s about geometry.

DICOM files contain not only pixels, but also the mathematical information that tells a viewer where those pixels belong in the patient’s body.

This information — stored in a few special orientation tags — determines whether your 3D reconstruction looks anatomically correct or completely inverted.

In this article, we’ll explore:

  • how DICOM defines spatial orientation,
  • what its key tags actually mean,
  • and how to verify them in Python.

By the end, you’ll understand why one missing minus sign can literally turn a patient upside down.

2. From Pixels to Space — How Medical Images Have Coordinates

When you view a CT slice, you’re looking at a 2D grid of numbers.

But in medicine, every pixel must correspond to a real point in space, measured in millimeters.

To achieve this, DICOM defines a patient-based coordinate system, called LPS:

L (Left) → x-axis positive toward the patient’s left

P (Posterior) → y-axis positive toward the back

S (Superior) → z-axis positive toward the head

So, instead of just rows and columns, every DICOM slice is a plane positioned in 3D, with its own origin, orientation, and scale.

Some research formats, such as NIfTI, use a different convention called RAS (Right–Anterior–Superior), where the X and Y axes are mirrored relative to DICOM’s LPS system.
For clinical DICOM images, however, all coordinates and orientation vectors are defined in the LPS frame, the only one used by PACS viewers and DICOM software.

3. DICOM Tags: How Geometry Is Stored

Every piece of information in a DICOM file is stored as a data element, identified by a tag.

Each data element has four key components:

FieldMeaningExample
Tag4-byte identifier (Group,Element) in hex(0020,0037)
VR (Value Representation)Data type (e.g., DS = Decimal String)DS
VM (Value Multiplicity)How many values (1, 2, 3, 6, …)6
ValueActual data stored as text or binary"1\\0\\0\\0\\-1\\0"

Together, these fields describe everything from patient name to scanner position — but for orientation, three particular tags define where and how each image plane exists in space.

4. The Geometry Trio: IPP, IOP, and PS

These three tags are the geometric foundation of every DICOM image:

TagNameVRVMPurposeExample
(0020,0032)ImagePositionPatient (IPP)DS33D coordinates (x, y, z) of the top-left pixel center (mm). Defines where the plane is."-121.7\\-23.7\\766.7"
(0020,0037)ImageOrientationPatient (IOP)DS6Two unit vectors describing row and column directions in patient coordinates. Defines how the plane is oriented."1\\0\\0\\0\\-1\\0"
(0028,0030)PixelSpacing (PS)DS2Physical distance (mm) between pixel centers along rows and columns. Defines scale."0.625\\0.625"

All coordinates are expressed in millimeters in the LPS frame.

5. How These Tags Define an Image Plane

Each DICOM image is not just a 2D grid of pixels — it’s a plane positioned in the 3D coordinate system of the patient.

To understand where each pixel lies in space, DICOM combines three pieces of information:

  1. ImagePositionPatient (IPP) → the 3D coordinates of the origin (the center of the top-left pixel).
  2. ImageOrientationPatient (IOP) → two unit vectors defining the row and column directions of the image plane.
  3. PixelSpacing (PS) → the physical distance between adjacent pixels, measured in millimeters.

Together, they define a simple but powerful equation that maps pixel indices (i, j) to their physical location (x, y, z) in the patient’s coordinate system (LPS).

Graphic illustration of Dicom spatial concepts

The DICOM Spatial Mapping Formula

According to the DICOM standard (Part 3, Section C.7.6.2.1-1):

P(i,j) = IPP + j · PS[1] · row + i · PS[0] · col

where:

SymbolMeaning
P(i, j)3D coordinates (x, y, z) of pixel (i, j) in the patient’s space
IPPImagePositionPatient — origin of the image plane (mm)
PS[0]PixelSpacing for rows (row spacing). It scales the column direction (col).
PS[1]PixelSpacing for columns (column spacing). It scales the row direction (row).
rowfirst three values of ImageOrientationPatient (direction cosines of image rows)
collast three values of ImageOrientationPatient (direction cosines of image columns)
i, jrow and column indices, starting from (0,0) in the top-left corner

Intuitive interpretation

  • Moving by +1 column (increasing j) shifts you along the row direction (row × PS[1] mm).
  • Moving by +1 row (increasing i) shifts you along the column direction (col × PS[0] mm).
  • The origin (0,0) is at the top-left pixel center, whose absolute coordinates are given by IPP.

The plane normal — the direction in which slices are stacked to form a 3D volume — is defined by the cross product:

normal = row × col

Practical insight

This simple affine relationship is what allows 3D reconstruction software (like 3D Slicer, OsiriX, or Weasis) to rebuild a consistent volume.

However, if the normal vector points in the wrong direction (for example, due to swapped axes or inconsistent slice order), the resulting volume will appear flipped — even though all the pixel data are numerically correct.

6. Example: Reading and Interpreting Real Tag Values

Let’s look at a real-world example taken from an actual DICOM header:

(0020,0032) ImagePositionPatient = -121.7\\-23.7\\766.7
(0020,0037) ImageOrientationPatient = 1\\0\\0\\0\\-1\\0
(0028,0030) PixelSpacing = 0.625\\0.625

From these values we can reconstruct the geometry of a single slice.

Step 1 – Extract the vectors

  • Row direction (first 3 values of IOP): row = [1, 0, 0] → points toward the patient’s left (L).
  • Column direction (last 3 values of IOP): col = [0, -1, 0] → points toward the patient’s posterior (P).
  • Normal vector (cross product): normal = row × col = [0, 0, -1] → points toward the inferior (feet).

This means the slices are physically stacked from superior to inferior (downward) along the patient’s body axis.

Step 2 – Understand the Pixel Spacing

PixelSpacing = [0.625, 0.625]

​These values represent the physical distance (in millimeters) between:

adjacent rows → along the column direction (PS[0]), and

adjacent columns → along the row direction (PS[1]).

So, moving one column to the right shifts the pixel 0.625 mm along row, and moving one row down shifts it 0.625 mm along col.

Step 3 – Compute any pixel’s real-world position

For pixel coordinates (i, j) (where i = row index, j = column index):

P(i,j) = IPP + j · PS[1] · row + i · PS[0] · col

Using the tag values:

P(i,j) = [-121.7, -23.7, 766.7] + j · 0.625 · [1, 0, 0] + i · 0.625 · [0, -1, 0]

This equation allows you to locate any pixel in absolute patient coordinates (LPS).

Step 4 – Analyze the slice orientation

Because the normal vector = [0, 0, -1], the Z-axis decreases as slice numbers increase — meaning that, in 3D, the next slice has a smaller Z value.

If your viewer assumes slices increase along +Z (superior direction), the reconstructed volume will appear upside down.

That’s why understanding the relationship between IOP, IPP, and slice order is essential for correct 3D visualization.

Summary

ConceptDefined byDirectionTypical interpretation
OriginImagePositionPatient(0,0) pixel center3D anchor point of slice
Row directionIOP[0:3]+X (Left)Horizontal axis on image
Column directionIOP[3:6]±Y (Posterior or Anterior, depending on IOP)Vertical axis on image
SpacingPixelSpacingPS[0] rows → along col • PS[1] cols → along rowPhysical scale
Normalrow × col+Z or –Z (depends on orientation)Slice stacking direction

In short, each DICOM slice is a mathematically defined plane in the patient’s body.

By combining ImagePositionPatient, ImageOrientationPatient, and PixelSpacing, you can reconstruct where every pixel lies in millimeter-accurate space — and explain exactly why a 3D volume looks “flipped” when these relationships are misunderstood.

7. Python Example — Read, Analyze, and Validate Orientation

The following script extracts and interprets the geometry of your DICOM files:

import numpy as np, pydicom
from glob import glob

def parse_floats(v):
    s = str(v).replace(',', '\\\\')
    return np.array([float(x) for x in s.split('\\\\') if x], dtype=float)

def read_geometry(ds):
    ipp = parse_floats(ds.ImagePositionPatient)
    iop = parse_floats(ds.ImageOrientationPatient)
    ps  = parse_floats(ds.PixelSpacing)
    row, col = iop[:3], iop[3:]
    row, col = row/np.linalg.norm(row), col/np.linalg.norm(col)
    normal = np.cross(row, col)
    return ipp, row, col, normal, ps

files = sorted(glob("DICOM_STACK/*.dcm"))
d1, d2 = map(pydicom.dcmread, files[:2])

ipp, row, col, normal, ps = read_geometry(d1)
print("IPP:", ipp)
print("Row:", row)
print("Column:", col)
print("Normal:", normal)
print("Pixel Spacing:", ps)

dz = [np.dot](<http://np.dot>)((read_geometry(d2)[0] - ipp), normal)
print("Δ along normal between slice #1 and #2 (mm):", dz)
if dz < 0:
    print("Warning: slices are stacked in the opposite direction.")

This lets you verify:

  • whether row/column vectors are orthogonal;
  • whether slices increase along the expected direction;
  • whether the viewer’s 3D reconstruction should appear upright.

8. Common Pitfalls and How to Avoid Them

❌ Assuming file order = anatomical order→ Always check the Z difference between consecutive ImagePositionPatient values.

❌ Mixing coordinate conventions→ DICOM uses LPS; some research tools use RAS (mirrored X/Y).

❌ Ignoring direction cosines→ The slice order alone doesn’t guarantee correct 3D orientation.

❌ Forgetting to normalize vectors→ Precision errors in floating-point values can distort 3D reconstructions.

9. References and further reading

  • DICOM Standard, Part 3, Section C.7.6.2 — Image Plane Module.[1]
  • SimpleITK documentation — orientation and DICOM conversion.[2]
  • MONAI documentation — spatial orientation and metadata.[3]
  • pydicom documentation — reading and writing headers and orientation tags.[4]

10. Conclusion

The DICOM format encodes geometry with precision — but that precision only helps if you understand it.

By reading and checking ImagePositionPatient, ImageOrientationPatient, and PixelSpacing, you can diagnose most orientation issues before they ruin your 3D visualization.

In medical imaging, orientation is anatomy — and a single misplaced sign can literally turn the patient upside down.

A glowing magnifying glass highlights “p < 0.05” beside a winding path toward a mountain marked “Large Effect Size,” while a roadside sign warns that statistical significance is not the same as practical importance.

Effect Size

Posted on by Michele Danilo Pierri

Effect Size: What It Is and Why It Matters More Than Statistical Significance

A result can be statistically significant — yet practically meaningless. Learn how effect size reveals the real-world impact of research findings.

Introduction: The Hidden Problem with p-values

You’ve probably seen headlines like:

“New Study Shows Coffee Improves Memory!”

But what if the improvement was just 0.3 points on a 100-point test?

Technically “significant” — but is it meaningful?

This is where effect size comes in.

While p-values tell us whether an effect exists, effect size tells us how large that effect is — a crucial distinction often overlooked in science, education, and media.

In this article, you’ll learn:

  • What effect size really means
  • How to calculate and interpret common measures (like Cohen’s d)
  • Why it’s essential for sound scientific reasoning
  • Best practices for reporting it in research

Let’s go beyond significance testing and focus on what truly matters: practical importance.

What Is Effect Size?

Effect size is a quantitative measure of the magnitude of a phenomenon. Unlike p-values, which depend heavily on sample size, effect size provides a standardized metric that reflects the strength of a relationship or difference — independent of how many people were studied.

In simple terms:

p-value: “Is there an effect?” → answers statistical significance

Effect size: “How big is the effect?” → answers practical significance

For example:

Two teaching methods differ by 5 points in average test scores.

With a small class, the difference might not be significant (high p-value).

With a huge sample, even a 0.5-point difference could be “significant” (low p-value).

But only effect size tells you whether those 5 (or 0.5) points matter in practice.

Effect Size vs. Hypothesis Testing: Key Differences

Featurep-value / Null Hypothesis TestingEffect Size
PurposeTest if an effect is likely due to chanceMeasure the strength of the effect
Depends on sample sizeYes — larger samples increase significanceNo — it’s independent of N
Tells youWhether an effect existsHow large the effect is
Common misuseMistaking statistical significance for importanceIgnoring it altogether

Key insight: A small effect can be highly significant with a large sample — but still too weak to justify policy changes, clinical use, or educational reform.

Side-by-side comparison showing how large samples can detect tiny, irrelevant effects.

How to Calculate Effect Size: Cohen’s d

One of the most widely used measures is Cohen’s d, ideal for comparing the means of two groups.

Formula:

d = \frac{\bar{X}_1 - \bar{X}_2}{s_{\text{pooled}}}

Where:

\bar{X}_1 and \bar{X}_2 are the means of the two groups

s_{\text{pooled}} is the pooled standard deviation:

    \[ s_{\text{pooled}} = \sqrt{\frac{(n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1 + n_2 - 2}} \]

If group sizes are equal, you can approximate s_{\text{pooled}} as the average of the two standard deviations

Example: Teaching Method Experiment

GroupMean ScoreSDN
New Method78.412.130
Traditional72.611.830
  1. Difference in means: 78.4 - 72.6 = 5.8
  2. Pooled SD ≈ \sqrt{\frac{12.1^2 + 11.8^2}{2}} \approx 11.95
  3. Cohen’s d = \frac{5.8}{11.95} \approx 0.48

Interpretation: d ≈ 0.48 → medium effect size

Even without knowing the p-value, we now know the intervention had a moderately strong impact.

Two-group comparison with mean values and confidence intervals

Interpreting Cohen’s d: Rules of Thumb

Jacob Cohen proposed general guidelines for interpreting d:

Cohen’s dInterpretation
0.2Small effect
0.5Medium effect
0.8Large effect

These are benchmarks, not strict rules. Context matters: in education, a d of 0.4 might be very meaningful, in medicine, even d = 0.3 could justify a new treatment if scalable.

Use them as starting points — not final judgments.

A horizontal scale from 0 to 1+ with labeled zones (small/medium/large) and real-world analogies

When Should You Report Effect Size?

Best practices recommend reporting effect size in all empirical studies, especially when:

  • Comparing groups (t-tests, ANOVA)
  • Measuring associations (correlations, regression)
  • Conducting meta-analyses
  • Evaluating interventions (education, psychology, health)

Major journals (APA, APA-style publications) require effect sizes alongside p-values.

Other Common Effect Size Measures

While Cohen’s d is great for mean differences, other contexts require different metrics:

TestEffect SizeRange
t-test (independent)Cohen’s d, Hedges’ g−∞ to +∞
ANOVAEta-squared (η²), Omega-squared (ω²)0 to 1
CorrelationPearson’s r−1 to +1
Chi-squareCramer’s V0 to 1
RegressionR², f²0 to 1

Why Effect Size Matters Researchers

Understanding effect size helps you:

  • Avoid overinterpreting statistically significant but trivial results
  • Compare findings across different studies and scales
  • Design better experiments (via power analysis)
  • Communicate results more honestly and transparently

Power analysis — used to determine required sample size — depends directly on expected effect size.

No effect size? You can’t plan a well-powered study.

Cosmic-scale artwork showing a balance between 'p-value' on one side and 'Effect Size' on the other, symbolizing the need to prioritize meaningful results in science.

Conclusion: Significance ≠ Importance

Let’s summarize the key takeaways:

  1. p-value answers: “Is the effect real?”
  2. Effect size answers: “How big is it?”
  3. A result can be significant but trivial — always check both.
  4. Cohen’s d is a powerful tool
  5. Interpret using benchmarks: 0.2 (small), 0.5 (medium), 0.8 (large) — but consider context.

Statistical significance tells you if you should pay attention. Effect size tells you how much.

Further Reading

Schober, Patrick MD, PhD, MMedStat*; Vetter, Thomas R. MD, MPH†. Effect Size Measures in Clinical Research. Anesthesia & Analgesia 130(4):p 869, April 2020. | DOI: 10.1213/ANE.0000000000004684

Kallogjeri D, Piccirillo JF. A Simple Guide to Effect Size Measures. JAMA Otolaryngol Head Neck Surg. 2023 May 1;149(5):447-451. doi: 10.1001/jamaoto.2023.0159. PMID: 36951858.

Aarts S, van den Akker M, Winkens B. The importance of effect sizes. Eur J Gen Pract. 2014 Mar;20(1):61-4. doi: 10.3109/13814788.2013.818655. Epub 2013 Aug 30. PMID: 23992128.

Paul Monsarrat, Jean-Noel Vergnes, The intriguing evolution of effect sizes in biomedical research over time: smaller but more often statistically significant, GigaScience, Volume 7, Issue 1, January 2018, gix121, https://doi.org/10.1093/gigascience/gix121

Pereira, T. V., Horwitz, R. I., & Ioannidis, J. P. A. (2012). Empirical evaluation of very large treatment effects in randomized controlled trials. JAMA, 308(16), 1689–1696. https://doi.org/10.1001/jama.2012.13444

“The primary product of a research inquiry is one or more measures of effect size, not p-values.” — Jacob Cohen (1994)

In an early-20th-century hospital room, two masked surgeons confer with a formally dressed man beside a reclining patient, surrounded by anatomical sketches, scientific diagrams, glass infusion bottles, and antique medical instruments.

Percutaneous Interventions Expand: Key Insights from 2025 ESC/EACTS Guidelines

Posted on by Michele Danilo Pierri

Introduction

The 2025 ESC/EACTS Guidelines on Valvular Heart Disease represent a pivotal shift in the balance between surgical and percutaneous interventions. While surgical aortic valve replacement (SAVR) remains the gold standard for younger, low-risk patients, the new guidelines expand the role of transcatheter aortic valve implantation (TAVI) and transcatheter edge-to-edge repair (TEER, MitraClip). This evolution reflects both the maturation of percutaneous technologies and mounting evidence supporting their safety and effectiveness across broader patient populations.

TAVI

1. Lower Age Threshold for TAVI

A significant change involves the age threshold. While the 2021 guidelines recommended TAVI primarily for patients ≥75 years, the 2025 update lowers this to

≥70 years (Class I, Level A)

. Surgical replacement is now indicated only for low-risk patients under 70, while patients aged 70–74 require Heart Team evaluation for personalized treatment decisions. This adjustment substantially increases the eligible population for TAVI.

2. Expanding Interventions for Asymptomatic Patients

Previously, the 2021 guidelines required asymptomatic patients with severe high-gradient aortic stenosis to meet multiple “trigger” conditions (such as very high Vmax, elevated BNP levels, or rapid progression) before intervention was justified. The 2025 guidelines now offer a simplified approach:

intervention may be considered (Class IIa, Level A)

in these patients, even without specific triggers, as long as the procedural risk remains low. This change represents a significant shift in approach, recognizing the potential dangers of waiting when patients have severe disease, even without symptoms.

3. Inclusion of Previously Excluded Valve Anatomies

The new guidelines formally approve TAVI for

bicuspid valves (Class IIb, Level B)

in selected high-risk patients with favorable anatomy. Even more significantly, they now permit

TAVI for pure native aortic regurgitation (Class IIb, Level B)

in symptomatic, inoperable patients with suitable anatomy. These conditions were either omitted or discouraged in the 2021 guidelines, representing a clear expansion of percutaneous treatment options.

4. Streamlined Diagnostic and Procedural Workflow

The 2025 document reflects the maturation of TAVI practice. Coronary computed tomography angiography (CCTA) is now accepted as an alternative to invasive angiography when image quality is adequate (Class IIa, Level B). Additionally, the guidelines have narrowed the criteria for percutaneous coronary intervention (PCI), which is now

restricted to lesions ≥90%

in vessels ≥2.5 mm, compared to the previous broader 70% threshold. These changes aim to reduce unnecessary procedures and minimize procedural burden for patients.

5. Antithrombotic Therapy: Simplified Approach

The new guidelines establish that

single antiplatelet therapy (SAPT) is recommended

as the standard approach after TAVI. Conversely, dual antiplatelet therapy (DAPT) is

not recommended

except when there are other specific indications. This straightforward guidance eliminates previous uncertainty and improves patient safety.

Clinical scenario2021 ESC/EACTS2025 ESC/EACTS
Age and risk (tricuspid)TAVI ≥75 y: I A; SAVR <75 y: I A; 70–75 y: Heart TeamTAVI ≥70 y: I A; SAVR <70 y: I A; 70–74 y: Heart Team
Severe symptomatic AS at high riskTAVI: I ATAVI: I A
Severe asymptomatic AS, high gradient, preserved EFIntervention: IIa A with specific triggers (BNP↑, Vmax ≥5.5 m/s, rapid progression)Intervention: IIa A even with just low procedural risk
Bicuspid valveNo formal recommendation / selected casesTAVI: IIb B if increased risk and favorable anatomy
Native aortic regurgitationTAVI not recommendedTAVI: IIb B if symptomatic, inoperable, favorable anatomy
Non-transfemoral accessIIb C (may)IIa B (should)
Coronary work-upICA standardAdequate CCTA can replace ICA (IIa B)
Associated PCIPCI to be considered if stenosis ≥70%PCI to be considered only if stenosis ≥90% (≥2.5 mm)
Post-TAVI antithrombotic therapySAPT preferred; DAPT acceptable in some casesSAPT recommended; DAPT not recommended unless other indications

Critical Appraisal

These changes collectively highlight a significant expansion of percutaneous therapies relative to surgical approaches. While SAVR remains essential for young, low-risk patients and complex anatomical cases, its overall role is diminishing. TAVI has become the predominant option for patients ≥70 years, certain asymptomatic individuals, and even previously excluded conditions such as bicuspid valves and native aortic regurgitation.

This shift represents more than a technological advancement—it marks a conceptual transformation. Percutaneous interventions have evolved from “alternative options” for inoperable cases into mainstream strategies across a broad clinical spectrum. The Heart Team’s role becomes increasingly vital in weighing factors such as long-term prosthesis durability, patient comorbidities, and individual preferences in this rapidly evolving landscape.

TEER

1. Primary Mitral Regurgitation (Degenerative)

In the 2021 guidelines, TEER for PMR was a Class IIb, Level B recommendation—meaning it “may be considered” in symptomatic patients deemed unsuitable for surgery.

2025 Update: TEER has been upgraded to Class IIa, Level B—”should be considered”—for symptomatic patients at high or prohibitive surgical risk with suitable anatomy.

This upgrade is more than just semantic; it establishes TEER as a standard treatment option for the high-risk degenerative population rather than an exceptional measure.

2. Secondary Mitral Regurgitation (Ventricular, “COAPT-like”)

This represents the most significant change in the guidelines. In 2021, TEER carried a Class IIa or IIb indication, varying based on patient selection and optimization of medical therapy.2025 Update: TEER for ventricular secondary MR has been upgraded to Class I, Level A—the strongest possible recommendation, supported by high-quality evidence. It is now recommended for reducing heart failure hospitalizations and improving quality of life in symptomatic patients with suitable anatomy, even after optimal guideline-directed medical therapy (GDMT) and, when indicated, cardiac resynchronization therapy (CRT).

This evolution reflects the enduring impact of the COAPT trial and establishes TEER as a first-line interventional strategy for a substantial subset of heart failure patients.

3. Secondary Mitral Regurgitation (Atrial Form)

The 2025 guidelines introduce a novel differentiation between ventricular SMR (caused by LV dilatation/dysfunction) and atrial SMR (resulting from atrial/annular dilatation with preserved LV function).

  • Surgery (mitral repair plus maze and/or LAA closure) is Class IIa, Level B for symptomatic patients who are surgical candidates.
  • TEER is Class IIb, Level B for symptomatic patients unsuitable for surgery, after rhythm optimization.

This represents an important conceptual refinement, recognizing that atrial SMR may differ from ventricular forms in both pathophysiology and treatment approach.

Critical Appraisal

These changes collectively mark a significant evolution of TEER from a last-resort option to a cornerstone therapy:

  • For PMR, TEER now occupies a well-defined position in the treatment of high-risk patients.
  • In ventricular SMR, its elevation to Class I, Level A standard of care transforms how clinicians manage heart failure patients with persistent symptoms.
  • For atrial SMR, the guidelines introduce a clear structure previously absent, delineating appropriate circumstances for surgical versus TEER approaches.

While surgical mitral repair continues to be the gold standard for younger, low-risk PMR patients, the reach of percutaneous intervention continues to broaden, mirroring the trajectory seen with TAVI in aortic valve treatment.

Clinical Scenario2021 ESC/EACTS2025 ESC/EACTS
Primary MR (degenerative, high/prohibitive surgical risk)IIb B – may be consideredIIa B – should be considered
Secondary MR – ventricular (functional, COAPT-like)IIa/IIb depending on selection and GDMTI A – recommended
Secondary MR – atrial formNot specifically differentiatedSurgery (repair+maze/LAAO): IIa B; TEER: IIb B if not surgical candidate, after rhythm optimization
Surgical mitral repair (PMR, low-risk)Gold standard, Class IGold standard, Class I

Critical Considerations

Durability. Long-term outcomes remain the primary limitation of percutaneous interventions. For TAVI, 5-year results in low-risk trials (e.g., PARTNER 3, Evolut Low Risk) demonstrate noninferiority to SAVR, but uncertainty exists beyond 8–10 years, particularly in younger patients. For TEER, COAPT shows sustained benefit at 5 years, though patient attrition reflects the advanced heart failure population.

Learning curve and volume. Outcomes for both TAVI and TEER strongly correlate with operator and center experience. High-volume centers achieve better procedural success rates and fewer adverse events, supporting the guidelines’ emphasis on Heart Valve Centers.

Patient selection. Appropriate patient selection remains essential. Lifetime management considerations, anatomic suitability, comorbidities, and future treatment options should guide clinical decisions. The 2025 guidelines emphasize the Heart Team’s central role in balancing surgical durability with percutaneous accessibility.

Conclusion

The 2025 ESC/EACTS Guidelines represent a decisive shift in valvular heart disease management, with TAVI and TEER evolving from selective alternatives to mainstream treatments. For aortic stenosis, TAVI’s scope has expanded beyond its traditional boundaries through a lower age threshold, extended indications for asymptomatic patients, and careful inclusion of bicuspid valves and native regurgitation cases. In mitral disease, TEER has progressed from a last-resort option to a standard of care, earning Class I, Level A endorsement for ventricular secondary MR and gaining prominence in high-risk primary MR.

Surgical approaches remain essential—especially for young, low-risk patients and complex anatomical cases—but the balance is clearly shifting toward less invasive interventions. The Heart Team now faces a new clinical landscape where transcatheter therapies serve as core components of evidence-based practice rather than mere alternatives. This transformation stems from technological advancement, stronger clinical evidence, and increasing emphasis on personalized, minimally invasive care.

References

Praz, F., Borger, M. A., Lanz, J., Marin-Cuartas, M., Abreu, A., Adamo, M., … Zamorano, J. L. (2025). 2025 ESC/EACTS Guidelines for the management of valvular heart disease. European Heart Journal. Advance online publication.

https://doi.org/10.1093/eurheartj/ehaf194

Vahanian, A., Beyersdorf, F., Praz, F., Milojevic, M., Baldus, S., Bauersachs, J., … Zamorano, J. L. (2021). 2021 ESC/EACTS Guidelines for the management of valvular heart disease. European Heart Journal, 43(7), 561–632.

https://doi.org/10.1093/eurheartj/ehab395

Mack, M. J., Thourani, V. H., Pibarot, P., Hahn, R. T., Genereux, P., Kodali, S. K., … Leon, M. B. (2023). Transcatheter aortic-valve replacement in low-risk patients at five years. New England Journal of Medicine, 389(21), 1949–1960.

https://doi.org/10.1056/NEJMoa2307447

Stone, G. W., Lindenfeld, J., Abraham, W. T., Kar, S., Lim, D. S., Mishell, J. M., … Mack, M. J.; COAPT Investigators. (2018). Transcatheter mitral-valve repair in patients with heart failure. New England Journal of Medicine, 379(24), 2307–2318.

https://doi.org/10.1056/NEJMoa1806640

Carroll, J. D., Vemulapalli, S., Dai, D., Matsouaka, R., Blackstone, E., Edwards, F., … Grover, F. (2017). Procedural experience for transcatheter aortic valve replacement and relation to outcomes: The STS/ACC TVT Registry. Journal of the American College of Cardiology, 70(1), 29–41.

https://doi.org/10.1016/j.jacc.2017.04.056

Sorajja, P., Vemulapalli, S., Feldman, T., Mack, M., Holmes, D. R., Stebbins, A., … Ailawadi, G. (2017). Outcomes with transcatheter mitral valve repair in the United States: An STS/ACC TVT Registry report. Journal of the American College of Cardiology, 70(19), 2315–2327.

https://doi.org/10.1016/j.jacc.2017.09.015

An elderly physician in a white coat carefully shapes a red calibration curve on a large graph comparing predicted and observed risk, in a warmly lit early 20th-century medical study rendered in muted sepia and ochre tones.

Calibration of Predictive Risk Models: A Guide for Clinicians

Posted on by Michele Danilo Pierri

Introduction: Understanding Calibration Challenges

Consider a thermometer that perfectly identifies when one temperature is higher or lower than another, but consistently reads 5 degrees too high. This thermometer has good discrimination (it correctly ranks temperatures), but poor calibration (its absolute values are inaccurate).

The same problem occurs with clinical risk models such as EuroSCORE: they can be excellent at distinguishing high-risk patients from low-risk ones, but the specific probabilities they provide might not accurately reflect the reality in your population.

Understanding Calibration

The calibration of a predictive model measures how accurately the predicted probabilities match the actual frequency of observed events.

In practical terms:

For optimal calibration, when a model predicts a 20% risk for a group of patients, exactly 20% of those patients should experience the event.

The model overestimates the risk when the event occurs in only 15% of cases.

Conversely, risk is underestimated when the event occurs in 25% of cases.

Why is Calibration Important?

Poor calibration can lead to significant clinical consequences when risk prediction models fail to accurately reflect real-world outcomes. These discrepancies between predicted and actual risk can substantially impact patient care decisions:

Overestimation of risk → Unnecessarily aggressive treatments and interventions that patients don’t need, potentially exposing them to adverse effects, complications, and increased healthcare costs without proportional clinical benefit

Underestimation of risk → Failure to provide appropriate treatments and inadequate attention to truly high-risk patients, potentially resulting in missed opportunities for preventive interventions, delayed recognition of deterioration, and suboptimal management strategies

Calibration Assessment Methods

Graphical methods

Calibration Plot

The most intuitive method to visualize calibration:

Patients are grouped according to their predicted risk levels (0-10%, 10-20%, etc.)

The actual proportion of events is calculated for each group

A graph displays predicted risk (X-axis) against observed risk (Y-axis)

Perfect calibration: all points align with the diagonal line

Points below the line: risk overestimation

Points above the line: risk underestimation

Calibration plot in case of overestimation or underestimation

Calibration Belt

An enhanced version of the calibration plot that incorporates confidence bands, helping clinicians distinguish between statistically significant deviations and random variations in the data.

Calibration belt for overestimation and underestimation

Classic Quantitative Methods

Calibration In The Large (CITL)

This represents the simplest calibration assessment method:

Compares the sum of all predicted events versus observed events

Ratio = Observed events / Predicted events

Value of 1: perfect calibration

Value < 1: global overestimationValue > 1: global underestimation

Hosmer-Lemeshow Test

The classic statistical test for calibration assessment, but use with caution as it has several important limitations:

Divides patients into predetermined risk groups (usually 10 equal-sized groups or deciles of risk)

Systematically compares predicted versus observed event rates within each of these risk groups

Produces a chi-square statistic with a corresponding P-value where values > 0.05 traditionally suggest acceptable calibration (failure to reject the null hypothesis of good fit)

Despite its widespread use in the literature, it has several critical limitations that significantly reduce its utility in modern predictive modeling:

Results depend arbitrarily on the number of groups chosen by the researcher, with different grouping strategies potentially yielding contradictory conclusions from the same dataset

Demonstrates inadequate sensitivity when applied to smaller sample sizes, potentially failing to detect meaningful calibration issues

Paradoxically becomes overly sensitive with very large samples, where clinically insignificant deviations may be flagged as statistically significant problems

Provides only a global assessment without indicating specific regions of the risk spectrum where calibration problems are most pronounced

Produces notably unstable results across different underlying risk distributions, limiting comparability between populations

Advanced Quantitative Methods

Integrated Calibration Index (ICI)

Measures the area between the observed calibration curve and the ideal one:

Values close to 0 indicate excellent calibration

Typically, an ICI < 0.01 is considered good

Provides a comprehensive global measure of calibration error

E50 and E90 (Error Percentiles)

E50 (median absolute error)**:

50% of patients have calibration errors lower than this value, providing a robust central tendency measure that is not influenced by extreme outliers

Indicates the “typical” quality of calibration in your patient population and serves as a reliable benchmark for routine clinical application

E90 (90th percentile of absolute error)**:

Only 10% of patients have calibration errors exceeding this threshold, making it a valuable indicator of the model’s worst-case performance

Indicates how severe the errors can be in the worst cases, which is particularly important when making high-stakes clinical decisions where safety margins are critical

Brier Score

Measures the average “distance” between predictions and actual outcomes:

Formula: Average of (Prediction – Reality)²

Range: 0 (perfect calibration) to 0.25 (random prediction)

Provides a comprehensive measure that combines both calibration and discrimination

Alternative Statistical Tests

Spiegelhalter’s Z-test

Provides a comprehensive single Z-score for overall calibration assessment, making it easier to interpret compared to multiple metrics

Statistical interpretation follows standard normal distribution where |Z| < 1.96 indicates acceptable calibration at the 95% confidence level

More robust than Hosmer-Lemeshow for large samples, as it doesn’t suffer from the same sensitivity issues when sample sizes increase

Particularly useful for comparing calibration across different predictive models applied to the same population

Maintains consistent performance regardless of the underlying risk distribution in your patient cohort

Le Cessie-Van Houwelingen Test

Does not require arbitrary division into groups, making it more consistent and reliable than methods that depend on subjective grouping choices

Particularly suitable for logistic regression models and provides a more sophisticated assessment of the goodness-of-fit for probability predictions in clinical settings

Less influenced by data distribution anomalies, offering more robust performance across diverse patient populations with varying risk profiles

Evaluates the overall calibration quality by examining the squared differences between observed outcomes and predicted probabilities

Maintains statistical power even with smaller sample sizes, making it valuable for specialized clinical applications with limited available data

Summary Table: Calibration Assessment Methods

MethodTypeInterpretationAdvantagesLimitations
Calibration PlotGraphPoints on diagonal = perfectVisual, intuitiveSubjective, depends on grouping
Calibration BeltGraphDiagonal + confidence bandsShows statistical significanceMore complex interpretation
CITLNumeric1 = perfect, <1 = overestimation, >1 = underestimationSimple single valueOnly overall assessment
Hosmer-Lemeshow testTestp>0.05 = accettable calibrationClassic, widely knownUnstable, poorly informative
ICINumeric0 = perfect, <0.01 = goodRobust global measureDoesn’t pinpoint issues
E50/E90NumericMedian/90th percentile of errorsClinically interpretableE90 sensitive to outliers
Brier ScoreNumeric0 = perfect, <0.25 = usefulComprehensive measureCombines multiple aspects
Spiegelhalter Z-testTest|Z| < 1.96 = acceptableSingle standardized scoreLess intuitive clinically
Le Cessie testTestp > 0.05 = acceptableNo arbitrary groupingMore complex than H-L

Correcting Calibration Errors

Platt Scaling

Platt Scaling represents the most extensively adopted and widely implemented methodology for correcting calibration issues in predictive models across various clinical domains:

How it works:

Applies a specialized form of logistic regression transformation directly to the original prediction outputs, effectively recalibrating them without altering their fundamental ranking properties

Formula: P_calibrated = 1 / (1 + exp(A × original_score + B)), where the equation creates a sigmoid-shaped adjustment that can correct both over-prediction and under-prediction issues across the risk spectrum

Parameters A and B are statistically derived using your local patient population data through maximum likelihood estimation, ensuring the calibration correction is specifically tailored to your clinical context

Advantages:

Simple to implement in most statistical packages and clinical decision support systems without requiring complex computational resources

Maintains the critical rank-order of predictions, ensuring the discrimination ability of the model (its capacity to separate high-risk from low-risk patients) remains completely unchanged

Effective for addressing most common calibration biases encountered in clinical predictive models, including systematic over-prediction and under-prediction patterns

Disadvantages:

Requires a separate, independent dataset for calibration parameter estimation and subsequent validation to avoid overfitting, which may be challenging in resource-limited settings

Assumes a sigmoidal relationship between the original scores and observed outcomes, which may not adequately correct more complex non-monotonic calibration issues that occasionally arise in heterogeneous patient populations

Temperature Scaling

Temperature Scaling represents a streamlined and computationally efficient variant of the more complex Platt Scaling methodology, providing a more accessible approach to model recalibration while maintaining essential corrective capabilities:

Employs a single correction parameter (temperature) that functions as a scaling factor applied uniformly across all predictions, significantly reducing computational complexity while still addressing systematic calibration issues

Offers enhanced implementation simplicity and reduced computational requirements, making it particularly suitable for resource-constrained clinical environments, though it provides less flexibility for correcting complex non-linear calibration patterns compared to multi-parameter approaches

Isotonic Regression

Does not assume a specific form of the relationship between predicted probabilities and actual outcomes, allowing for more flexible correction of complex calibration errors across different risk ranges

Particularly useful when the bias is very irregular or non-monotonic, such as when a model simultaneously overestimates risk in some patients while underestimating it in others, depending on where they fall in the risk spectrum

More complex to implement than simpler methods like Platt scaling, requiring specialized algorithms and additional computational resources, though the improved calibration accuracy often justifies this increased complexity in high-stakes clinical applications

Practical Advice

Begin calibration analysis with “raw” data to evaluate a model’s actual clinical performance in your specific patient population. This initial assessment often provides sufficient insight into the model’s direct applicability to your clinical practice.

If needed, proceed with calibration. This process isolates the model’s intrinsic predictive ability and enables fairer comparisons between models developed across different populations.

Follow this sequence: conduct raw data analysis first, then calibrate if necessary, and finally evaluate post-calibration performance.

You can validate calibration using several approaches:

Temporal split: Use historical data for calibration and recent data for validation

Cross-validation: Divide your dataset into multiple parts for more robust validation

Geographic split: Calibrate using data from one center and validate with data from others

Goal-oriented approach

For Daily Clinical Use

Incorporate calibration plots alongside ICI and E50/E90 metrics as your primary assessment tools to efficiently evaluate model performance in routine clinical practice

These methodologies provide an optimal balance between statistical rigor and intuitive interpretation, allowing healthcare providers to quickly grasp calibration quality without requiring extensive statistical expertise

The visual nature of calibration plots coupled with the numerical precision of ICI and E50/E90 metrics creates a comprehensive assessment framework that can be readily communicated to clinical teams during decision-making processes

For Scientific Research

Implement a multi-dimensional approach by strategically combining complementary graphical representations and quantitative metrics to capture the full spectrum of calibration characteristics

Establish methodological robustness by utilizing multiple independent indicators that collectively provide a comprehensive assessment of model calibration across different dimensions of performance

Exercise caution regarding over-reliance on the Hosmer-Lemeshow test due to its known limitations with large sample sizes and sensitivity to arbitrary grouping decisions

Maintain a dual focus on both statistical significance and clinical relevance by contextualizing calibration findings within the specific medical domain and intended application, recognizing that statistically significant deviations may not always translate to clinically meaningful differences in patient outcomes

For Developing New Models

Integrate calibration assessment into the earliest phases of model development rather than treating it as an afterthought, establishing it as a core design consideration alongside discrimination metrics

Implement rigorous methodological standards by systematically setting aside truly independent validation datasets that remain completely untouched during model development and initial calibration phases

Ensure comprehensive documentation of all calibration methodologies, decision thresholds, and statistical approaches used throughout the development process to enhance transparency and facilitate proper implementation by other researchers and clinicians

Consider the temporal stability of calibration by designing periodic reassessment protocols that can identify calibration drift as patient populations and clinical practices evolve over time

Best Practices

Always use multiple methods – combine graphical, quantitative, and (when appropriate) statistical test approaches

Avoid using Hosmer-Lemeshow test – opt for more robust alternatives

Look beyond statistical significance – prioritize clinical relevance

Reserve independent data for validation when applying calibration corrections

Document thoroughly all methods used and defined acceptability thresholds

1. ALWAYS start with: Calibration plot + CITL
↓ 
2. If CITL ≠ 1 → Global calibration problem 
↓ 
3. Add: ICI + E50 + E90 to quantify 
↓ 
4. If statistical test needed: - Small sample (<500): Le Cessie test - Large sample (>500):      Spiegelhalter test - NEVER use Hosmer-Lemeshow alone 
↓
5. If correction needed: - Simple: Platt Scaling - Complex: Isotonic Regression 
↓ 
6. ALWAYS validate on independent data

​

​

References and Suggested Literature

Hosmer-Lemeshow test

Hosmer, D.W., Lemeshow, S. (1980). A goodness-of-fit tests for the multiple logistic regression model. Communications in Statistics, 10, 1043-1069[1]

Kramer, A.A., Zimmerman, J.E. (2007). Assessing the calibration of mortality benchmarks in critical care: The Hosmer-Lemeshow test revisited. Critical Care Medicine, 35(9), 2052-2056[2]

Integrated Calibration Index (ICI)

Austin, P.C., Steyerberg, E.W. (2019). The Integrated Calibration Index (ICI) and related metrics for quantifying the calibration of logistic regression models. Statistics in Medicine, 38(21), 4051-4065[3]

Calibration Plot

Austin, P.C., Steyerberg, E.W. (2014). Graphical assessment of internal and external calibration of logistic regression models by using loess smoothers. Statistics in Medicine, 33(3), 517-535[4]

Calibration Belt

Finazzi, S., Poole, D., Luciani, D., Cogo, P.E., Bertolini, G. (2011). Calibration belt for quality-of-care assessment based on dichotomous outcomes. PLOS ONE, 6(2), e16110[6]

Brier Score

Brier, G.W. (1950). Verification of forecasts expressed in terms of probability. Monthly Weather Review, 78(1), 1-3

Le Cessie-van Houwelingen test

le Cessie, S., van Houwelingen, J.C. (1991). A goodness-of-fit test for binary regression models, based on smoothing methods. Biometrics, 47(4), 1267-1282[5]

Platt Scaling

Platt, J. (1999). Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Advances in Large Margin Classifiers, 61-74

Review

Van Calster, B., McLernon, D.J., van Smeden, M., Wynants, L., Steyerberg, E.W. (2019). Calibration: the Achilles heel of predictive analytics. BMC Medicine, 17, 230[7]

Huang, Y., Li, W., Macheret, F., Gabriel, R.A., Ohno-Machado, L. (2020). A tutorial on calibration measurements and calibration models for clinical prediction models. Journal of Biomedical and Health Informatics, 24(4), 1079-1090[8]

Conclusions

Calibration isn’t merely a technical consideration but a fundamental requirement for the safe and effective use of predictive models in medicine. A well-calibrated model accurately identifies not just who faces risk but precisely how much risk they face, enabling clinicians to make more informed and appropriate decisions.