Propensity score matching

What this page answers

The cohort study chapter took the most direct route available: put the confounders into the Cox model together. That works well when the covariates are few and the events are many, but it has two limits.

First, it estimates “who gets treated” and “who has an event” inside the same model. Once there are many covariates you no longer know which patients are holding up the treatment coefficient. Two groups may not overlap at all on some covariate, and the model will still extrapolate a number for you along its linearity assumption, without warning.

Second, the number of events decides how many variables you may fit. The Cox model page gives the rule of thumb of about ten events per variable. This dataset has 47 events, which by that standard allows four or five variables — while the list of clinically worthwhile covariates runs to a dozen.

The propensity score addresses exactly that squeeze: compress a dozen covariates into a single number, so that the complexity of “adjustment” is no longer limited by the number of events, only by the number of treated patients.

What a propensity score is

The definition is one sentence: the probability that a person receives the exposure, given their covariates.

In practice it is estimated with logistic regression in which the outcome variable is treatment, not the clinical outcome:

It is useful because of a property that is not at all obvious: if a set of covariates suffices to remove confounding, then balancing this one score is enough to balance that entire set across the groups. A balancing problem in 12 dimensions collapses to one dimension. This is the balancing score property.

The worked example

We stay with the cohort from the cohort study chapter: patients undergoing radical prostatectomy who were transfused during surgery, where the exposure is the storage age of the transfused red cells (the oldest group versus short and medium storage) and the outcome is biochemical recurrence. Reusing the cohort is deliberate — one clinical question, answered two ways.

The dataset holds 316 patients. After dropping the 29 with a missing covariate, 287 enter the analysis (99 exposed, 188 controls), of whom 47 had a biochemical recurrence. The propensity score model contains 12 covariates.

Matching by hand, once

The code below is written from scratch in base R: logistic regression for the score, greedy nearest-neighbour matching, SMD computed from its definition. In practice MatchIt does all of this in one line (that line appears later on this page), but someone who has only ever seen that line does not know whom the caliper threw out.

library(medicaldata); library(survival)
data(blood_storage, package = "medicaldata")
set.seed(20260822)

d <- blood_storage
d$treat <- as.integer(d$RBC.Age.Group == 3)      # 1 = longest storage
covs <- c("Age", "AA", "FamHx", "PVol", "TVol", "T.Stage", "bGS",
          "PreopPSA", "PreopTherapy", "Units", "sGS", "AnyAdjTherapy")
cc <- d[complete.cases(d[, c("TimeToRecurrence", "Recurrence", "treat", covs)]), ]

# -- 1. Propensity score: the outcome variable is TREATMENT, not the endpoint --
ps_fit <- glm(reformulate(covs, "treat"), data = cc, family = binomial())
cc$ps  <- fitted(ps_fit)
cc$lps <- qlogis(cc$ps)          # match on the logit scale so the tails do not pile up

# -- 2. Greedy nearest neighbour, 1:1, without replacement, caliper = 0.2 logit SD --
cal  <- 0.2 * sd(cc$lps)
tr   <- sample(which(cc$treat == 1))   # random order: greedy matching is order-sensitive
pool <- which(cc$treat == 0)
pairs <- list()
for (i in tr) {
  if (!length(pool)) break
  dist <- abs(cc$lps[pool] - cc$lps[i])
  m <- which.min(dist)
  if (dist[m] > cal) next              # caliper blocks it: no match, so no pair
  pairs[[length(pairs) + 1]] <- c(i, pool[m])
  pool <- pool[-m]                     # without replacement: a used control leaves the pool
}
matched <- cc[unlist(pairs), ]
matched$pair <- rep(seq_along(pairs), each = 2)

# -- 3. SMD from its definition; the denominator stays the PRE-matching pooled SD --
smd_denom <- sapply(covs, function(v) {
  x <- cc[[v]]; t <- cc$treat == 1
  sqrt((var(x[t]) + var(x[!t])) / 2)
})
smd <- function(dat) sapply(covs, function(v) {
  x <- dat[[v]]; t <- dat$treat == 1
  (mean(x[t]) - mean(x[!t])) / smd_denom[[v]]
})
round(cbind(before = smd(cc), after = smd(matched)), 3)

# -- 4. Effect estimate, with the variance acknowledging the pairs --
coxph(Surv(TimeToRecurrence, Recurrence) ~ treat + cluster(pair), data = matched)

import numpy as np, pandas as pd, statsmodels.api as sm
from lifelines import CoxPHFitter

