7.3 倾向得分估计

Rosenbaum & Rubin 的降维奇迹:从高维匹配到单一标量

倾向得分定理:降维的魔法

Rosenbaum & Rubin (1983) 的伟大贡献

定义: 倾向得分(Propensity Score)

直觉: 给定特征,个体接受处理的概率。

Propensity Score Theorem

定理1: 如果

则

惊人结论:

不需要在高维上匹配,只需要在单一标量上匹配!

证明草图:

最后一步用了: 是的函数,给定,的分布是确定的。

定理2 (平衡性):

含义: 给定倾向得分,处理组和对照组的协变量分布相同(平衡)!

核心假设

假设1: 条件独立性(CIA)

强无混淆假设(Strongly Ignorable Treatment Assignment):

白话: 给定观测到的协变量,处理分配与潜在结果独立。

等价表述:

含义: 控制后,没有未观测的混淆变量。

关键: CIA是不可检验的假设!依赖领域知识。

假设2: 共同支撑(Common Support)

含义:

每个个体都有正概率接受或不接受处理
排除"确定性处理"的情况

检查: 倾向得分分布重叠图

python

plt.hist(e[D==1], alpha=0.5, label='处理组', bins=30)
plt.hist(e[D==0], alpha=0.5, label='对照组', bins=30)
plt.xlabel('倾向得分')
plt.ylabel('频数')
plt.legend()

如果违反:

处理组和对照组在某些区域完全不重叠
匹配/加权可能失效

倾向得分匹配(PSM)

算法流程

步骤1: 估计倾向得分

使用逻辑回归(Logit):

或Probit:

Python:

python

from sklearn.linear_model import LogisticRegression

# 估计倾向得分
ps_model = LogisticRegression(penalty='l2', C=1.0)
ps_model.fit(X, D)
e_hat = ps_model.predict_proba(X)[:, 1]

步骤2: 匹配

对于每个处理个体,找到倾向得分最接近的对照个体:

匹配方法

A. 最近邻匹配(Nearest Neighbor):

1:1匹配(每个处理个体配1个对照)
1:匹配(每个处理个体配个对照)

B. 卡尺匹配(Caliper Matching):

其中是卡尺宽度(如)

优点: 避免质量很差的匹配缺点: 部分处理个体可能找不到匹配对象

C. 核匹配(Kernel Matching):

其中是核函数(如Epanechnikov核),是带宽。

优点: 使用所有对照个体的信息(加权) 缺点: 带宽选择影响结果

步骤3: 估计ATT

其中是匹配得到的反事实结果。

Python实现:PSM

python

import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import NearestNeighbors

def estimate_psm(X, D, Y, n_neighbors=1, caliper=None):
    """
    倾向得分匹配(PSM)

    参数:
    ------
    X: (N, K) 协变量
    D: (N,) 处理指示器
    Y: (N,) 结果变量
    n_neighbors: 匹配邻居数
    caliper: 卡尺宽度(None表示无卡尺)

    返回:
    ------
    ATT: 处理组平均处理效应
    matched_data: 匹配后数据
    """
    # 第1步:估计倾向得分
    ps_model = LogisticRegression(max_iter=1000)
    ps_model.fit(X, D)
    e_hat = ps_model.predict_proba(X)[:, 1]

    # 第2步:匹配
    treated_idx = np.where(D == 1)[0]
    control_idx = np.where(D == 0)[0]

    # 用最近邻匹配
    nbrs = NearestNeighbors(n_neighbors=n_neighbors, metric='euclidean')
    nbrs.fit(e_hat[control_idx].reshape(-1, 1))

    distances, indices = nbrs.kneighbors(e_hat[treated_idx].reshape(-1, 1))

    # 应用卡尺
    matched_pairs = []
    for i, (dist, idx) in enumerate(zip(distances, indices)):
        for d, j in zip(dist, idx):
            if caliper is None or d < caliper:
                matched_pairs.append({
                    'treated_idx': treated_idx[i],
                    'control_idx': control_idx[j],
                    'distance': d
                })

    # 第3步:计算ATT
    ATT_list = []
    for pair in matched_pairs:
        i_t = pair['treated_idx']
        i_c = pair['control_idx']
        ATT_list.append(Y[i_t] - Y[i_c])

    ATT = np.mean(ATT_list)

    # 构造匹配后数据
    matched_data = pd.DataFrame({
        'treated_idx': [p['treated_idx'] for p in matched_pairs],
        'control_idx': [p['control_idx'] for p in matched_pairs],
        'Y_treated': [Y[p['treated_idx']] for p in matched_pairs],
        'Y_control': [Y[p['control_idx']] for p in matched_pairs],
        'effect': ATT_list
    })

    return ATT, matched_data, e_hat

# 示例
np.random.seed(42)
N = 1000
X = np.random.randn(N, 5)
true_e = 1 / (1 + np.exp(-X[:, 0] - 0.5*X[:, 1]))
D = (np.random.rand(N) < true_e).astype(int)
Y0 = X[:, 0] + X[:, 1] + np.random.randn(N)
Y1 = Y0 + 2  # 真实ATE=2
Y = D * Y1 + (1-D) * Y0

ATT, matched_data, e_hat = estimate_psm(X, D, Y, n_neighbors=1)
print(f"PSM估计的ATT: {ATT:.3f}")

平衡性检验

为什么重要?

倾向得分定理:

实践: 需要检验匹配/加权后协变量是否平衡!

标准化均值差(SMD)

定义:

其中:

: 处理组第个协变量均值
: 对照组第个协变量均值
: 方差

阈值:

优:
可接受:
差:

Love Plot

可视化: 所有协变量的SMD

python

