因果推断

2024-06-17 12:48| 来源: 网络整理| 查看: 265

因果推断 | 倾向得分匹配（PSM）// Propensity Score Matching 发表于 2022-03-23 更新于 2023-09-22 分类于 Casual Inference 因果推断阅读次数： Valine：倾向评分匹配

倾向评分匹配/倾向性得分匹配（PSM，Propensity Score Matching）

通过观测变量对实验组的预测概率进行匹配处理观察性研究（observational study）的经典方法可以解决内生性问题

适用情形：

在观察研究中，对照组与实验组中可直接比较的个体数量很少由于衡量个体特征的参数很多，所以想从对照组中选出一个跟实验组在各项参数上都相同或相近的子集作对比变得非常困难对照组和实验组之间的差异可能是选择性偏差导致的，不一定是干预策略导致的

Picture from PROPENSITY SCORE MATCHING（PSM图解）

倾向性得分：一个用户/样本属于实验组的“倾向性”（给定$x$的情况下，样本进入实验组的条件概率）。$$e(x)=P\left(T=1|X=x\right)$$其中，

$T$：干预$X$：样本属性/协变量$Y$：观测结果/结果变量

Theorem 1: 平衡性质（Balancing Property）:$$T_i \perp X_i \mid e(X_i)$$

Theorem 2: 无混淆（Unconfoundedness）：$$ T_i \perp Y_{i0}, Y_{i1} \mid e(X_i)$$

倾向性得分是一种“Balancing score” 对于倾向性得分$e(X_i)$相同的一群用户$i$：干预$T_i$和特征$X_i$是独立的干预$T_i$和潜在结果$Y_{i0}, Y_{i1} $是独立的当协变量$X$不影响干预变量$T$、但影响结果变量$Y$时，引入协变量可以提高干预精度当遗漏了重要的协变量时，会导致偏差前提/假设

PSM需满足以下两个前提/假设：

条件独立假设：控制了协变量后，样本是否接收干预与潜在结果相互独立共同支撑条件：若匹配结果较为理想，匹配后的两条核密度曲线较为相近匹配方法

PSM通过统计模型对每个样本计算倾向性得分，再按照倾向得分是否接近进行匹配。以分组变量为解释变量（分类变量），样本属性（可能会影响到结果的变量）作为被解释变量，进行逻辑回归（Logistic Regression），从而计算每个样本的倾向性得分。

一对一匹配：在控制组中选择与实验组样本最相近的单个样本进行匹配优点：匹配的样本之间较为相近，偏差较小缺点：样本匹配量小、估计的方差较大一对多匹配：在控制组中选择多个与实验组样本相近的样本进行匹配缺点：在控制组中找到的第2个、第3个、……样本与实验组样本的近似性会减弱，会增大估计的偏差

控制组样本量较多时，可以考虑使用一对多匹配

匹配方法特点偏差方差最近邻匹配增加近邻数，可重复使用半径匹配容忍度增加核匹配带宽增加卡尺匹配容忍度增加分块匹配法最近邻匹配

最近邻匹配（Nearest-neighbor Matching）：选择控制组中倾向性得分最接近的样本作为匹配样本

干预样本$i$与无干预样本$j$相匹配：$$\begin{equation}\mid p_i - p_j \mid = \min_{k \in { T=0 } } \left\lbrace \mid p_i - p_k \mid \right\rbrace\end{equation}$$

半径匹配

半径匹配（Radius Matching）：

核匹配

核匹配（Kernel Matching）

卡尺匹配

卡尺匹配（Caliper Matching）：给定误差 $\delta > 0$，干预样本$i$与无干预样本$j$相匹配$$\begin{equation}\delta > \mid p_i - p_j \mid = \min_{k\in { T=0 } } \left\lbrace \mid p_i - p_k \mid \right\rbrace\end{equation}$$

可能会有样本找不到匹配对象分块匹配法

分块匹配法（Stratification）

匹配效果检验

应用因果推断方法尽可能保证干预组与控制组的分布较为相近（尽可能在已观测到的变量上同质），尽可能去“模拟”随机试验的情况，需要结合一定检验方法来检验“模拟”效果，以保障因果推断评估结果的稳健性（Robust）。

共同支撑检验

共同支撑检验（Common Support Test）：为了提高匹配质量，仅保留倾向得分重叠部分的个体（会损失一定的样本量）。

方法：

比较匹配前后的核密度图绘制条形图，显示倾向得分的共同取值范围平衡性检验

平衡性检验（Balance Diagnose）：观察匹配后的各个变量的均值是否有显著差异。计算每个变量的干预组和虚拟控制组的差异（SMD，Standardized Mean Difference）。

$$\begin{equation}\mathrm{ SMD } = \frac{ \bar{ X }{ \mathrm{ Treat } } - \bar{ X }{ \mathrm{ Control } } }{ \sqrt{ \frac{1}{2} \left( s^2_{ \mathrm{ Treat } } + s^2_{ \mathrm{ Control } } \right) } }\end{equation}$$