bs = sm.datasets.get_rdataset("blood_storage", "medicaldata").data
bs["treat"] = (bs["RBC.Age.Group"] == 3).astype(int)
covs = ["Age", "AA", "FamHx", "PVol", "TVol", "T.Stage", "bGS",
        "PreopPSA", "PreopTherapy", "Units", "sGS", "AnyAdjTherapy"]
cc = bs[["TimeToRecurrence", "Recurrence", "treat"] + covs].dropna().copy()

ps = sm.Logit(cc["treat"], sm.add_constant(cc[covs])).fit(disp=0).predict()
cc["lps"] = np.log(ps / (1 - ps))

rng = np.random.default_rng(20260822)
cal = 0.2 * cc["lps"].std()
tr = rng.permutation(cc.index[cc.treat == 1])
pool = list(cc.index[cc.treat == 0])
pairs = []
for i in tr:
    if not pool:
        break
    dist = (cc.loc[pool, "lps"] - cc.loc[i, "lps"]).abs()
    j = dist.idxmin()
    if dist[j] > cal:
        continue
    pairs.append((i, j)); pool.remove(j)

idx = [k for p in pairs for k in p]
matched = cc.loc[idx].copy()
matched["pair"] = np.repeat(np.arange(len(pairs)), 2)
CoxPHFitter().fit(matched[["TimeToRecurrence", "Recurrence", "treat", "pair"]],
                  duration_col="TimeToRecurrence", event_col="Recurrence",
                  cluster_col="pair").print_summary()

What comes out: 91 pairs, 182 people in total, with 8 exposed patients excluded because no control lay close enough inside the caliper.

Overlap: the figure to look at before matching

Before asking how well the matching went, there is a prior question: do the two groups’ propensity scores overlap at all? Where they do not overlap, no statistical method can supply a comparable control — any estimate there rests on the model extrapolating.

In this dataset the propensity scores run from 0.125 to 0.613, and the region of common support — the interval where both groups have people — is 0.233 to 0.609. Anyone outside common support is discarded automatically by matching, and that is a genuine advantage of matching over regression adjustment: regression will not tell you it is extrapolating; matching hands you a list of the people it could not match.

Judging balance: SMD, not p-values

Once matching is done, the one thing to check is balance. The criterion is the standardised mean difference (SMD), and by convention |SMD| below 0.1 counts as balanced. The definition of the SMD, why it does not inflate with sample size, and where the 0.1 line came from are covered in full on Table 1 and the standardised mean difference, so they are not repeated here.

One thing is worth stressing: do not run between-group significance tests before and after matching. Matching shrinks the sample, so p-values rise on their own — the sentence “after matching, the difference was no longer significant” holds even when balance did not improve at all. What that test measures is sample size, not balance.

Before matching, 5 covariates had |SMD| above 0.1, the largest being 0.156. After matching, 0 remain above the threshold, the maximum falls to 0.067, and the mean falls from 0.067 to 0.029.

Double robustness: why adjusting again is not double-counting

The previous section said that an important covariate still imbalanced after matching should be adjusted for again in the post-matching model. That advice has a name: double robustness.

Matching and weighting need the propensity score model to be right; the post-matching regression needs the outcome model to be right. Either can be wrong. The doubly robust property is that the estimate stays consistent as long as at least one of the two is correct — you do not have to bet on which one, you hold two tickets and need only one to win.

This is why “do not put covariates back in after matching, that would be adjusting twice” is wrong. Adjusting again costs a few degrees of freedom and buys a second chance at getting it right. In practice it is the default, not a remedy for special cases.

The three choices inside “matching”

Matching is not one action; it is three independent decisions. Here are three configurations on the same data:

What each decision buys:

• 1:1 or 1:k. When controls are plentiful, 1:k raises the number of events and therefore the precision. But the k-th nearest neighbour is always further away than the first, so the larger k is, the worse the balance — which is what the third row shows: 1:2 pushes the maximum |SMD| past 0.1.
• The caliper (the matching radius). The convention is 0.2 standard deviations of the logit propensity score (here, 0.067). What it blocks are the exposed cases for whom nobody similar enough exists. The first row discards 8 of them and buys better balance; the second sets no caliper, discards nobody, and ends up with a larger maximum |SMD|. This is a direct trade of bias against sample size.
• With or without replacement. The code above matches without replacement: a control that has been used leaves the pool. Allowing replacement lets one control serve several exposed patients, which usually improves balance, but the same person is then counted more than once and the variance has to be handled specially — and if a handful of controls are reused repeatedly, the effective sample size falls far below the nominal count.