def balance_table(X, D, w=None, varnames=None):
    """
    计算平衡性统计

    参数:
    ------
    X: (N, K) 协变量
    D: (N,) 处理指示器
    w: (N,) 权重(None表示无加权)
    varnames: 变量名列表

    返回:
    ------
    df: 平衡性统计表
    """
    K = X.shape[1]
    if varnames is None:
        varnames = [f'X{k+1}' for k in range(K)]

    results = []

    for k in range(K):
        x = X[:, k]

        if w is None:
            # 无加权
            mean_t = x[D==1].mean()
            mean_c = x[D==0].mean()
            std_t = x[D==1].std()
            std_c = x[D==0].std()
        else:
            # 加权
            mean_t = np.average(x[D==1], weights=w[D==1])
            mean_c = np.average(x[D==0], weights=w[D==0])
            std_t = np.sqrt(np.average((x[D==1] - mean_t)**2, weights=w[D==1]))
            std_c = np.sqrt(np.average((x[D==0] - mean_c)**2, weights=w[D==0]))

        # SMD
        smd = (mean_t - mean_c) / np.sqrt((std_t**2 + std_c**2) / 2) * 100

        results.append({
            'Variable': varnames[k],
            'Mean_Treated': mean_t,
            'Mean_Control': mean_c,
            'SMD': smd
        })

    return pd.DataFrame(results)

def plot_love(balance_before, balance_after):
    """绘制Love Plot"""
    fig, ax = plt.subplots(figsize=(10, 8))

    y = np.arange(len(balance_before))

    ax.scatter(balance_before['SMD'], y, s=100, marker='o',
               color='red', label='匹配前', alpha=0.7)
    ax.scatter(balance_after['SMD'], y, s=100, marker='s',
               color='blue', label='匹配后', alpha=0.7)

    ax.axvline(x=-10, color='gray', linestyle='--', linewidth=1, alpha=0.5)
    ax.axvline(x=10, color='gray', linestyle='--', linewidth=1, alpha=0.5)
    ax.axvline(x=0, color='black', linestyle='-', linewidth=1)

    ax.set_yticks(y)
    ax.set_yticklabels(balance_before['Variable'])
    ax.set_xlabel('标准化均值差 (%)', fontsize=12, fontweight='bold')
    ax.set_title('Love Plot: 协变量平衡性', fontsize=14, fontweight='bold')
    ax.legend(fontsize=11)
    ax.grid(True, axis='x', alpha=0.3)

    plt.tight_layout()
    return fig

敏感性分析

Rosenbaum Bounds

问题: CIA假设无法检验!

问题: 如果存在未观测混淆,结论如何改变?

Rosenbaum (2002) 的敏感性分析:

设定: 两个匹配个体,倾向得分的几率比:

解释:

: 完美匹配(无隐藏偏差)
: 匹配个体的接受处理几率相差2倍
: 完全的选择性偏差

分析: 对于不同的,计算效应的上下界

Python:

python

def rosenbaum_bounds(matched_data, gamma_range):
    """
    Rosenbaum边界敏感性分析

    参数:
    ------
    matched_data: 匹配数据(包含treated和control结果)
    gamma_range: Gamma值列表

    返回:
    ------
    bounds: 上下界
    """
    effects = matched_data['Y_treated'] - matched_data['Y_control']
    n = len(effects)

    results = []
    for gamma in gamma_range:
        # 计算Wilcoxon统计量的上下界
        # (简化版,完整版需要更复杂的计算)
        p_plus = gamma / (1 + gamma)
        p_minus = 1 / (1 + gamma)

        # 上界(最不利情况)
        upper = np.mean(effects) + np.std(effects) / np.sqrt(n) * 1.96

        # 下界
        lower = np.mean(effects) - np.std(effects) / np.sqrt(n) * 1.96

        results.append({
            'Gamma': gamma,
            'Lower': lower,
            'Upper': upper
        })

    return pd.DataFrame(results)

# 使用
gammas = [1.0, 1.5, 2.0, 2.5, 3.0]
bounds = rosenbaum_bounds(matched_data, gammas)
print(bounds)

本节小结

核心要点

1. 倾向得分定理: 降维奇迹——从高维到单维

2. PSM 算法:

步骤1: 用Logit/Probit估计
步骤2: 在上匹配(最近邻/卡尺/核)
步骤3: 估计ATT

3. 核心假设:

CIA(条件独立性): 不可检验!
共同支撑: 可检验,必须满足

4. 必做检验:

平衡性检验(Love Plot, SMD)
敏感性分析(Rosenbaum Bounds)

关键公式

倾向得分:

PSM-ATT:

上一节: 7.2 匹配方法基础 | 下一节: 7.4 IPW与双重稳健估计

7.3 倾向得分估计 ​

倾向得分定理:降维的魔法 ​

Rosenbaum & Rubin (1983) 的伟大贡献 ​

Propensity Score Theorem ​

核心假设 ​

假设1: 条件独立性(CIA) ​

假设2: 共同支撑(Common Support) ​

倾向得分匹配(PSM) ​

算法流程 ​

匹配方法 ​

Python实现:PSM ​

平衡性检验 ​

为什么重要? ​

标准化均值差(SMD) ​

Love Plot ​

敏感性分析 ​

Rosenbaum Bounds ​

本节小结 ​

核心要点 ​

关键公式 ​

7.3 倾向得分估计

倾向得分定理:降维的魔法

Rosenbaum & Rubin (1983) 的伟大贡献

Propensity Score Theorem

核心假设

假设1: 条件独立性(CIA)

假设2: 共同支撑(Common Support)

倾向得分匹配(PSM)

算法流程

匹配方法

Python实现:PSM

平衡性检验

为什么重要?

标准化均值差(SMD)

Love Plot

敏感性分析

Rosenbaum Bounds

本节小结

核心要点

关键公式