倾向得分匹配：原理、步骤与Python实现

!Propensity-Score-Matching

在数据科学和因果推断的领域里，我们常常面临一个核心挑战：如何从充满了偏差的观测数据中提取出可靠的因果关系？特别是在2026年的今天，随着企业对数据驱动决策的要求从“是什么”转向“为什么”，传统的简单对比已经无法满足需求。在这篇文章中，我们将深入探讨 倾向得分匹配 这一强大的统计工具，并结合我们最新的工程实践和AI辅助开发流程，带你看看如何在现代技术栈中高效地实现它。

倾向得分的核心概念

让我们先回到基础。倾向得分是指在给定一组观测协变量的情况下，个体接受治疗（处理）的概率。从数学上讲，倾向得分定义如下：

> e(X) = P(T = 1 \mid X)

其中：

• T = 1 表示该个体在治疗组中。
• X 是观测协变量的向量。
• e(X) 是倾向得分。

观测数据中的隐性陷阱

在我们最近的一个用户增长项目中，我们遇到了一个非常典型的问题。在随机实验（A/B测试）中，治疗通常是随机分配的，这保证了组间的公平性。但在观测数据中，治疗往往不是随机的——用户可能会根据年龄、收入、活跃度等因素进行自我选择。

这种情况会导致混淆，从而使治疗效果估计产生偏差。比如，我们发现使用新功能（T=1）的用户留存率更高，但这可能是因为这些用户本身就是高活跃度的忠实用户，而不是因为功能本身有效。

解决方案：基于倾向得分进行匹配

我们不需要在每个具体的协变量（如年龄、收入、所在城市）上精确地匹配个体——这在变量维度很高时不仅计算量大，而且几乎是不可能做到的。相反，我们可以根据倾向得分——一个汇总了所有协变量信息的单一数值——来进行匹配。这就是 Rosenbaum 和 Rubin 提出的伟大之处。

现代工程实战：构建生产级 PSM 流水线

在2026年的开发环境中，我们编写代码的方式已经发生了深刻的变化。当我们着手实现 PSM 时，我们不再仅仅是写一个脚本，而是在构建一个可维护、可解释的因果推断流水线。

1. 倾向得分的现代估计方法

虽然逻辑回归是最经典的方法，但在处理复杂非线性关系时，它可能显得力不从心。在我们的工具箱中，梯度提升机 和 深度学习 模型正成为主力。

让我们看看如何使用 Python 结合现代工具链执行一个更健壮的 PSM。

import pandas as pd
import numpy as np
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.neighbors import NearestNeighbors
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

# 模拟更接近真实生产环境的数据（包含非线性关系）
np.random.seed(42)
n_samples = 2000
data = pd.DataFrame({
    ‘age‘: np.random.normal(35, 10, n_samples),
    ‘income‘: np.random.normal(60000, 20000, n_samples),
    ‘credit_score‘: np.random.normal(650, 100, n_samples),
    ‘tenure_months‘: np.random.randint(0, 120, n_samples)
})

# 这里的接受治疗概率是非线性的，仅用逻辑回归可能欠拟合
def treatment_assignment(row):
    base_prob = 0.3
    if row[‘age‘] > 30: base_prob += 0.2
    if row[‘income‘] > 50000: base_prob += 0.1
    if row[‘credit_score‘] > 700: base_prob += 0.2
    return 1 if np.random.binomial(1, min(base_prob, 0.99)) == 1 else 0

data[‘treatment‘] = data.apply(treatment_assignment, axis=1)

# 结果变量：包含异质性处理效应
data[‘outcome‘] = (0.5 * data[‘age‘] + 
                   0.0002 * data[‘income‘] + 
                   3.0 * data[‘treatment‘] +  # 真实的 ATE 是 3.0
                   0.5 * data[‘treatment‘] * (data[‘age‘] > 40).astype(int) + # 交互效应
                   np.random.normal(0, 5, n_samples))

# --- 步骤 1: 使用现代 ML 模型估计倾向得分 ---
features = [‘age‘, ‘income‘, ‘credit_score‘, ‘tenure_months‘]
X = data[features]
y = data[‘treatment‘]

# 使用 GBM 代替简单的 Logistic Regression，捕捉非线性特征
ps_model = GradientBoostingClassifier(random_state=42)
ps_model.fit(X, y)
data[‘propensity_score‘] = ps_model.predict_proba(X)[:, 1]

print("倾向得分统计:")
print(data[‘propensity_score‘].describe())

2. 优化匹配策略：处理边界情况

在上述代码中，你可能已经注意到，我们使用了更复杂的模型来生成得分。接下来是匹配环节，也是我们最容易踩坑的地方。在2026年的工程实践中，我们强烈建议使用 卡尺匹配 结合 有放回的最近邻匹配，以防止引入糟糕的匹配。

# --- 步骤 2: 高级匹配策略 ---

