Propensity Scores

32.1 One number that summarises selection (40 min)

Concept. The propensity score is each customer's probability of joining given their covariates: e(X) = P(joined = 1 | X). Rosenbaum and Rubin (1983) proved two things:

1. Balancing property: among customers with the same propensity score, the covariates have the same distribution in joiners and non-joiners. Comparing them is like a small randomised experiment.
2. Sufficiency: if adjusting for X removes confounding, so does adjusting for e(X) alone.

So instead of matching on six covariates, we can match on one number. The catch: e(X) must be estimated, and the result is only as good as the model and the covariates in it.

Two conditions make it work, and you should be able to name them:

32.2 Estimating the propensity score (50 min)

Type PS_FORMULA, propensity and smd into estimators.py. smd is the standardised mean difference from Day 23's balance table, now with optional weights.

Create day32/estimate_ps.py:

# day32/estimate_ps.py   (python -m day32.estimate_ps)
"""The propensity score: each customer's probability of joining, given what we know about them."""
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.formula.api as smf

from causalkit.estimators import PS_FORMULA, propensity
from causalkit.loyalty import public, world

full = world()
df = public(full).to_pandas()

model = smf.logit(PS_FORMULA, df).fit(disp=False)
print(model.summary().tables[1])

df["ps"] = propensity(df)
print("\nPropensity score by group (quantiles 1%, 25%, 50%, 75%, 99%):")
for g, label in [(1, "joined"), (0, "did not join")]:
    print(f"  {label:<13}", np.round(np.quantile(df.ps[df.joined == g], [0.01, 0.25, 0.5, 0.75, 0.99]), 3))

fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(df.ps[df.joined == 0], bins=50, alpha=0.6, label="did not join", color="#d95926")
ax.hist(df.ps[df.joined == 1], bins=50, alpha=0.6, label="joined", color="#3987e5")
ax.set_xlabel("Estimated propensity to join"); ax.set_ylabel("Customers"); ax.legend()
ax.set_title("Overlap: both groups appear at every propensity level")
fig.tight_layout(); fig.savefig("figures/day32_propensity_overlap.png", dpi=150)
print("Saved figures/day32_propensity_overlap.png")

python -m day32.estimate_ps

                        coef    std err          z      P>|z|      [0.025      0.975]
Intercept            -7.4894      0.249    -30.090      0.000      -7.977      -7.002
np.log(pre_spend)     1.0603      0.043     24.910      0.000       0.977       1.144
pre_orders            0.0002      0.004      0.048      0.962      -0.008       0.008
engaged               0.8460      0.033     25.954      0.000       0.782       0.910
tenure_years          0.1525      0.007     21.252      0.000       0.138       0.167
age                  -0.0089      0.001     -7.650      0.000      -0.011      -0.007
mobile                0.0275      0.032      0.852      0.394      -0.036       0.091

Propensity score by group (quantiles 1%, 25%, 50%, 75%, 99%):
  joined        [0.136 0.438 0.603 0.741 0.924]
  did not join  [0.061 0.228 0.366 0.529 0.844]
Saved figures/day32_propensity_overlap.png

Open the figure. The two histograms are shifted, which is selection, but they overlap over almost the whole range: for nearly every joiner there are non-joiners with a similar score. That's what makes the comparisons on Days 32–36 possible.

Don't read the logit coefficients causally, and don't drop "insignificant" variables like mobile. The propensity model's only job is to balance confounders, not to explain joining.

32.3 The balancing property (55 min)

Create day32/balancing_property.py:

# day32/balancing_property.py   (python -m day32.balancing_property)
"""
The balancing property: customers with the SAME propensity score look alike, whether they joined or not.
Check it within propensity-score quintiles.
"""
import pandas as pd
import polars as pl

from causalkit.estimators import COVARIATES, propensity, smd
from causalkit.loyalty import public, world

df = public(world()).to_pandas()
df["ps"] = propensity(df)
df["quintile"] = pd.qcut(df.ps, 5, labels=False)

rows = []
for col in COVARIATES:
    row = {"covariate": col,
           "overall": round(smd(df[col][df.joined == 1], df[col][df.joined == 0]), 3)}
    for q in range(5):
        s = df[df.quintile == q]
        row[f"Q{q + 1}"] = round(smd(s[col][s.joined == 1], s[col][s.joined == 0]), 3)
    rows.append(row)
print("Standardised mean differences, joiners vs non-joiners:")
print(pl.DataFrame(rows))

python -m day32.balancing_property

Standardised mean differences, joiners vs non-joiners:
│ covariate    ┆ overall ┆ Q1     ┆ Q2     ┆ Q3     ┆ Q4     ┆ Q5     │
│ pre_spend    ┆ 0.649   ┆ 0.189  ┆ 0.054  ┆ -0.006 ┆ 0.025  ┆ 0.32   │
│ pre_orders   ┆ 0.595   ┆ 0.1    ┆ 0.066  ┆ 0.017  ┆ -0.02  ┆ 0.271  │
│ engaged      ┆ 0.465   ┆ 0.071  ┆ -0.001 ┆ 0.015  ┆ 0.073  ┆ -0.009 │
│ tenure_years ┆ 0.553   ┆ 0.147  ┆ 0.006  ┆ 0.069  ┆ -0.033 ┆ 0.071  │
│ age          ┆ -0.097  ┆ -0.026 ┆ -0.0   ┆ 0.012  ┆ 0.014  ┆ -0.074 │
│ mobile       ┆ 0.017   ┆ 0.051  ┆ -0.025 ┆ -0.024 ┆ 0.009  ┆ 0.035  │

