Propensity Score Matching

Production engineering guide for propensity score matching covering patterns, implementation strategies, and operational best practices.

Propensity Score Matching (PSM) is a causal inference technique used to estimate the effect of a treatment on an outcome when randomization is not feasible. It balances covariates between treated and control groups by matching units based on their propensity scores—predicted probabilities of receiving treatment derived from observed covariates. PSM matters when observational data is used to answer causal questions, such as evaluating a new onboarding flow on user retention, a marketing campaign on conversion, or a feature rollout on engagement—where treatment assignment is non-random and confounded by measurable characteristics.

Estimating Propensity Scores

To begin PSM, compute the probability of treatment assignment using a model trained on pre-treatment covariates. The model must be well-calibrated and include all relevant confounders.

import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler


data = pd.read_csv("user_engagement.csv")

# Define covariates (X) and treatment (T)
X = data[["age", "signup_source", "device_type", "first_login_days", "avg_session_length"]]
T = data["treatment"]  # 1 if treated, 0 if control

# Fit logistic regression to estimate propensity scores
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

model = LogisticRegression(max_iter=1000)
model.fit(X_scaled, T)

# Predict propensity scores
data["propensity_score"] = model.predict_proba(X_scaled)[:, 1]

The propensity score is the fitted probability of receiving treatment. It is critical that the model includes all variables that affect both treatment assignment and outcome. Missing a key covariate leads to residual confounding, which PSM cannot correct.

Matching Strategy Selection

After estimating propensity scores, assign each treated unit to one or more control units based on proximity in score space. The choice of matching algorithm affects performance, memory, and stability.

Nearest Neighbor Matching

The most common approach: for each treated unit, find the control unit with the closest propensity score.

import numpy as np

def nearest_neighbor_match(treated, control, max_distance=None):
    """Match treated units to control units using nearest neighbor in propensity score space."""
    matched_pairs = []
    unmatched_control = control.copy()

    for idx, t_row in treated.iterrows():
        scores = np.abs(unmatched_control["propensity_score"] - t_row["propensity_score"])
        min_idx = np.argmin(scores)
        min_score = scores[min_idx]

        if max_distance is not None and min_score > max_distance:
            continue  # Skip if distance exceeds threshold

        matched_pairs.append({
            "treated_id": t_row.name,
            "control_id": unmatched_control.iloc[min_idx].name,
            "distance": min_score
        })

        unmatched_control = unmatched_control.drop(unmatched_control.iloc[min_idx].name)

    return matched_pairs, unmatched_control

The max_distance parameter—commonly defined as a fixed value (e.g., 0.2) or as a fraction of the standard deviation of the propensity scores—is critical. Without it, outliers in the propensity space can be matched poorly, leading to biased estimates.

Caliper Matching

To improve robustness, apply a caliper: a maximum allowed difference in propensity scores. Units are matched only if their scores differ by no more than the caliper width.

caliper = 0.2 * data["propensity_score"].std()
matched_pairs, remaining_control = nearest_neighbor_match(treated, control, max_distance=caliper)

A caliper of 0.2 standard deviations is common, but the optimal value depends on the data. Too narrow a caliper excludes many controls and increases variance; too wide, and the match loses precision.

Full Matching

When a single control unit can match multiple treated units, or vice versa, use full matching. This allows for variable numbers of matched units per group and reduces bias by using all available data.

# Assign each treated unit to all control units within caliper
def full_match(treated, control, caliper):
    matched_pairs = []
    control_used = set()

    for t_idx, t_row in treated.iterrows():
        candidates = control[
            (control["propensity_score"] >= t_row["propensity_score"] - caliper) &
            (control["propensity_score"] <= t_row["propensity_score"] + caliper)
        ]

        for c_idx, c_row in candidates.iterrows():
            matched_pairs.append({
                "treated_id": t_idx,
                "control_id": c_idx,
                "distance": abs(c_row["propensity_score"] - t_row["propensity_score"])
            })
            control_used.add(c_idx)

    return matched_pairs, control_used

Full matching increases computational cost and is more sensitive to overlapping propensity scores. It is especially useful when treatment is rare and control units are abundant.

Assessing Balance

After matching, verify that covariates are balanced between treated and control groups. A failure to balance is a silent but critical flaw.

Pre- and Post-Matching Balance Checks

Compute standardized mean differences (SMDs) before and after matching. SMDs > 0.1 indicate meaningful imbalance.

