Propensity Score Methods (PSM, IPTW)

In plain language

A propensity score is a single number — the probability that a patient would receive the study drug, given everything measurable about them at baseline (age, diagnoses, prior medications, etc.). Once every patient has that score, you can pair treated patients with comparators who had nearly the same probability of being treated, or you can reweight the groups so their measured characteristics look alike. The goal is to make the two treatment arms resemble each other the way randomization would, so that any difference in outcomes can be attributed to the drug rather than to the kinds of patients who happened to receive it. This approach can only account for differences you actually measured — confounders hidden in the data, like disease severity not recorded in claims, are not removed.

The propensity score is the conditional probability of receiving the treatment of interest given measured baseline covariates, e(X)=Pr(A=1|X). Under the identification assumptions of conditional exchangeability (no unmeasured confounding given X), positivity (0<e(X)<1 for every covariate pattern), and consistency, conditioning on the scalar e(X) recovers the same covariate balance as conditioning on the full vector X (the balancing property of Rosenbaum and Rubin). In pharmacoepidemiology the PS is the standard tool for making initiators of two treatments comparable on measured confounders. It is an analysis step layered on top of design — it presupposes that time zero, eligibility, new-user status, the active comparator, and the outcome/censoring rules are already correct. A perfectly balanced PS analysis built on a wrong time zero is still biased; the PS cannot repair immortal time, depletion of susceptibles, or a comparator that itself causes the outcome.

Core estimand distinction

. The choice of PS method is the choice of estimand, and the four common methods target four different populations and answer four different questions: - ATT (average treatment effect in the treated) — what would happen to the treated population if untreated. Targeted by 1:1 PS matching (matched-set ATT) and by SMR weighting (treated weight 1, comparator weight e(X)/(1-e(X))). This is usually the policy-relevant estimand for a safety signal of a drug already in use. - ATE (average treatment effect) — what would happen to the whole eligible population if all vs none were treated. Targeted by (stabilized) IPTW: treated weight Pr(A=1)/e(X), comparator weight Pr(A=0)/(1-e(X)). Most sensitive to near-positivity violations because tail subjects with e(X)→0 or 1 receive enormous weights. - ATO (average treatment effect in the overlap population) — the effect in patients with genuine clinical equipoise. Targeted by overlap weighting (treated weight 1-e(X), comparator weight e(X)); exactly balances every covariate mean, yields the smallest variance among balancing weights, and never produces extreme weights. The estimand is a smoothly-weighted population, not a fixed clinical group, which must be stated. - ATM (matching weights) — a weighting analogue of pair matching targeting the "matchable" subpopulation. None of these is a free choice cosmetically: ATT, ATE, and ATO can differ in magnitude and even sign when effect modification by the PS is present, so the estimand must be pre-specified, not chosen after seeing results.

Pros, cons, and trade-offs

- vs covariate adjustment / outcome regression on X: PS reduces the confounder set to one dimension, decouples the "design" stage (fit and balance the PS, blinded to outcomes) from the "analysis" stage, and lets you assess overlap and balance explicitly before touching the outcome — discipline that pure outcome regression hides. Cost: outcome regression can be more efficient when the outcome model is correct and is not derailed by positivity violations. Prefer PS when you want an outcome-blind design stage and transparent balance diagnostics in a regulated study. - vs doubly-robust estimation (AIPW / TMLE; see g-estimation and MSM entries): single-PS methods are simpler to specify, explain to reviewers, and audit. Cost: they are consistent only if the PS model is correct — a single line of defense. Doubly-robust estimators stay consistent if either the PS or the outcome model is correct and achieve smaller asymptotic variance. Prefer plain PS for transparency and when the confounder structure is well understood; escalate to doubly-robust/g-methods for time-varying confounding affected by prior treatment, which PS cannot handle. - vs high-dimensional PS (hdPS): an investigator-specified PS is interpretable and easy to defend clause-by-clause. Cost: it can miss empirical proxy confounders captured by hundreds of claims codes. Prefer hdPS (or hdPS as a sensitivity analysis) in large claims data where unmeasured confounding by proxy is plausible. - matching vs weighting, head to head: matching yields a concrete, inspectable matched cohort and is robust to extreme e(X) (poor-overlap subjects are simply unmatched), but it discards data, can depend on greedy sort order, and shifts the estimand to the matched region. Weighting uses all data and targets a clean estimand, but inherits the full danger of extreme weights (IPTW) unless you use stabilized weights, truncation, or overlap weights. Overlap weighting is the modern default when positivity is shaky.

