---
title: Propensity Score
date: 2026-06-07T06:03:09Z
modified: 2026-06-07T06:03:13Z
permalink: "https://www.micheledpierri.com/statistics/propensity-score/"
type: page
status: publish
excerpt: ""
wpid: 2370
featured_image: "https://www.micheledpierri.com/wp-content/uploads/2026/06/Propensity-score_2.png"
featured_image_alt: a young gardener prunes a row of uneven trees, shaping them into uniform ones
timestamp: 2026-06-07T06:03:13Z
tags: []
---
# 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:
-  for treated patients,  for controls
-  the vector of baseline covariates (age, sex, comorbidities, risk scores…)
then the propensity score is:
In practice,  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 (
):
- For controls (
):
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:
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
```
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
})
```
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
```
# 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
```
# 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()
```
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
```
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))
```
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()
```
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()
```
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
```
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})")
```
# 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})")
```
# 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)
```
# 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)
```
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/](https://stats.idre.ucla.edu/)
- **EpiLunch / Causal inference tutorials**: [https://www.theeffectbook.net/](https://www.theeffectbook.net/)
- **Hernán & Robins, full book (free)**: [https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/](https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/)
- **Love plots in R (MatchIt, cobalt)**: [https://cran.r-project.org/web/packages/cobalt/](https://cran.r-project.org/web/packages/cobalt/)
- **Python IPTW and PS tutorials**: [https://www.pymc.io/projects/docs/en/stable/](https://www.pymc.io/projects/docs/en/stable/)