def compute_balance(data, matched_pairs, covariates):
    # Create matched dataset
    matched_data = data.loc[
        [p["treated_id"] for p in matched_pairs] +
        [p["control_id"] for p in matched_pairs]
    ].copy()

    matched_data["group"] = "treated"
    matched_data.loc[matched_data.index.isin([p["control_id"] for p in matched_pairs]), "group"] = "control"

    balance_results = {}
    for cov in covariates:
        treated_mean = matched_data[matched_data["group"] == "treated"][cov].mean()
        control_mean = matched_data[matched_data["group"] == "control"][cov].mean()
        pooled_sd = np.sqrt(
            ((len(matched_data[matched_data["group"] == "treated"]) - 1) * matched_data[matched_data["group"] == "treated"][cov].var() +
             (len(matched_data[matched_data["group"] == "control"]) - 1) * matched_data[matched_data["group"] == "control"][cov].var()) /
            (len(matched_data) - 2)
        )

        smd = (treated_mean - control_mean) / pooled_sd
        balance_results[cov] = smd

    return balance_results

Common failure modes:

The most common cause of imbalance is a mismatched set of covariates. For example, if avg_session_length is included in the propensity model but not in the outcome model, bias persists.

Estimating Treatment Effects

With matched pairs, compute the average treatment effect (ATE) or average treatment effect on the treated (ATT).

Simple Mean Difference

# Assume outcome column is 'retention_30d'
matched_data = pd.read_csv("matched_pairs.csv")

# Compute ATT
att = (
    matched_data[matched_data["group"] == "treated"]["retention_30d"].mean() -
    matched_data[matched_data["group"] == "control"]["retention_30d"].mean()
)

Weighted Estimation

Use inverse probability weighting (IPW) to account for multiple matches or varying quality. Each matched pair contributes a weight based on the number of times a unit appears in the matched sample.

# Count how many times each unit appears in matched pairs
treated_counts = pd.Series(
    [p["treated_id"] for p in matched_pairs]
).value_counts()
control_counts = pd.Series(
    [p["control_id"] for p in matched_pairs]
).value_counts()

# Compute weighted ATE
weighted_ate = 0
total_weight = 0

for pair in matched_pairs:
    t_id = pair["treated_id"]
    c_id = pair["control_id"]
    outcome_diff = (
        data.loc[t_id, "retention_30d"] -
        data.loc[c_id, "retention_30d"]
    )
    weight = 1 / (treated_counts[t_id] * control_counts[c_id])

    weighted_ate += outcome_diff * weight
    total_weight += weight

weighted_ate /= total_weight

This method reduces bias when units are matched to multiple partners and is especially robust when caliper width is varied across groups.

Common Failure Modes

Validating the PSM Pipeline

To ensure robustness, validate the entire PSM pipeline end-to-end.

# Validate data pipeline: check for missingness, data types, and distribution
python -c "
import pandas as pd
data = pd.read_csv('user_engagement.csv')
print('Missing values:', data.isnull().sum().sum())
print('Covariate distributions:', data[['age', 'signup_source']].describe())
"

# Validate model: check AUC on training and holdout set
python -c "
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split

X = data[['age', 'signup_source', 'device_type', 'first_login_days', 'avg_session_length']]
T = data['treatment']
X_train, X_test, T_train, T_test = train_test_split(X, T, test_size=0.2, random_state=42)
model = LogisticRegression().fit(X_train, T_train)
preds = model.predict_proba(X_test)[:, 1]
print('AUC:', roc_auc_score(T_test, preds))
"

# Validate matching: check caliper, balance, and output
python -c "
import pandas as pd
pairs = pd.read_csv('matched_pairs.csv')
data = pd.read_csv('user_engagement.csv')

# Compute balance
covariates = ['age', 'signup_source', 'device_type', 'first_login_days', 'avg_session_length']
balance = compute_balance(data, pairs.to_dict('records'), covariates)
print('Standardized Mean Differences:', balance)
"

The validation should be part of the CI/CD pipeline. A failed balance check (e.g., SMD > 0.2 for any covariate) should trigger a warning or block deployment.

When to Use PSM

Use PSM when:

Do not use PSM when:

PSM is not a black box—it is a process. The success of a PSM analysis depends on the clarity of the matching algorithm, the rigor of the balance checks, and the consistency of the data and model pipeline.

This page was rewritten on 17 October 2026. It replaced a templated version that was not specific to its own title. How these pages are written and checked. If anything here is wrong, tell us at [email protected] and we will correct it.