micheledpierri.com

  • HOME
    • Python
    • Statistics
    • Data Analysis
    • Machine Learning
  • WRITINGS
  • VISIONS
  • ABOUT
Home / Statistics / Propensity Score
a young gardener prunes a row of uneven trees, shaping them into uniform ones

Propensity Score

Propensity Score in Medical Research

A Practical Introduction with Python

Last update: June 2026

Author: Michele D. Pierri

Reading time: 20 minutes

1. Why do we need propensity scores?

In clinical research, randomized controlled trials (RCTs) represent the methodological reference point for estimating causal effects, because randomization (at least in principle) distributes both measured and unmeasured confounders across treatment groups. This is the core strength of the RCT design, and it is why observational data requires a fundamentally different approach.

In observational studies (retrospective cohorts, registries, administrative databases), treatment assignment is not randomized. Sicker patients tend to receive one intervention over another; surgeons choose techniques based on age, frailty, and comorbidities; hospital protocols vary in ways that are rarely fully documented. All of this means the groups being compared are not exchangeable at baseline.

A naïve comparison of raw outcomes between treatment A and treatment B ends up confounding the effect of treatment with the effect of different baseline risks. This is confounding by indication, and it is one of the most persistent problems in observational clinical research.

The propensity score is a tool designed to address this problem, at least partially. That last qualifier matters.


2. What is a propensity score?

Formally, the propensity score (PS) is:

The probability of receiving the treatment, given the observed covariates.

If we call:

  • T = 1 for treated patients, T = 0 for controls
  • X the vector of baseline covariates (age, sex, comorbidities, risk scores…)

then the propensity score is:

e(X) = P(T = 1 \mid X)

In practice, e(X) is usually estimated with a logistic regression (or alternative models: random forest, gradient boosting) where the dependent variable is the treatment assignment, not the clinical outcome. This distinction sometimes trips people up.


3. What is it used for?

Once we have an estimated propensity score for each patient, it can be used to create more comparable groups in several ways:

  • Matching: pair treated and control patients with similar propensity scores.
  • Stratification (subclassification): divide patients into PS strata (e.g., quintiles) and compare outcomes within each stratum.
  • Weighting: use the inverse probability of treatment (IPTW) to re-weight the sample and generate a pseudo-population where treatment is independent of covariates.
  • Covariate adjustment: include the propensity score (or its logit) directly as a covariate in an outcome regression model.
  • Doubly robust methods: combine PS weighting with an outcome regression model.

These approaches differ in how they balance the trade-off between sample size, precision, and covariate balance. The choice depends on the study question, the available sample, and what you are trying to estimate.


4. When should we use propensity scores?

The clearest use cases involve retrospective cohort studies built on hospital records or registries, comparative effectiveness research (e.g., surgery vs. PCI, valve A vs. valve B), and any situation where randomization is impossible or unethical but some approximation of causal inference is still needed.

Propensity scores tend to work well when the number of covariates is large relative to the number of events, when there is a clear treatment-versus-control dichotomy, and when the goal is to emulate a hypothetical randomized trial using observational data. That said, they are not universally appropriate.

They are less useful, and can be misleading, when unmeasured confounding is severe and the critical variables are absent from the dataset. In cardiac surgery, for example, frailty and functional status are rarely captured in administrative records, yet they often drive both treatment selection and outcomes. No PS model can correct for what was never measured. Poor overlap between groups is another limiting condition: if very young patients consistently receive treatment A and very old patients consistently receive treatment B, there is no adequate basis for comparison regardless of the analytical method applied.


5. Ways to use the propensity score

5.1 Matching

Propensity score matching (PSM) creates pairs (or sets) of treated and control patients with similar PS values.

  • 1:1 nearest-neighbor matching: each treated patient is matched to the closest control (in PS space).
  • Optional caliper: only matches within a defined PS distance are accepted (e.g., 0.2 SD of the logit PS).
  • Can be done with or without replacement.

Matching is intuitive: you are literally comparing similar patients. The matched cohort resembles an RCT design, which makes results easier to communicate to a clinical audience. The main cost is the loss of unmatched patients, which reduces precision and can introduce selection effects depending on who is discarded. Diagnostics matter here. Caliper choice and balance checking are not optional steps.


5.2 Stratification (subclassification)

Here you:

  1. Divide the sample into PS strata (e.g., quintiles).
  2. Compare treated vs. control outcomes within each stratum.
  3. Combine the stratum-specific estimates (e.g., weighted by stratum size).

This approach is relatively simple to implement and transparent in presentation. It is also useful for exploring effect heterogeneity across risk strata. On the other hand, it is less precise than matching or weighting, and residual imbalance within strata is possible, particularly in smaller datasets.