The last two columns give each configuration’s HR and confidence interval (all with pair-robust variance). They are listed not so that you can pick one, but because this is what the “report all of them” rule demanded in the next section actually looks like: three point estimates from 0.83 to 1.03, with heavily overlapping intervals, all of which cross 1. The one that sits furthest from 1 is 1:2, caliper 0.2 SD (0.83, 0.40–1.74), and it is simultaneously the configuration that discarded the most exposed patients (23) and achieved the worst balance (maximum |SMD| 0.102; covariates past the line: 1) — the configuration whose point estimate moved most is also the one whose sample was replaced most.

The way to read a sensitivity analysis is to see how far the numbers spread, not to choose one of them. Here they do not spread far, but all three intervals are wide enough that a moderate effect cannot be excluded, so the narrow spread does not establish much either.

The matched sample is no longer the population you started with

This is the point most often skipped, and it directly determines how the conclusion may be phrased.

1:1 matching, especially with a caliper, keeps every exposed case for whom a control could be found, plus the controls most like them. So the matched sample represents “people who look like the exposed group”, and the quantity being estimated is the ATT (average treatment effect on the treated), not the ATE for the whole population.

When the effect varies with patient characteristics (effect modification) these three differ numerically, and they answer different clinical questions. ATT asks “were these people, already on the drug, right to be on it?”; ATE asks “what would happen if we rolled this out to everyone?”. The inverse probability weighting page shows all three estimands computed on this same dataset.

The model on this page discarded 8 exposed cases, so strictly speaking it is not even the ATT — it is the ATT among the exposed patients who could be matched. That sentence belongs in the limitations, together with the number of people dropped.

Variance after matching

The two rows of a matched pair are not independent — they were put together because their propensity scores were close. In principle that means the variance estimate has to account for the pairing, or the confidence interval will be off. The way to do it is to add cluster(pair) to the Cox model, which switches to a robust (sandwich) variance.

The four rows of this table are not the same estimand, so do not use the size of the HR to decide which method is better. “Multivariable Cox (all patients)” is a conditional HR in the full sample, given the covariates; “after 1:1 matching” is an HR in the matched sample targeting the ATT. The estimation target and the target population both differ — and on top of that, the Cox HR is non-collapsible, so even with no residual confounding at all, the conditional and marginal values need not agree. A difference between the numbers is not in itself evidence that matching improved the estimate.

The one-line version, with MatchIt

The thirty hand-written lines above are, in practice, this:

library(MatchIt)
m <- matchit(treat ~ Age + AA + FamHx + PVol + TVol + T.Stage + bGS +
               PreopPSA + PreopTherapy + Units + sGS + AnyAdjTherapy,
             data = cc, method = "nearest", distance = "glm",
             caliper = 0.2, std.caliper = TRUE, replace = FALSE)
summary(m)                 # balance table, with SMDs before and after
md <- match.data(m)        # matched data, carrying subclass and weights columns
coxph(Surv(TimeToRecurrence, Recurrence) ~ treat + cluster(subclass), data = md)

MatchIt 4.7.2 produces 92 pairs and an HR of 1.04 (0.51–2.13), against the hand-written version’s 91 pairs and HR of 1.02; the difference is nothing but the random order in which greedy matching processed the exposed cases. The two not agreeing exactly is expected — greedy matching is sensitive to who is handled first, which is why a matching script must always set a seed.

In real work, use the package: it handles optimal matching, full matching, weights for matching with replacement, and a complete set of balance diagnostics. The point of the hand-written version is to let you see what the package is doing, and to know where to look when it gives you something strange.

How to write up the result

All three models put the point estimate near 1, with every confidence interval crossing it. The correct wording is: in this cohort, no association was detected between red-cell storage age and biochemical recurrence.

It may not be written as “storage age is unrelated to recurrence”, and the matched HR of 1.02 may not be described as “an increase in risk” — the interval 0.51–2.03 accommodates substantial benefit and substantial harm at the same time, and with this many events (32 after matching) the data cannot exclude a moderate effect.

Common misuses

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B6-02-psm.R

Watch next

▶資料小探 – 傾向分數配對法

▶Propensity scores: Everything you need to know in 5min

▶How Propensity Scores Work | NEJM Evidence

▶Propensity score matching: an introduction

▶Outcome research: Causal inference & Propensity score II 傾向分數

Report a problem on this page