相关论文：Zhang, Z., Kim, H. J., Lonjon, G., & Zhu, Y. (2019). Balance diagnostics after propensity score matching. Annals of translational medicine, 7(1). 经验表明，如果一个变量的SMD不超过0.1，一般认为这个变量的匹配质量可以接受；否则，需要结合经验来判断该变量是否可以剔除。反驳式检验

反驳式检验（Refute）：

添加随机混杂因子：添加随机的混杂因子（控制特征），如果因果效应变化较小，表明PSM结果较为稳健安慰剂检验：将干预替换为随机变量，如果因果效应接近于0，表明PSM结果较为稳健

反驳式检验本质上是弱检验，检验通过表明PSM估计结果在已有的控制变量上能够保持较好的稳健性；检验不通过，则认为PSM估计结果是不可用的，需要重新评估。

敏感性分析

敏感性分析（Sensitivity Analysis）：TBD

Pros & Cons优势 Strengths 在计算倾向性得分时，考虑了交互项（interaction terms）可用于二分变量（dichotomous variables）和连续变量（continuous variables）较少依赖于p值或其他特定于模型的假设局限性 Limitations 需要大样本（PSA works best in large samples to obtain a good balance of covariates.）要求控制组和实验组有较大的共同取值范围（Group overlap must be substantial to enable appropriate matching.）要尽可能控制可观测的变量，如果存在不可观测的协变量，会导致“隐形偏差”（The most serious limitation is that PSA only controls for measured covariates. Matching on observed covariates may open backdoor paths in unobserved covariates and exacerbate hidden bias.）没有考虑样本的聚集性，邻域层级的研究可能会存在问题（Does not take into account clustering (problematic for neighborhood-level research).） Pythonpymatch GitHub地址：https://github.com/benmiroglio/pymatch 不要使用pip install pymatch安装，应下载修正版至pymatch同名文件夹，直接加载使用 compare_continuous会报错，相应解决方案可参考https://github.com/benmiroglio/pymatch/issues/36 可使用causalinference包计算ATE（参考Causal Inference fro the Brave and True - Propensity Score - Propensity Score matching） causalinference 官方地址：https://causalinferenceinpython.org/ GitHub地址：https://github.com/laurencium/Causalinference 计算ATE：https://matheusfacure.github.io/python-causality-handbook/11-Propensity-Score.html#propensity-score-matching DoWhy 官方地址：https://www.pywhy.org/dowhy

相关资料：