5.3 IPTW (Inverse Probability of Treatment Weighting)

IPTW uses the inverse of the propensity score as a weight:

  • For treated patients (T = 1):

w_i = \frac{1}{e(X_i)}

  • For controls (T = 0):

w_i = \frac{1}{1 - e(X_i)}

In the re-weighted pseudo-population, covariates should be independent of treatment assignment, mimicking what randomization would have achieved. In practice, stabilized weights reduce variance:

w_i = \begin{cases} \frac{P(T=1)}{e(X_i)} & \text{if } T=1 \\ \frac{P(T=0)}{1-e(X_i)} & \text{if } T=0 \end{cases}

The main advantage of IPTW is that it uses all available patients, without discarding anyone. It integrates well with outcome models (weighted regression, Cox models) and is widely reported in the literature. The main vulnerability is sensitivity to extreme PS values: a PS near 0 or 1 generates very large weights, which can inflate variance substantially and produce unstable estimates. This requires careful diagnostics.


5.4 Covariate adjustment and doubly robust methods

  • Covariate adjustment: Include PS directly as a covariate in a regression (e.g., logistic regression of outcome on treatment + PS, or on treatment + logit(PS)).
  • Doubly robust methods: Combine IPTW with an outcome model. If either the PS model or the outcome model is correctly specified (not necessarily both), the estimator remains consistent. This built-in protection against misspecification is the key appeal of doubly robust approaches.

6. Advantages and limitations of propensity scores

Advantages

  • Make assumptions explicit: you must decide which covariates to include, forcing a deliberate analytical choice.
  • Can substantially improve covariate balance compared to naïve regression adjustment alone.
  • Particularly helpful when events are rare, covariates are many, or when you want a design that is more legible to clinical readers.

Limitations

  • Only adjusts for measured confounders. Missing a key prognostic variable from the dataset means the PS cannot account for it. This is not a limitation of the method per se; it is a limitation of the data.
  • Relies on reasonable model specification (or at least a plausible approximation).
  • Requires careful diagnostics: overlap of PS distributions, covariate balance after matching or weighting, and stability of weights.
  • Can be misused as a black box that lends false rigor to an analysis that was not designed with causal inference in mind.

7. A step-by-step Python example with a simulated clinical dataset

In this section, we will:

  1. Simulate a simple clinical dataset (cardiac surgery patients).
  2. Show the biased naïve comparison between treatments.
  3. Estimate propensity scores using logistic regression.
  4. Apply IPTW and compare the estimated treatment effect.

The code below is an educational template, not a production-ready pipeline. The goal is transparency, not brevity.

7.1 Setup

import numpy as np
import pandas as pd

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
import statsmodels.api as sm

If you prefer to avoid statsmodels, you can do everything with sklearn and sample_weight, but statsmodels makes it easier to extract standard errors and odds ratios directly.


7.2 Simulate a clinical dataset

We simulate:

  • N = 2000 patients
  • Baseline covariates:
    • age (years)
    • sex (0 = female, 1 = male)
    • charlson (comorbidity index 0–6)
    • euroscore (simplified risk score 1–10)
  • Treatment: T = 1 (new surgical technique) vs T = 0 (standard technique)
  • Outcome: in-hospital mortality (0/1)

We introduce confounding: high-risk patients are more likely to receive the new technique and also more likely to die, regardless of which technique they receive.

np.random.seed(42)
N = 2000

# Baseline covariates
age = np.random.normal(loc=70, scale=8, size=N)             # years
sex = np.random.binomial(1, 0.6, size=N)                    # 1 = male
charlson = np.random.poisson(lam=2, size=N)                 # comorbidity index
charlson = np.clip(charlson, 0, 6)
euroscore = np.random.uniform(1, 10, size=N)                # simplified risk

df = pd.DataFrame({
    "age": age,
    "sex": sex,
    "charlson": charlson,
    "euroscore": euroscore
})

Now we generate the probability of receiving the new technique as a function of risk factors:

# True treatment assignment model (logistic in the background)
logit_treatment = (
    -5
    + 0.05 * (age - 70)       # older patients slightly more likely to get new technique
    + 0.4 * sex               # males more likely
    + 0.3 * charlson          # higher comorbidity
    + 0.25 * euroscore        # higher risk score
)
p_treatment = 1 / (1 + np.exp(-logit_treatment))

treatment = np.random.binomial(1, p_treatment, size=N)
df["treatment"] = treatment

Now we simulate the true outcome model. The new technique has a mildly beneficial effect (lower mortality), but this will not be apparent in the unadjusted analysis:

