从‘智商与收入’到‘审计质量’：用Python和Sklearn轻松复现PSM完整流程（附代码）

用Python实现倾向得分匹配从理论到电商场景实战在互联网数据分析领域我们经常需要评估某个干预措施如营销活动、产品改版的真实效果。但简单的对比实验组和对照组均值往往会受到选择偏差的干扰——那些收到优惠券的用户可能本来就是高价值用户。倾向得分匹配PSM正是解决这类问题的有力工具本文将带你用Python完整实现这一计量经济学方法。传统上PSM主要在Stata等统计软件中讨论这让很多熟悉Python的数据科学家望而却步。实际上借助sklearn和statsmodels等库我们完全可以在Python生态中构建更灵活、可扩展的PSM工作流。下面我将通过一个电商优惠券对用户留存影响的案例展示如何用Python实现从倾向得分估计到效应评估的全流程。1. 倾向得分匹配的核心原理倾向得分匹配的本质是为处理组中的每个个体找到双胞胎——在控制组中特征相似但未接受处理的个体。通过比较这些匹配对的结果差异我们能够更准确地估计处理效应。三个关键假设需要满足条件独立假设给定协变量X处理分配T与潜在结果独立共同支撑假设处理组和对照组的倾向得分分布有足够重叠平衡性假设匹配后协变量在处理组和对照组间应达到平衡注意PSM只能解决可观测变量带来的偏差无法处理未观测变量的影响在电商场景中假设我们想评估发送优惠券处理变量对用户7日留存率结果变量的影响。可能影响分配和结果的协变量包括用户历史购买频次最近一次购买金额注册时长设备类型所在城市等级2. 数据准备与探索性分析首先加载必要的Python库并准备模拟数据import pandas as pd import numpy as np from sklearn.ensemble import GradientBoostingClassifier from sklearn.preprocessing import StandardScaler import matplotlib.pyplot as plt import seaborn as sns # 生成模拟数据 np.random.seed(42) n_samples 5000 data pd.DataFrame({ history_freq: np.random.poisson(3, n_samples), last_purchase: np.random.exponential(50, n_samples), member_days: np.random.randint(30, 365, n_samples), device: np.random.choice([iOS,Android], n_samples), city_tier: np.random.choice([1,2,3], n_samples, p[0.2,0.5,0.3]), coupon: np.random.binomial(1, 0.4, n_samples) # 40%用户收到优惠券 }) # 模拟留存率受优惠券和用户特征影响 data[retention] np.where( data[coupon] 1, np.random.binomial(1, 0.3 0.1*data[history_freq]/5 - 0.002*data[last_purchase]), np.random.binomial(1, 0.2 0.08*data[history_freq]/5 - 0.0015*data[last_purchase]) )进行初步的数据探索print(data.groupby(coupon).agg({retention:mean, history_freq:mean, last_purchase:mean}))输出结果可能显示couponretentionhistory_freqlast_purchase00.242.9849.810.323.0548.7初步看收到优惠券的用户留存率更高但这些用户本身购买频次也略高存在选择偏差。3. 倾向得分估计与模型选择倾向得分是给定协变量下个体接受处理的概率通常使用逻辑回归或机器学习模型估计from sklearn.model_selection import train_test_split from sklearn.metrics import roc_auc_score # 准备特征矩阵 X pd.get_dummies(data.drop([coupon,retention], axis1)) y data[coupon] # 划分训练测试集 X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state42) # 使用GBDT估计倾向得分 ps_model GradientBoostingClassifier(n_estimators100, learning_rate0.1, max_depth3) ps_model.fit(X_train, y_train) # 评估模型 train_auc roc_auc_score(y_train, ps_model.predict_proba(X_train)[:,1]) test_auc roc_auc_score(y_test, ps_model.predict_proba(X_test)[:,1]) print(fTrain AUC: {train_auc:.3f}, Test AUC: {test_auc:.3f}) # 为所有样本添加倾向得分 data[pscore] ps_model.predict_proba(X)[:,1]模型选择建议逻辑回归可解释性强但可能欠拟合随机森林/GBDT捕捉非线性关系但可能过拟合正则化模型当特征维度高时推荐使用评估倾向得分模型的质量主要看预测准确性AUC协变量平衡性匹配后检验4. 匹配算法实现与比较Python中有多种实现匹配的方法下面介绍三种主流方法4.1 最近邻匹配k-NNfrom sklearn.neighbors import NearestNeighbors # 分离处理组和对照组 treated data[data[coupon]1] control data[data[coupon]0] # 1:1最近邻匹配 nbrs NearestNeighbors(n_neighbors1).fit(control[pscore].values.reshape(-1,1)) distances, indices nbrs.kneighbors(treated[pscore].values.reshape(-1,1)) # 构建匹配后数据集 matched_control control.iloc[indices.flatten()].copy() matched_data pd.concat([treated, matched_control])4.2 半径匹配Caliper# 设置半径通常取倾向得分标准差的0.2倍 caliper data[pscore].std() * 0.2 matched_pairs [] for i, (idx, row) in enumerate(treated.iterrows()): pscore row[pscore] # 在对照组中寻找倾向得分在(pscore-caliper, pscorecaliper)范围内的样本 candidates control[(control[pscore] pscore - caliper) (control[pscore] pscore caliper)] if len(candidates) 0: # 选择倾向得分最接近的一个 best_match candidates.iloc[(candidates[pscore]-pscore).abs().argsort()[:1]] matched_pairs.append(best_match.index[0]) matched_control_caliper control.loc[matched_pairs].copy()4.3 核匹配Kernelfrom sklearn.neighbors import KernelDensity from scipy.optimize import minimize_scalar # 定义核函数高斯核 def kernel_match(treated_pscore, control_pscore, bandwidth0.1): kde KernelDensity(kernelgaussian, bandwidthbandwidth).fit(control_pscore.reshape(-1,1)) weights np.exp(kde.score_samples(treated_pscore.reshape(-1,1))) return weights / weights.sum() weights kernel_match(treated[pscore].values, control[pscore].values)三种方法对比方法优点缺点适用场景最近邻匹配简单直观匹配质量高可能找不到合适匹配对照组样本充足时半径匹配避免糟糕匹配提高平衡性可能损失大量样本倾向得分分布重叠好时核匹配利用所有样本减少方差计算复杂需要选择合适带宽大样本情况下5. 匹配质量评估匹配后必须检查协变量是否达到了平衡。常用方法包括5.1 标准化偏差计算def calculate_std_diff(treated, control, var): mean_treated treated[var].mean() mean_control control[var].mean() std_treated treated[var].std() std_pooled np.sqrt((std_treated**2 control[var].std()**2)/2) return 100 * (mean_treated - mean_control) / std_pooled # 计算匹配前后各变量的标准化偏差 variables [history_freq, last_purchase, member_days] std_diff_before {var: calculate_std_diff(treated, control, var) for var in variables} std_diff_after {var: calculate_std_diff(treated, matched_control, var) for var in variables} pd.DataFrame({Before Matching: std_diff_before, After Matching: std_diff_after}).T理想情况下匹配后的标准化偏差应小于5%。5.2 倾向得分分布可视化plt.figure(figsize(10,6)) sns.kdeplot(treated[pscore], labelTreated, shadeTrue) sns.kdeplot(matched_control[pscore], labelMatched Control, shadeTrue) plt.title(Propensity Score Distribution After Matching) plt.legend() plt.show()5.3 t检验评估from scipy.stats import ttest_ind for var in variables: t_stat, p_val ttest_ind(treated[var], matched_control[var]) print(f{var}: t-stat{t_stat:.3f}, p-value{p_val:.4f})6. 处理效应估计与结果解释匹配完成后我们可以通过多种方法估计平均处理效应ATE6.1 简单均值比较ate_naive treated[retention].mean() - control[retention].mean() ate_matched treated[retention].mean() - matched_control[retention].mean() print(fNaive ATE: {ate_naive:.3f}) print(fMatched ATE: {ate_matched:.3f})6.2 回归调整import statsmodels.formula.api as smf # 在匹配样本上进行回归调整 model smf.ols(retention ~ coupon history_freq last_purchase member_days, datamatched_data).fit() print(model.summary())6.3 双重稳健估计结合倾向得分和结果模型即使其中一个模型错误设定仍能得到一致估计# 为所有样本预测潜在结果 data[y1] np.where(data[coupon]1, data[retention], model.predict(data.assign(coupon1))) data[y0] np.where(data[coupon]0, data[retention], model.predict(data.assign(coupon0))) # 计算双重稳健估计 dr_ate (data[y1] - data[y0]).mean() print(fDoubly Robust ATE: {dr_ate:.3f})7. 实战建议与常见陷阱在实际项目中应用PSM时有几个关键注意事项样本量要求处理组和对照组原始样本量都应足够大匹配后会损失部分样本确保剩余样本仍有统计功效协变量选择包含所有同时影响处理分配和结果的变量避免包含仅影响结果或仅影响处理的变量考虑加入高阶项和交互项提高倾向得分模型灵活性常见问题排查如果匹配后样本量过少尝试放宽匹配标准或使用核匹配如果平衡性不理想检查倾向得分模型或尝试不同匹配算法对关键结论进行敏感性分析验证结果稳健性与其他方法结合当存在不可观测混杂时考虑结合工具变量法对于面板数据可结合双重差分法(DID)高维数据下可考虑LASSO等变量选择方法在电商优惠券案例中我们最终可能发现原始分析高估了优惠券效果约30%优惠券对低频用户效果更显著在高端机型用户中效果差异更大这些洞察能帮助运营团队更精准地设计优惠策略避免资源浪费。