Overall, joiners differ by 0.5–0.65 standard deviations on four covariates. Within the middle quintiles of the propensity score, almost every difference falls under 0.1, the usual threshold. The edge quintiles (Q1, Q5) are wider bins, with a lot of variation inside them, so balance there is worse. Finer bins fix that, as the next section shows.

32.4 Subclassification on the propensity score (40 min)

Concept. Cut the score into strata, compute the joiner/non-joiner difference within each, and average the differences weighted by each stratum's share of everyone (ATE) or of joiners (ATT). It's Day 23's stratified estimator, with strata built from the score.

Create day32/ps_stratification.py:

# day32/ps_stratification.py   (python -m day32.ps_stratification)
"""
Mini-project: subclassification on the propensity score. Within each stratum compare joiners
with non-joiners, then average the stratum effects (weighted by stratum size for the ATE,
by number of joiners for the ATT). How many strata are enough?
"""
import pandas as pd
import polars as pl

from causalkit.estimators import propensity
from causalkit.loyalty import public, truth, world

full = world()
df = public(full).to_pandas()
t = truth(full)
df["ps"] = propensity(df)


def subclassify(df: pd.DataFrame, k: int) -> dict:
    s = df.assign(stratum=pd.qcut(df.ps, k, labels=False))
    g = s.groupby(["stratum", "joined"]).spend.mean().unstack()
    effect = g[1] - g[0]
    n_all = s.groupby("stratum").size()
    n_treated = s[s.joined == 1].groupby("stratum").size()
    return {"ATE": (effect * n_all).sum() / n_all.sum(), "ATT": (effect * n_treated).sum() / n_treated.sum()}


rows = [{"strata": k, **{e: round(v, 1) for e, v in subclassify(df, k).items()}} for k in [1, 5, 10, 20, 50]]
print(pl.DataFrame(rows))
print(f"Truth: ATE {t['ATE']:.1f}   ATT {t['ATT']:.1f}")

python -m day32.ps_stratification

│ strata ┆ ATE   ┆ ATT   │
│ 1      ┆ 593.0 ┆ 593.0 │
│ 5      ┆ 224.1 ┆ 281.4 │
│ 10     ┆ 194.5 ┆ 244.7 │
│ 20     ┆ 181.4 ┆ 227.1 │
│ 50     ┆ 173.5 ┆ 214.9 │
Truth: ATE 158.4   ATT 186.0

The textbook claim is that five strata remove about 90% of the bias. Here five strata remove about 85% (593 → 224 against a truth of 158), and even 50 strata still overstate the ATE by $15. Why? Within a stratum, joiners still have somewhat higher pre_spend, and spending is convex in pre_spend, so small leftover differences become large outcome differences. Subclassification is a good diagnostic and a rough first estimate, not a final one.

32.5 Day 32 mini-project: what a missing confounder does (65 min)

Suppose the CRM team forgot to export engaged (opens marketing emails). Re-estimate the propensity score without it and repeat the 20-strata ATE. Create day32/omit_engaged.py.

# day32/omit_engaged.py   (python -m day32.omit_engaged)
import pandas as pd

from causalkit.estimators import propensity
from causalkit.loyalty import public, truth, world

full = world()
df = public(full).to_pandas()
for name, f in [("full model", None),
                ("engaged omitted", "joined ~ np.log(pre_spend) + pre_orders + tenure_years + age + mobile")]:
    ps = propensity(df) if f is None else propensity(df, f)
    s = df.assign(stratum=pd.qcut(ps, 20, labels=False))
    g = s.groupby(["stratum", "joined"]).spend.mean().unstack()
    n = s.groupby("stratum").size()
    print(f"{name:<16} ATE with 20 strata: {((g[1] - g[0]) * n).sum() / n.sum():.1f}")
print(f"Truth: ATE {truth(full)['ATE']:.1f}")

full model       ATE with 20 strata: 181.4
engaged omitted  ATE with 20 strata: 205.9
Truth: ATE 158.4

Add a short "What we assume" section to day32/NOTES.md: list the confounders, and for each one say why it affects both joining and spending.

32.6 Save your work with Git (20 min)

git add .
git commit -m "Day 32: propensity score, balancing property, subclassification"
git push

Day 32 checklist

• I can define the propensity score and state the balancing property
• I can name unconfoundedness and overlap, and say which one the data can check
• I can estimate a propensity score and plot its overlap by group
• I can estimate ATE and ATT by subclassification, and explain why leftover imbalance remains
• I know that balance on measured covariates says nothing about unmeasured ones