# True outcome model
logit_outcome = (
    -6
    + 0.06 * (age - 70)
    + 0.5 * sex
    + 0.4 * charlson
    + 0.35 * euroscore
    - 0.5 * treatment      # new technique reduces log-odds of death
)
p_death = 1 / (1 + np.exp(-logit_outcome))

death = np.random.binomial(1, p_death, size=N)
df["death"] = death

Check crude outcome rates:

df.groupby("treatment")["death"].mean()

You will typically see something like treatment = 0 at 7–8% mortality and treatment = 1 at 9–10%, even though by construction the new technique is protective. This is confounding by indication in action.


7.3 Naïve (biased) comparison

import statsmodels.formula.api as smf

model_naive = smf.logit("death ~ treatment", data=df).fit(disp=False)
print(model_naive.summary())
np.exp(model_naive.params)  # odds ratios

The naïve model typically yields an odds ratio > 1 for the new technique, suggesting it is harmful. We know from the simulation that this is wrong. It is purely an artifact of confounding.


7.4 Estimating the propensity score

X = df[["age", "sex", "charlson", "euroscore"]]
y = df["treatment"]

ps_model = LogisticRegression(
    solver="lbfgs",
    max_iter=1000
)
ps_model.fit(X, y)

# Propensity scores
ps = ps_model.predict_proba(X)[:, 1]
df["ps"] = ps

print("AUC of PS model:", roc_auc_score(y, ps))

A high AUC here means that baseline covariates strongly predict treatment assignment, which is direct evidence that confounding is present.


7.5 Checking overlap

In practice, you must inspect the PS distributions in treated vs. controls (histograms or KDE plots). The basic requirement is overlap: in each region of PS, you should have both treated and control patients. When one group has PS values clustered near 0 and the other near 1, the positivity assumption is violated and causal inference becomes unreliable.

import matplotlib.pyplot as plt

plt.hist(df.loc[df["treatment"] == 1, "ps"], bins=30, alpha=0.5, label="Treated")
plt.hist(df.loc[df["treatment"] == 0, "ps"], bins=30, alpha=0.5, label="Control")
plt.xlabel("Propensity score")
plt.ylabel("Count")
plt.legend()
plt.show()

7.6 Building IPTW weights

p_treat_overall = df["treatment"].mean()
p_control_overall = 1 - p_treat_overall

df["weight"] = np.where(
    df["treatment"] == 1,
    p_treat_overall / df["ps"],
    p_control_overall / (1 - df["ps"])
)
df["weight"].describe()

Check that weights are not extremely large. In real data, truncating at the 1st–99th percentiles before fitting the outcome model is often warranted.


7.7 Weighted outcome model (IPTW)

model_iptw = smf.glm(
    "death ~ treatment",
    data=df,
    family=sm.families.Binomial(),
    freq_weights=df["weight"]
).fit()

print(model_iptw.summary())
np.exp(model_iptw.params)  # odds ratios

7.7.1 Uncertainty estimation: robust SEs and 95% confidence intervals

IPTW changes the variance structure of the data, so standard errors from unweighted models are not appropriate. Use robust (sandwich) SEs:

# Robust (sandwich) SEs and 95% CI for IPTW GLM
model_iptw = smf.glm(
    "death ~ treatment",
    data=df,
    family=sm.families.Binomial(),
    freq_weights=df["weight"]
).fit(cov_type="HC3")   # HC3 is a common robust choice

import numpy as np

coef = model_iptw.params["treatment"]
se   = model_iptw.bse["treatment"]
OR   = np.exp(coef)
CI_l = np.exp(coef - 1.96 * se)
CI_u = np.exp(coef + 1.96 * se)
print(f"IPTW OR={OR:.3f} (95% CI {CI_l:.3f}–{CI_u:.3f})")

If the data have a clustered structure (e.g., patients nested within hospitals), replace HC3 with clustered SEs: cov_type="cluster", cov_kwds={"groups": df["hospital_id"]}.

7.7.2 Trimming extreme weights and sensitivity of CI

# Trim weights at the 1st–99th percentiles
lo, hi = np.percentile(df["weight"], [1, 99])
df["w_trim"] = df["weight"].clip(lo, hi)

model_iptw_trim = smf.glm(
    "death ~ treatment",
    data=df,
    family=sm.families.Binomial(),
    freq_weights=df["w_trim"]
).fit(cov_type="HC3")