def perform_matching(df, treatment_col=‘treatment‘, ps_col=‘propensity_score‘, caliper=0.05):
    treated = df[df[treatment_col] == 1].copy()
    control = df[df[treatment_col] == 0].copy()
    
    # 使用 NearestNeighbors 进行快速近似搜索
    # 注意：在大规模数据集上，可以考虑使用 Faiss 等向量检索引擎加速
    nn = NearestNeighbors(n_neighbors=1)
    nn.fit(control[[ps_col]])
    
    distances, indices = nn.kneighbors(treated[[ps_col]])
    
    # 关键工程实践： flatten 数组并应用卡尺阈值
    matched_indices = indices.flatten()
    distances_flat = distances.flatten()
    
    # 只有当距离在卡尺范围内时，才认为是有效匹配
    # 0.05 的标准差意味着我们只匹配非常相似的个体
    valid_match_mask = distances_flat <= caliper
    
    valid_treated = treated[valid_match_mask].copy()
    valid_control = control.iloc[matched_indices[valid_match_mask]].copy()
    
    # 为匹配对生成唯一的 PairID，方便后续追踪和调试
    valid_treated['match_id'] = range(len(valid_treated))
    valid_control['match_id'] = range(len(valid_control))
    
    matched_data = pd.concat([valid_treated, valid_control])
    return matched_data, valid_match_mask.sum()

matched_df, n_pairs = perform_matching(data)
print(f"成功匹配了对数: {n_pairs}")

3. 评估与诊断：不要相信你的直觉

在匹配完成后，我们的工作并没有结束。我们必须验证“平衡性”。这是新手经常忽略的一步。

# --- 步骤 3: 评估匹配质量 (SMD) ---
def calculate_smd(df, treatment_col, feature_cols):
    treated = df[df[treatment_col] == 1]
    control = df[df[treatment_col] == 0]
    smds = {}
    for col in feature_cols:
        t_mean = treated[col].mean()
        c_mean = control[col].mean()
        t_std = treated[col].std()
        c_std = control[col].std()
        # Pooled Standard Deviation
        pooled_std = np.sqrt((t_std**2 + c_std**2) / 2)
        smd = abs((t_mean - c_mean) / pooled_std) if pooled_std > 0 else 0
        smds[col] = smd
    return pd.Series(smds)

# 对比匹配前后的 SMD
raw_smd = calculate_smd(data, ‘treatment‘, features)
matched_smd = calculate_smd(matched_df, ‘treatment‘, features)

smd_comparison = pd.DataFrame({‘Raw SMD‘: raw_smd, ‘Matched SMD‘: matched_smd})
print("
标准化均值差异 比较:")
print(smd_comparison)

# 可视化检查 (Love Plot 的简化版逻辑)
# 如果 Matched SMD > 0.1，说明匹配效果不佳，需要重新审视模型或特征

4. 估计最终效果

一旦平衡性达标（通常要求 SMD < 0.1），我们就可以自信地计算平均处理效应（ATE）。

# --- 步骤 4: 计算处理效果 ---
treated_outcome = matched_df[matched_df[‘treatment‘] == 1][‘outcome‘].mean()
control_outcome = matched_df[matched_df[‘treatment‘] == 0][‘outcome‘].mean()

ate = treated_outcome - control_outcome
print(f"
估计的 ATE (Average Treatment Effect): {ate:.4f}")
print(f"真实 ATE 设定约为 3.0 (加上交互效应偏差)")

AI 辅助开发：2026年的最佳实践

作为开发者，我们不仅要写代码，还要善用工具。在实现上述逻辑回归或复杂的 GBM 匹配时，“氛围编程” 是我们的秘密武器。

使用 LLM 驱动的调试

你可能会遇到这样的情况：你的 SMD 始终降不下来。在传统模式下，你可能会花费数小时去查阅统计学文献。但在 2026 年，我们可以这样与 AI 结对编程：

• 上下文注入：将 smd_comparison 的结果直接粘贴给 AI IDE（如 Cursor 或 GitHub Copilot）。
• 询问：“为什么 tenure_months 的 SMD 在匹配后依然很高？请基于我的代码逻辑分析原因。”
• 迭代：AI 可能会指出你混淆了连续变量与分类变量，或者建议你尝试不同的核函数。

这种多模态的开发方式——结合代码、图表和自然语言交互——极大地提高了我们的工作效率。

进阶视角：什么情况下不使用 PSM？

虽然 PSM 很强大，但在我们最近的企业级架构重构中，我们意识到它不是万能药。

• 当样本量极小时：PSM 会丢弃大量无法匹配的数据（尤其是在严格卡尺下），导致统计功效下降。此时，逆概率加权 可能是更好的选择。
• 当存在未观测的混淆因子时：这是 PSM 的阿喀琉斯之踵。如果用户决定接受治疗是因为某种我们数据里没有的情绪（X），那么无论我们如何匹配现有的 X，结果都是有偏的。这时，我们可能需要考虑工具变量法（IV）或者断点回归（RDD）。
• 性能瓶颈：在拥有数百万行数据的实时推荐系统中，进行全局最近邻搜索是非常昂贵的。我们可能会转向近似算法，或者直接使用双鲁棒估计。

总结与展望

从数学定义到 Python 代码实现，倾向得分匹配为我们提供了一把从观测数据中挖掘因果关系的钥匙。在未来的几年里，随着因果机器学习的兴起，PSM 将不再仅仅是一个统计学的步骤，而是融入到 AutoML 和强化学习（RL）的基础设施中。

我们希望这篇文章不仅能帮助你理解 PSM 的原理，更能让你在实际项目中自信地应用它。记住，代码只是工具，对业务逻辑的深刻理解和对偏差的敏锐嗅觉，才是我们作为数据工程师的核心竞争力。

让我们保持好奇心，继续在数据的海洋中探索真相吧！