# Python实战：从零构建倾向得分匹配系统（附完整代码与工业级避坑指南） 在数据科学和策略评估领域，我们常常需要回答"如果"类问题：如果用户接受了某项干预（如营销活动、产品改版），他们的行为会发生怎样的变化？当随机实验不可行时，倾向得分匹配（Propensity Score Matching, PSM）成为了观察性研究中模拟随机对照试验的黄金标准。本文将手把手带你用Python实现工业级的PSM系统，解决实际业务中的因果推断难题。 ## 1. 核心概念与业务场景 在电商平台的优惠券发放场景中，运营团队面临一个经典难题：领取优惠券的用户本身消费意愿就更强，如何剥离用户自身特性影响，准确评估优惠券的真实效果？PSM通过构建"平行宇宙"解决了这个问题——为每个领券用户找到未领券的"双胞胎"，确保比较是在相似用户间进行。 **倾向得分**的本质是用户受到干预的概率估计，数学表示为： ``` P(X) = P(T=1|X) ``` 其中X是用户特征，T=1表示受到干预。通过逻辑回归、梯度提升树等模型，我们可以计算每个用户的倾向得分，然后在得分相近的用户间进行匹配。这种方法依赖于两个关键假设： 1. **条件独立假设**：给定X，潜在结果与干预分配独立 2. **共同支撑假设**：干预组和对照组的倾向得分分布存在重叠区域 > 工业实践中常见的误区是将所有特征盲目纳入倾向得分模型。实际上，应该排除那些可能被干预影响的变量（如干预后的用户行为），否则会引入新的偏差。 ## 2. 数据准备与特征工程 我们模拟一个教育干预场景，评估参加培训课程对学生成绩的影响。首先生成合成数据： ```python import numpy as np import pandas as pd from sklearn.preprocessing import StandardScaler np.random.seed(42) # 生成10000个样本 n_samples = 10000 data = pd.DataFrame({ 'family_income': np.random.normal(50, 15, n_samples), 'parent_education': np.random.choice([1,2,3,4], size=n_samples, p=[0.2,0.3,0.3,0.2]), 'baseline_score': np.random.normal(70, 10, n_samples), 'extracurricular': np.random.poisson(2, n_samples) }) # 生成倾向得分（真实PS模型） true_ps = 1/(1+np.exp(-(0.3*data['family_income']/10 + 0.5*data['parent_education'] - 0.2*data['baseline_score']/10 + 0.4*data['extracurricular']))) data['treatment'] = np.random.binomial(1, true_ps) data['final_score'] = (data['baseline_score'] + 2*data['treatment'] + 0.5*data['family_income']/10 + np.random.normal(0, 5, n_samples)) # 标准化连续变量 scaler = StandardScaler() cont_vars = ['family_income', 'baseline_score'] data[cont_vars] = scaler.fit_transform(data[cont_vars]) ``` 特征工程的关键步骤： 1. **处理缺失值**：对于小于5%缺失的特征采用中位数填充，大于5%考虑添加缺失指示变量 2. **特征变换**：对偏态分布变量进行对数变换 3. **交互特征**：创建父母教育与家庭收入的交互项 4. **平衡检查**：确保干预组和对照组在关键特征上分布均衡 ## 3. 倾向得分建模实战 我们对比三种主流建模方法： | 方法 | 优点 | 缺点 | 适用场景 | |------|------|------|---------| | 逻辑回归 | 可解释性强，计算高效 | 难以捕捉非线性 | 特征间交互较少时 | | 随机森林 | 自动处理非线性关系 | 可能过拟合 | 高维特征空间 | | XGBoost | 处理缺失值，表现优异 | 需要调参 | 大规模数据集 | ```python from sklearn.linear_model import LogisticRegression from sklearn.ensemble import RandomForestClassifier from xgboost import XGBClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import roc_auc_score features = ['family_income', 'parent_education', 'baseline_score', 'extracurricular'] X = data[features] y = data['treatment'] # 划分训练/测试集 X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42) # 初始化模型 models = { "Logistic": LogisticRegression(C=1e6, max_iter=1000), "RandomForest": RandomForestClassifier(n_estimators=100, max_depth=3), "XGBoost": XGBClassifier(n_estimators=100, max_depth=3) } # 训练并评估 results = {} for name, model in models.items(): model.fit(X_train, y_train) pred = model.predict_proba(X_test)[:,1] auc = roc_auc_score(y_test, pred) results[name] = auc data[f'ps_{name}'] = model.predict_proba(X)[:,1] # 存储全量PS print("模型AUC比较:", results) ``` > 模型选择时不应盲目追求AUC最高。在实际业务中，我们更关注匹配质量而非PS预测精度。有时简单模型的匹配效果反而更好。 ## 4. 匹配算法实现与优化 匹配是PSM的核心环节，常见方法对比如下： ```python from sklearn.neighbors import NearestNeighbors def psm_match(data, ps_col, treatment_col, ratio=1, caliper=None): """ 实现带卡尺的最近邻匹配 :param data: 包含PS和干预状态的数据框 :param ps_col: PS列名 :param treatment_col: 干预列名 :param ratio: 匹配比例 (1:1, 1:2等) :param caliper: 卡尺宽度 (PS标准差倍数) :return: 匹配后的数据框 """ treated = data[data[treatment_col]==1].copy() control = data[data[treatment_col]==0].copy() # 计算卡尺 if caliper is not None: ps_std = data[ps_col].std() caliper_val = caliper * ps_std else: caliper_val = None # 最近邻匹配 nbrs = NearestNeighbors(n_neighbors=ratio, metric='euclidean') nbrs.fit(control[[ps_col]].values) matched_control = [] for idx, row in treated.iterrows(): ps = row[ps_col] distances, indices = nbrs.kneighbors([[ps]], n_neighbors=ratio) # 应用卡尺过滤 if caliper_val is not None: mask = distances.flatten() <= caliper_val indices = indices.flatten()[mask] else: indices = indices.flatten() for i in indices: matched_control.append(control.iloc[i].to_dict()) # 合并匹配结果 matched_control_df = pd.DataFrame(matched_control) matched_all = pd.concat([treated, matched_control_df]) matched_all['matched_group'] = np.where(matched_all[treatment_col]==1, 'treated', 'control') return matched_all # 执行1:1卡尺匹配 matched_data = psm_match(data, 'ps_Logistic', 'treatment', ratio=1, caliper=0.2) ``` 匹配算法优化技巧： 1. **分层匹配**：先在关键特征（如性别）上精确匹配，再在子群体内进行PS匹配 2. **核匹配**：使用核函数加权多个邻近样本，减少方差 3. **全局优化**：采用匈牙利算法最小化整体匹配距离 4. **替换策略**：有放回匹配提高质量，无放回匹配保持样本独立性 ## 5. 匹配质量评估体系 匹配后必须系统评估质量，我们构建多维评估指标： ```python def evaluate_balance(data, covariates, treatment_col, before_data=None): """ 评估匹配前后协变量平衡性 :return: 平衡性报告DataFrame """ report = [] for var in covariates: if before_data is not None: # 匹配前 t_before = before_data[before_data[treatment_col]==1][var] c_before = before_data[before_data[treatment_col]==0][var] mean_diff_before = t_before.mean() - c_before.mean() std_before = np.sqrt(t_before.var() + c_before.var()) smd_before = mean_diff_before / std_before # 匹配后 t_after = data[data['matched_group']=='treated'][var] c_after = data[data['matched_group']=='control'][var] mean_diff_after = t_after.mean() - c_after.mean() std_after = np.sqrt(t_after.var() + c_after.var()) smd_after = mean_diff_after / std_after if before_data is not None: report.append({ 'variable': var, 'smd_before': smd_before, 'smd_after': smd_after, 'reduction': 1 - abs(smd_after)/abs(smd_before) if smd_before!=0 else np.nan }) else: report.append({ 'variable': var, 'smd': smd_after }) return pd.DataFrame(report) # 执行评估 balance_report = evaluate_balance(matched_data, features, 'treatment', before_data=data) print("平衡性改善报告:") print(balance_report.sort_values('smd_after', key=abs)) ``` 可视化诊断工具： ```python import matplotlib.pyplot as plt import seaborn as sns def plot_balance_diagnostics(balance_report): plt.figure(figsize=(10, 6)) sns.scatterplot(data=balance_report, x='smd_before', y='smd_after', hue='variable') plt.axhline(0, color='grey', linestyle='--') plt.axvline(0, color='grey', linestyle='--') plt.axhspan(-0.1, 0.1, alpha=0.1, color='green') plt.xlabel('Standardized Mean Difference (Before)') plt.ylabel('Standardized Mean Difference (After)') plt.title('Balance Diagnostics Plot') plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') plt.show() plot_balance_diagnostics(balance_report) ``` ## 6. 处理效应估计与敏感性分析 获得平衡样本后，我们通过双重差分法计算处理效应： ```python def estimate_effect(matched_data, outcome_var): """ 使用匹配后样本估计处理效应 :return: 效应大小及置信区间 """ treated = matched_data[matched_data['matched_group']=='treated'] control = matched_data[matched_data['matched_group']=='control'] effect = treated[outcome_var].mean() - control[outcome_var].mean() # 计算标准误 n_treated = len(treated) n_control = len(control) var_treated = treated[outcome_var].var(ddof=1) var_control = control[outcome_var].var(ddof=1) se = np.sqrt(var_treated/n_treated + var_control/n_control) # 95%置信区间 ci_lower = effect - 1.96 * se ci_upper = effect + 1.96 * se return { 'effect': effect, 'se': se, 'ci_lower': ci_lower, 'ci_upper': ci_upper, 'n_treated': n_treated, 'n_control': n_control } effect_result = estimate_effect(matched_data, 'final_score') print(f"处理效应: {effect_result['effect']:.2f} (95% CI: [{effect_result['ci_lower']:.2f}, {effect_result['ci_upper']:.2f}])") ``` 敏感性分析检查结果稳健性： ```python def sensitivity_analysis(effect, se, gamma_range=np.arange(1, 3.1, 0.2)): """ Rosenbaum敏感性分析 :param gamma: 隐藏偏差的强度 :return: 不同gamma下的p值变化 """ results = [] for gamma in gamma_range: sigma = gamma / (1 + gamma) z = effect / se p_upper = 1 - stats.norm.cdf(z / np.sqrt(sigma)) p_lower = 1 - stats.norm.cdf(z * np.sqrt(sigma)) results.append({'gamma': gamma, 'p_upper': p_upper, 'p_lower': p_lower}) return pd.DataFrame(results) sensitivity_results = sensitivity_analysis(effect_result['effect'], effect_result['se']) print("敏感性分析结果:") print(sensitivity_results[sensitivity_results['gamma'].isin([1, 1.5, 2, 2.5, 3])]) ``` ## 7. 工业级PSM管道实现 将上述模块整合为可复用的PSM管道： ```python class PSMipeline: def __init__(self, treatment_var, outcome_var, covariates, ps_method='logistic', matching_method='nearest', caliper=0.2, ratio=1): """ 初始化PSM管道 :param treatment_var: 干预变量名 :param outcome_var: 结果变量名 :param covariates: 协变量列表 :param ps_method: PS模型类型 (logistic/rf/xgb) :param matching_method: 匹配方法 (nearest/optimal/caliper) :param caliper: 卡尺宽度 (None表示无卡尺) :param ratio: 匹配比例 """ self.treatment_var = treatment_var self.outcome_var = outcome_var self.covariates = covariates self.ps_method = ps_method self.matching_method = matching_method self.caliper = caliper self.ratio = ratio self.ps_model = None self.matched_data = None def fit(self, data): """拟合PS模型并执行匹配""" # 拟合PS模型 if self.ps_method == 'logistic': self.ps_model = LogisticRegression(C=1e6, max_iter=1000) elif self.ps_method == 'rf': self.ps_model = RandomForestClassifier(n_estimators=100, max_depth=3) elif self.ps_method == 'xgb': self.ps_model = XGBClassifier(n_estimators=100, max_depth=3) self.ps_model.fit(data[self.covariates], data[self.treatment_var]) data['ps'] = self.ps_model.predict_proba(data[self.covariates])[:,1] # 执行匹配 self.matched_data = psm_match( data, 'ps', self.treatment_var, ratio=self.ratio, caliper=self.caliper ) return self def check_balance(self, plot=True): """检查匹配质量""" balance_report = evaluate_balance( self.matched_data, self.covariates, self.treatment_var, before_data=data ) if plot: plot_balance_diagnostics(balance_report) return balance_report def estimate_effect(self, outcome_model='simple'): """估计处理效应""" return estimate_effect(self.matched_data, self.outcome_var) def run_full_analysis(self, data): """运行完整分析流程""" self.fit(data) balance_report = self.check_balance() effect_result = self.estimate_effect() sensitivity_results = sensitivity_analysis( effect_result['effect'], effect_result['se'] ) return { 'balance_report': balance_report, 'effect_result': effect_result, 'sensitivity_results': sensitivity_results, 'matched_data': self.matched_data } # 使用示例 ps_pipeline = PSMipeline( treatment_var='treatment', outcome_var='final_score', covariates=features, ps_method='logistic', matching_method='nearest', caliper=0.2, ratio=1 ) results = ps_pipeline.run_full_analysis(data) ``` ## 8. 实战中的陷阱与解决方案 **陷阱1：模型过拟合导致糟糕匹配** - *现象*：PS模型AUC很高但匹配后平衡性差 - *解决*：简化模型，使用正则化或降低树模型深度 **陷阱2：共同支撑区域不足** - *现象*：匹配后样本量大幅减少 - *解决*：放宽卡尺限制或考虑核匹配方法 **陷阱3：隐藏偏差影响结果** - *现象*：敏感性分析显示结果对gamma敏感 - *解决*：收集更多协变量或考虑工具变量法 **陷阱4：时间趋势干扰** - *现象*：干预前后存在外部因素变化 - *解决*：结合双重差分法（PSM-DID） ```python # PSM-DID实现示例 def psm_did_effect(matched_data, outcome_var, time_var, treatment_var): """ 计算PSM-DID处理效应 """ # 预处理期 pre_treated = matched_data[(matched_data['matched_group']=='treated') & (matched_data[time_var]==0)] pre_control = matched_data[(matched_data['matched_group']=='control') & (matched_data[time_var]==0)] # 后处理期 post_treated = matched_data[(matched_data['matched_group']=='treated') & (matched_data[time_var]==1)] post_control = matched_data[(matched_data['matched_group']=='control') & (matched_data[time_var]==1)] # DID计算 effect = (post_treated[outcome_var].mean() - pre_treated[outcome_var].mean()) - \ (post_control[outcome_var].mean() - pre_control[outcome_var].mean()) return effect ``` 在真实业务场景中，我们发现PSM最适用于以下情况： - 干预组和对照组存在明显选择偏差 - 关键混淆变量可观测且可测量 - 样本量足够大（干预组至少数百样本） 对于小样本场景，考虑使用贝叶斯倾向得分方法；当存在不可观测混淆时，工具变量或断点回归可能是更好的选择。