coef = model_iptw_trim.params["treatment"]
se   = model_iptw_trim.bse["treatment"]
OR_t = np.exp(coef)
CI_l = np.exp(coef - 1.96 * se)
CI_u = np.exp(coef + 1.96 * se)
print(f"IPTW (trimmed) OR={OR_t:.3f} (95% CI {CI_l:.3f}–{CI_u:.3f})")

Report both sets of results. If trimming barely changes the point estimate or CI, the result is more stable.

7.7.3 Naïve vs IPTW vs PS-adjusted: a compact comparison

# Naïve (unweighted) logistic regression with robust SEs
model_naive = smf.glm(
    "death ~ treatment",
    data=df,
    family=sm.families.Binomial()
).fit(cov_type="HC3")

# PS-adjusted regression (include PS as covariate)
model_ps_adj = smf.glm(
    "death ~ treatment + ps",
    data=df,
    family=sm.families.Binomial()
).fit(cov_type="HC3")

def or_ci(m):
    b  = m.params["treatment"]
    se = m.bse["treatment"]
    OR = np.exp(b)
    lo = np.exp(b - 1.96*se)
    hi = np.exp(b + 1.96*se)
    return OR, lo, hi

rows = []
for name, m in [
    ("Naive", model_naive),
    ("IPTW", model_iptw),
    ("IPTW trimmed", model_iptw_trim),
    ("PS-adjusted", model_ps_adj),
]:
    OR, lo, hi = or_ci(m)
    rows.append({"Model": name, "OR": OR, "CI": f"{lo:.3f}–{hi:.3f}"})

pd.DataFrame(rows)

How to read this table: if the naïve OR is > 1 but IPTW and PS-adjusted show OR < 1 with reasonable CIs, this indicates the original comparison was biased by confounding. The CI width reflects residual uncertainty. Stability across trimmed and untrimmed weights is a useful internal consistency check.

7.7.4 Reporting checklist for uncertainty

  • Specify the covariance type used (e.g., HC3, cluster-robust) and the rationale.
  • Report the exact weight formulation (stabilized ATE) and any trimming rule applied.
  • Provide OR with 95% CI for each method you compare.
  • Consider a forest plot to visualize point estimates and CIs side by side.

7.8 Optional: PS as covariate or matching

Adjusting for PS in a regression

model_ps_adjusted = smf.logit(
    "death ~ treatment + ps",
    data=df
).fit(disp=False)

print(model_ps_adjusted.summary())
np.exp(model_ps_adjusted.params)

This is simpler than IPTW, though results may differ. In some settings, the two approaches yield consistent estimates; in others, the differences can themselves be informative.

Basic 1:1 nearest-neighbor matching

A full matching implementation requires more code, but the logic is:

  1. Split treated and controls.
  2. For each treated patient, find the control with the closest PS.
  3. Keep only matched pairs and analyze the outcome in the matched sample (paired analysis or logistic regression with robust SEs).

8. How to translate this to real clinical data

In real projects:

  1. Define the question clearly. For example: “In adult patients undergoing isolated CABG, is off-pump surgery associated with lower in-hospital mortality compared to on-pump surgery?”
  2. Specify the target estimand. Average Treatment Effect (ATE) vs. Average Treatment Effect on the Treated (ATT): this choice is not trivial and influences which method is most appropriate.
  3. Select covariates. Include all variables that are plausible confounders (affecting both treatment choice and outcome). Do not include mediators or post-treatment variables in the PS model.
  4. Fit the PS model. Start with logistic regression; consider more flexible models (splines, tree-based methods) if the relationship between covariates and treatment assignment is likely non-linear.
  5. Check overlap. Inspect PS distributions. Poor overlap is a reason to restrict the analysis, not to proceed regardless.
  6. Assess covariate balance after PS adjustment. Standardized mean differences (SMD) before and after adjustment, with a declared threshold (typically SMD < 0.1). Failure to achieve balance is not a reason to skip the report; it is a reason to re-specify the PS model.
  7. Estimate the treatment effect. Use the outcome model appropriate to the endpoint (logistic, Cox, linear), with robust SEs if weighting was applied.
  8. Sensitivity analyses. Different PS specifications, alternative estimands, trimming decisions, and a frank discussion of the potential impact of unmeasured confounding.

9. Conclusion

Propensity score methods provide a structured approach to confounding in observational medical studies. When applied and reported carefully, they make analytic assumptions visible, improve group comparability, and help readers evaluate the credibility of a comparison. But they are not a panacea, and the literature contains too many examples of PS analyses used to create an appearance of rigor without the substance behind it.

The Python template above is a starting point. Replace the simulated covariates with your real variables, build and validate the PS model with appropriate diagnostics, and invest time in the balance assessment before moving to the outcome model. The reporting checklist in Section 10 can serve as a minimum standard for what to include in a manuscript.