[ZhiHu]因果推断框架 DoWhy入门实例伪Demo（pymatch）123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197# 载入需要的库import pandas as pdimport numpy as npimport osimport warningswarnings.filterwarnings('ignore')from pymatch.Matcher import Matcherfrom causalinference.causal import CausalModelfrom xgboost import XGBClassifier, plot_importancefrom sklearn.model_selection import train_test_splitfrom sklearn.metrics import roc_auc_score, accuracy_score, precision_score, recall_score, f1_scorefrom scipy import statsfrom statsmodels.stats.power import GofChisquarePowerfrom statsmodels.stats.proportion import proportion_effectsize as esimport matplotlib.pyplot as pltfrom IPython.display import display# 自定义检验函数def chi2_test(df1, df2, col): ''' 卡方检验 ''' n_table = [] ## 存储列联表 ## 待检验指标的类别枚举值 vals1 = set(df1[col].unique()) vals2 = set(df2[col].unique()) col_values = set(vals1 + vals2) for val in col_values: if val in vals1: a = len(df1[df1[col] == val]) else: a = 0 if val in vals2: b = len(df2[df2[col] == val]) else: b = 0 n_table.append([a, b]) stat, pval, d_f, expected = stats.chi2_contingency(n_table) return pvaldef continuous_test(df1, df2, col): ''' 连续指标的显著性检验 ''' result = [] if len(df1[col]) < 30 or len(df2[col]) < 30: result = [-1.0, '样本量过少'] else: ## 1. 正态性检验 norm_stat1, norm_pval1 = stats.kstest(df1[col], 'norm') norm_stat2, norm_pval2 = stats.kstest(df2[col], 'norm') if norm_pval1 > 0.05 and norm_pval2 > 0.05: ## 正态性检验通过 ## 2.1 方差齐性检验 F = df1[col].var() / df2[col].var() n1 = len(df1[col]) - 1 n2 = len(df2[col]) - 1 F_pval = 1 - stats.f.cdf(F, n1, n2) if F_pval < 0.05: ## 方差齐性检验不通过 ## 3.1 不等方差T检验（Welch T Test） stat, pval = stats.ttest_ind(df1[col], df2[col], equal_var=False) test_method = 'Welch-t检验' else: ## 方差齐性检验通过 ## 3.2 等方差T检验（Student T Test） stat, pval = stats.ttest_ind(df1[col], df2[col], equal_var=True) test_method = 'Student-t检验' else: ## 2.2 MWU非参数检验 stat, pval = stats.mannwhitneyu(df1[col], df2[col], alternative='two-sided') test_method = 'Mann-Whitney-U检验' result = [pval, test_method] return result# 数据准备df = pd.read_csv("xxx.csv", encoding='utf_8_sig')df.head()len(df)## 一些数据处理## 1. 缺失值填补cols_fn0 = [ ## 一些需要用0填补缺失值的列名]for col in cols_fn0: df[col] = df[col].fillna(0)## 2. 剔除多余的列cols_drop = [ ## 需要剔除的列]df = df.drop(cols_drop, axis=1)# 特征重要性（筛选变量）cols_x = [ ## 特征]cols_y = [ ## 因变量]X = df[cols_x]y = df[cols_y]## 如果特征的量纲较大，应先进行标准化arr_mean = np.mean(X)arr_std = np.std(X, ddof=1)newX = (X - arr_mean) / arr_stdX = newXprint("XGBoost Start!")# 训练模型model = XGBoostClassifier()model.fit(X, y)df_ip = pd.DataFrame(columns=['Feature', 'Importance'])df_ip['Feature'] = X.columnsdf_ip['Importance'] = model.feature_importances_df_ip = df_ip.sort_values(by=['Importance'], ascending=False)df_ip = df_ip.reset_index(dro=p=True)df_ip.to_csv('xxxxx_xgboost_feature_importance.csv', index=False, encoding='utf_8_sig')cols_ip = df_ip.Feature[:10].to_list() ## 选取TOP 10特征最为PSM的控制特征display(df_ip)print("XGBoost Done!\n")cols_cate = [ ## Category特征]for col in cols_cate: if col in cols_ip: df[col] = df[col].astype('category')# PSMprint("PSM Start!")df1 = df[cols_ip + ['new_index', 'group_name']]experiment = df1[df1.group_name == 'experiment'].reset_index(drop=True)control = df1[df1.group_name == 'control'].reset_index(drop=True)experiment['treatment'] = 1control['treatment'] = 0cols_exclude = [ ## 不参与PSM的特征]m = Matcher(experiment, control, yvar='treatment', exclude=cols_exclude)np.random.seed(2022) ## 设置随机seed，以便后续结果能够复现m.fit_scores(balance=True, nmodels=100)print("1. Fit Scores Done!")m.predict_scores()fig = plt.figure(figsize=(10, 5), dpi=200)m.plot_scores()fig.savefig('xxx_PSM_scores_plot.jpg')print("2. Scores Plot Done!")fig = plt.figure(figsize=(10, 5), dpi=200)m.tune_threshold(method='random')fig.savefig('xxx_PSM_tune_threshold_plot.jpg')print("3. Tune Threshold Plot Done!")m.match(method='min', nmatches=1, threshold=0.0001)print("4. Match Done!")df_freq = m.record_frequency()df_freq.to_csv('xxx_PSM_record_frequency.csv', index=False, encoding='utf_8_sig')print("5. Record Frequency Done!")df_matched = m.matched_data.sort_values('match_id') ## 匹配后的数据df_matched.to_csv('xxx_PSM_matched_data.csv', index=False, encoding='utf_8_sig')print("6. Matched Data Saving Done!")df_original = m.data ## 源数据（含scores）df_original.to_csv('xxx_PSM_original_data.csv', index=False, encoding='utf_8_sig')print("7. Original Data Saving Done!")if len(set(cols_cate).intersection(set(cols_ip))) < len(cols_ip): cc_continuous = m.compare_continuous(return_table=True) cc_continuous.to_csv('xxx_PSM_compare_continuous.csv', index=False, encoding='utf_8_sig') print("8.1 Compare Continuous Done!")elif len(set(cols_cate).intersection(set(cols_ip))) > 0: cc_categorical = m.compare_categorical(return_table=True) cc_categorical.to_csv('xxx_PSM_compare_categorical.csv', index=False, encoding='utf_8_sig') print("8.2 Compare Categorical Done!")print("All Complete! ^ _ ^\n") 计算显著性差异所需的最小样本量：Post not found: 数据分析-ABtest「样本量」小节 RMatching 使用Matching包实现PSM

参见：R: Sekhon J S . Multivariate and Propensity Score Matching Software with Automated Balance Optimization: The Matching package for R [J]. Journal of statistical software, 2011, 42(i07).

12345678910## 安装Matching包 / install Matching packageinstall.packages("Matching", dependencies=TRUE)install.packages("rgenoud")## 加载Matching包 / library Matching packagelibrary("Matching")## 加载lalonde数据集 / load lalonde datasetdata("lalonde")attach(lalonde) 匹配123456789Y

【本文地址】

因果推断

因果推断

今日新闻

推荐新闻