When to use

. Comparative effectiveness or safety of two (or more) treatments where confounding is by measured baseline covariates, follow-up starts at a correctly defined time zero, and you can demonstrate overlap. PS is the standard balancing step after an active-comparator new-user cohort is built. Use matching or SMR weighting when the question is about the treated (a safety signal for current users); stabilized IPTW for a population (ATE) effect with acceptable overlap; overlap weighting when overlap is limited or extreme weights threaten IPTW; stratification or PS as a covariate only as quick checks, because they balance less completely and obscure overlap.

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

- Unmeasured or mismeasured confounding dominates. PS balances only what is in X. In claims, frailty, disease severity, BMI, smoking, and over-the-counter use are typically unmeasured; a beautifully balanced PS table provides false reassurance. Pair PS with negative-control outcomes, E-values, or quantitative bias analysis — do not present a balance table as proof of no confounding. - Positivity is violated. If one drug is reserved for sicker or renally-impaired patients, e(X) piles up near 0/1, IPTW weights explode, matching discards most of one arm, and the surviving estimand no longer maps to any clinically meaningful population. Trimming/restriction or overlap weighting changes the question; say so. - Time-varying confounding affected by prior treatment. Conditioning on post-baseline confounders that are also mediators induces collider/over-adjustment bias; a baseline PS cannot fix this. Use marginal structural models with IPTW over time or g-estimation instead. - Treating the PS as a substitute for design. Using post-index variables in the PS, or applying PS to a prevalent-user cohort with immortal time, launders a biased cohort into a balanced-looking one — the most dangerous failure mode because the diagnostics look clean. - Few outcome events. With sparse events, weighting inflates variance and instability; matched or stratified analyses with firth-penalized or exact methods are safer.

Data-source operational depth

- Claims (FFS vs MA vs commercial): Covariates are built from demographics, calendar time, baseline diagnoses/procedures, prior drug classes, utilization intensity, and cost in the lookback window. Failure mode: Medicare Advantage person-time lacks fee-for-service encounter claims, so confounders measured from claims are systematically under-captured for MA enrollees — the PS is fit on degraded covariates and balance looks fine while real confounding persists. Restrict to enrollees with complete Parts A/B/D (or commercial medical+pharmacy) and exclude MA-only person-time. Failure mode: in elderly cohorts, death is a competing risk that often differs by exposure; a PS-weighted Kaplan-Meier or cause-specific hazard answers a different question than a subdistribution (Fine-Gray) analysis, and the two diverge — pre-specify which (see competing-risks entry). Failure mode: procedure-anchored index dates can smuggle in immortal time that no PS will remove. - EHR: Labs, vitals, smoking, and BMI sharpen the PS but are missing not-at-random — captured only when ordered. A missingness indicator vs multiple imputation are different estimands of e(X); site and encounter-frequency are themselves confounders (sicker, more-monitored patients have more covariate data). Include site/visit-intensity terms and treat loss to follow-up as potentially informative. - Registry: Strong on clinical confounders (stage, biomarker, ECOG, device details) that improve exchangeability, but weak on complete pharmacy exposure — link to claims for fills and to a death index for censoring before fitting the PS. - Linked claims–EHR–vital records: The ideal substrate, but order vs fill vs service-date discrepancies must be reconciled before fixing time zero, or baseline covariates leak across the index boundary and contaminate e(X).

Worked claims example