10. Reporting Checklist: What a Good Propensity Score Study Must Show

This compact checklist can be included in any manuscript, preprint or blog article. It helps readers quickly verify whether the analysis was conducted rigorously.

1. Study Design & Question

  • Clear research question (“Does treatment X reduce outcome Y in population Z?”).
  • Target estimand stated (ATE, ATT, ATC).
  • Inclusion/exclusion criteria defined.

2. Covariate Selection

  • Covariates chosen based on clinical rationale + literature review.
  • Only pre-treatment variables included.
  • No mediators or post-treatment variables used in the PS model.
  • Variables influencing treatment and outcome included.

3. Propensity Score Model

  • Type of model specified (logistic regression, boosting, random forest, etc).
  • Full model equation or variables listed.
  • Handling of non-linearities (splines, interactions) explained.
  • Missing data approach described (multiple imputation, complete case).

4. Overlap/Positivity

  • Visual inspection of PS distributions in both groups (histograms/KDE).
  • Overlap documented.
  • Trimming/truncation decisions reported if necessary.

5. Balance Assessment

  • Standardized Mean Differences (SMD) before/after adjustment reported.
  • Balance threshold declared (e.g., SMD < 0.1).
  • Love plots or tables included in supplement.
  • If balance not achieved: re-specify PS model.

6. Method Used

  • Method clearly stated: matching (caliper? ratio? replacement?), stratification, IPTW (stabilized? truncated?), covariate adjustment, or doubly robust (IPTW + regression).

7. Outcome Analysis

  • Correct model for outcome (logistic, Cox, linear).
  • Robust/clustered SEs if weighting or matching used.
  • Sensitivity analyses: model variations, trimming, alternative PS estimation.

8. Results & Interpretation

  • Effect estimate with CI and p-value reported.
  • Clinical relevance discussed, not only statistical significance.
  • Limitations highlighted: unmeasured confounding, positivity violations, sample size.

9. Reproducibility

  • Clear code (R or Python) shared if possible.
  • Workflow reproducible.
  • Dataset origin and cleaning pipeline described.

10. Transparent Limitations

  • Acknowledge PS is not magic: works only for observed confounders, does not replace randomization.

Suggested Bibliography (Key References)

Core papers

  • Rosenbaum PR, Rubin DB. The central role of the propensity score in observational studies for causal effects. Biometrika. 1983.
  • Austin PC. An Introduction to Propensity Score Methods for Reducing the Effects of Confounding. Multivariate Behavioral Research. 2011.
  • Austin PC, Stuart EA. Moving towards best practice when using inverse probability of treatment weighting. Stat Med. 2015.

Practice & tutorials

  • Stuart EA. Matching methods for causal inference: A review and a look forward. Stat Sci. 2010.
  • Hernán MA, Robins JM. Causal Inference: What If. Chapman & Hall/CRC. 2020 (free PDF online).

Diagnostics & reporting

  • Franklin JM, et al. Accuracy in Claims-Based Studies. Epidemiology. 2019.
  • Wang Y, et al. Love plots: visual analytics for PS balance. Stat Med. 2017.

Web Resources

  • UCLA Statistical Consulting — practical examples in Stata/R/SPSS: https://stats.idre.ucla.edu/
  • EpiLunch / Causal inference tutorials: https://www.theeffectbook.net/
  • Hernán & Robins, full book (free): https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/
  • Love plots in R (MatchIt, cobalt): https://cran.r-project.org/web/packages/cobalt/
  • Python IPTW and PS tutorials: https://www.pymc.io/projects/docs/en/stable/

Cite this article

Pierri, M. D. (2026). Propensity Score. micheledpierri.com. Permalink

Share:Email·LinkedIn
Previous← Survival Analysis
Statistics
  1. Statistics in Medicine
  2. Data Type
  3. Statistical Distributions
  4. Central Limit Theorem
  5. Descriptive Statistics
  6. Study Design in Medicine
  7. Inferential Statistics: Hypothesis Testing
  8. Choose the Right Statistical Test
  9. T test
  10. One-way ANOVA
  11. Repeated measure ANOVA
  12. Correlations
  13. Linear Regression
  14. Logistic Regression 1
  15. Logistic Regression 2
  16. Chi Square
  17. ANCOVA
  18. Bayes’ Theorem
  19. Monte Carlo Simulation
  20. Nonparametric statistics
  21. Survival Analysis
  22. Propensity Score

Leave a Reply Cancel reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.

© 2024–2026 micheledpierri.com · Privacy Policy · Impressum