Question: risk of hospitalized heart failure with an SGLT2 inhibitor vs a DPP-4 inhibitor in a commercial + Medicare FFS database (ATT, because the signal concerns SGLT2i initiators). (1) Build the cohort: adults with type 2 diabetes, 365 days of continuous medical+pharmacy FFS enrollment before the first qualifying `fill_date`, and no SGLT2i or DPP-4i fill in that 365-day washout (incident users of both arms); index_date = that first fill; arm = the NDC dispensed on index_date; exclude MA-only person-time. (2) Measure covariates only in `[index_date-365, index_date]`: age, sex, baseline HbA1c proxies, prior insulin, CKD/HF/MI diagnosis codes, prior-year inpatient count, and total paid cost. (3) Fit e(X) by logistic regression (or hdPS). (4) Because the estimand is ATT, use SMR weighting (treated=1, comparator=e/(1-e)) or 1:1 caliper matching on the logit-PS (caliper 0.2 SD); confirm every post-balancing standardized mean difference <0.1 and inspect the e(X) overlap plot. (5) Define follow-up from index_date to first validated HF hospitalization, censoring at disenrollment, death, end of data, and — for an as-treated estimand — the end of the last `days_supply` plus a 30-day grace period or a switch to the other class. (6) Estimate the effect with a weighted Cox model using robust (sandwich) standard errors, report the effective sample size, and run sensitivity analyses on washout length, weight truncation at the 1st/99th percentile, and a negative-control outcome to probe residual confounding. Reporting must include the PS model specification, pre/post balance, the overlap plot, the weight distribution and ESS, the trimming rule, the targeted estimand, and whether outcome modeling was additionally used (never use post-index variables in the PS unless you are explicitly fitting a longitudinal g-method).

Interpreting the output

Consider a PSM analysis of Drug A versus Drug B — the readmission scenario above — reporting HR = 0.82 (95% CI 0.71–0.95) for 30-day hospital readmission after matching on 28 covariates, with post-match SMDs all below 0.10.

Formal interpretation: Among matched Drug A initiators — the subset of treated patients for whom a comparably similar Drug B patient existed in the data — the instantaneous readmission rate was 18% lower (HR 0.82, 95% CI 0.71–0.95). This is the average treatment effect in the treated (ATT) within the matched population, not among all Drug A initiators. Under IPTW targeting the full target population, the analogous estimate applies to the inverse-probability-weighted pseudo-population (ATE) or the ATT-weighted pseudo-population, depending on the weight specification. Both are valid causal estimates only under two untestable conditions: no unmeasured confounding — every covariate that jointly predicts drug choice and readmission is captured and balanced — and positivity — every treated patient had a realistic probability of receiving the comparator. The 95% CI means results this extreme would arise in fewer than 5% of identically designed studies if the true HR were 1.0.

Practical interpretation: Drug A initiators who closely resembled Drug B initiators had readmissions arriving at an 18% slower rate during follow-up. This is not a statement that 18% fewer total events occurred — the HR is not a risk ratio or a cumulative probability. Under PSM, patients without a suitable match are excluded; under IPTW, they are down-weighted. The finding speaks to the matched or reweighted population, not all Drug A users. Unmeasured factors — such as disease severity not captured in claims — remain uncontrolled and could bias the estimate in either direction.

Worked example

Scenario

A researcher wants to compare 30-day hospital readmission rates between patients who started Drug A (a newer diabetes medication) versus Drug B (an older one). Before comparing outcomes, they need to check whether the two groups look similar on key baseline characteristics. The raw data show that Drug A patients tend to be older and have more diabetes complications — a direct reflection of prescribing patterns, not drug effects. The researcher fits a propensity score model (logistic regression predicting Drug A vs. Drug B from age, diabetes diagnosis, and other baseline covariates) and then either matches each Drug A patient to the most similar Drug B patient, or reweights the sample so the groups balance. The table below shows covariate balance before and after that adjustment.

Dataset

Balance table: two covariates before and after propensity-score matching/weighting. SMD = standardized mean difference; values below 0.1 indicate acceptable balance.

Steps

• Before adjustment, Drug A patients are on average 6.2 years older than Drug B patients (67.4 vs 61.2), giving an SMD of 0.52 — a large imbalance by any standard.

• The diabetes diagnosis rate is also much higher in Drug A patients (78% vs 54%), SMD 0.50, reflecting that sicker patients with more comorbidities were more likely to be prescribed the newer drug.

• A logistic regression model is fit predicting Drug A vs Drug B from all baseline covariates; each patient's predicted probability from that model is their propensity score.

• Patients are then matched (each Drug A patient is paired with the Drug B patient whose propensity score is closest) or reweighted (each patient receives a weight inversely proportional to their probability of receiving the arm they actually received).

• After adjustment, mean age in the two arms is 64.8 vs 64.3 years (SMD = 0.04) and diabetes diagnosis rates are 71% vs 69% (SMD = 0.04) — both well below the 0.1 threshold.

• SMD arithmetic check for age after adjustment: difference = 64.8 - 64.3 = 0.5 years; pooled SD estimated at approximately 12.5 years; SMD = 0.5 / 12.5 = 0.04, consistent with the table.

Result

After propensity-score adjustment both covariates achieve SMD < 0.1, indicating the treated and comparator arms are now balanced on these measured characteristics. The researcher can proceed to compare readmission rates in the balanced sample, knowing that age and diabetes burden are no longer driving any observed difference. Note: unmeasured characteristics — such as disease severity not coded in claims — are still uncontrolled.

Runnable example

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from lifelines import CoxPHFitter

xcols = ["age", "sex", "cci", "ckd", "prior_hf", "prior_insulin", "prior_hosp", "prior_cost"]

# --- Stage 1 (design, outcome-blind): estimate e(X) = Pr(treated | X) ---
ps = LogisticRegression(max_iter=2000, C=1e6)  # near-unpenalized logistic PS
ps.fit(df[xcols], df["treated"])
df["e"] = np.clip(ps.predict_proba(df[xcols])[:, 1], 1e-6, 1 - 1e-6)

p = df["treated"].mean()
# Stabilized IPTW -> ATE (average treatment effect in the whole eligible population)
df["w_ate"] = np.where(df["treated"] == 1, p / df["e"], (1 - p) / (1 - df["e"]))
df["w_ate"] = df["w_ate"].clip(upper=df["w_ate"].quantile(0.99))  # truncate extreme weights
# Overlap weights -> ATO (clinical-equipoise population; no extreme weights by construction)
df["w_ato"] = np.where(df["treated"] == 1, 1 - df["e"], df["e"])
# SMR weights -> ATT (effect in the treated; comparator reweighted to look like treated)
df["w_att"] = np.where(df["treated"] == 1, 1.0, df["e"] / (1 - df["e"]))

def weighted_smd(x, t, w):
    m1 = np.average(x[t == 1], weights=w[t == 1]); m0 = np.average(x[t == 0], weights=w[t == 0])
    v1 = np.average((x[t == 1] - m1) ** 2, weights=w[t == 1])
    v0 = np.average((x[t == 0] - m0) ** 2, weights=w[t == 0])
    return (m1 - m0) / np.sqrt((v1 + v0) / 2)

wcol = "w_ate"  # pick the estimand here; balance below should be checked on the SAME weights
bal = {c: weighted_smd(df[c].values, df["treated"].values, df[wcol].values) for c in xcols}
ess = df[wcol].sum() ** 2 / (df[wcol] ** 2).sum()  # effective sample size after weighting
print("max |SMD| =", max(abs(v) for v in bal.values()), "| ESS =", round(ess, 1))
assert max(abs(v) for v in bal.values()) < 0.1, "covariate imbalance remains; revisit the PS model"

# --- Stage 2 (analysis): marginal weighted Cox with robust (sandwich) SEs ---
cox = CoxPHFitter()
cox.fit(df[["time", "event", "treated", wcol]], duration_col="time", event_col="event",
        weights_col=wcol, robust=True)  # robust=True accounts for the estimated weights
print(cox.summary[["coef", "exp(coef)", "se(coef)", "p"]])

library(WeightIt); library(cobalt); library(survey)

xform <- treated ~ age + sex + cci + ckd + prior_hf + prior_insulin + prior_hosp + prior_cost

# --- Estimate the PS and weights; estimand sets the target (ATE shown; use "ATT" or "ATO" to switch) ---
w <- weightit(xform, data = df, method = "ps", estimand = "ATE", stabilize = TRUE)

# --- Balance + overlap diagnostics on the weighted sample (must precede any outcome look) ---
bal.tab(w, un = TRUE, thresholds = c(m = .1))   # standardized mean differences pre/post weighting
love.plot(w, abs = TRUE, thresholds = c(m = .1))
cat("ESS:", summary(w)$effective.sample.size, "\n")

# --- Marginal weighted Cox model; survey design gives robust variance for the weighted estimator ---
df$ipw <- w$weights
des <- svydesign(ids = ~1, weights = ~ipw, data = df)
fit <- svycoxph(Surv(time, event) ~ treated, design = des)
summary(fit)   # exp(coef) = weighted hazard ratio for the chosen estimand