因果推断笔记——python 倾向性匹配PSM实现示例（三）

比眉伴天荒 2022-09-07 15:45 1843阅读 0赞

[因果推断笔记—— 相关理论：Rubin Potential、Pearl、倾向性得分、与机器学习异同（二）][_Rubin Potential_Pearl]  
[因果推断笔记——因果图建模之微软开源的dowhy（一）][dowhy]

--------------------

### 文章目录 ###

*  1 从matching -> 倾向性得分
 *   *  1.1 精准匹配
     *  1.2 倾向性得分 / 倾向评分匹配 由来
     *  1.3 精准匹配 -> 倾向评分匹配
 *  2 基于倾向性评分法的因果推断
 *   *  2.1 因果效应估计三种方法
     *  2.2 倾向性评分法的均衡性检验
     *  2.3 反驳
 *  3 倾向性得分案例解读一（无代码）：就业与收入的影响
 *   *  3.1 第一步:使用倾向性评分法估计因果效应
     *  3.2 第二步：评估各倾向性评分方法的均衡性
     *  3.3 第三步：反驳
 *  4 倾向性匹配案例二（无代码）：是否有诊所与死亡率的因果关系
 *   *  4.1 PSM Step1：计算Propensity Score
     *  4.2 PSM Step2：Matching
     *  4.3 PSM Step3：实验组 VS 新对照组 评估建立健康诊所对新生儿死亡率的影响
 *  5 倾向性匹配案例（含代码）案例三
 *   *  5.1 倾向性得分计算
     *  5.2 Matching
 *  6 dowhy实现PSM（含代码）案例四

--------------------

# 1 从matching -> 倾向性得分 #

具体从[因果推断笔记—— 相关理论：Rubin Potential、Pearl、倾向性得分、与机器学习异同（二）][_Rubin Potential_Pearl]摘录的相关的章节

随机化试验那部分我们讲到了：

> T⊥(Y(1), Y(0))(⊥表示独立性)

这个公式其实包含了较强的可忽略性（Ignorability）假定，但我们之前说了，这种方式比较“贵”，所以通常我们会希望收集足够多的X，使得：  
![在这里插入图片描述][b91edadb726942eba5d60c93401a0a21.png]

## 1.1 精准匹配 ##

Conditional Independence Assumption (CIA)假定：![在这里插入图片描述][3833e2ef924a4da6b927e8a6e32f0a49.png]  
那么这个方法实际操作起来通常就是称之为Exact matching（精准匹配），也就是对于每一个 T=1的用户，我们从 T=0的分组里找一个变量 一模一样的用户，在这个X的条件下，满足T与Y相互独立，此时直接比较两组用户观察结果 的差异就可以得到结论，但是，“匹配用户的变量 完全相等” 这个要求过于严格，不实用。

## 1.2 倾向性得分 / 倾向评分匹配 由来 ##

由此衍生了用倾向性得分进行匹配的方式**Propensity Score Matching**，简称 PSM：  
在 PSM 方法中，我们首先对每一个用户计算一个倾向性得分（propensity score），定义为 e(x)=Pr(T=1 | X=x) ，这里面的假定由T⊥(Y(1),Y(0)) | X变为了T⊥(Y(1),Y(0)) | e(x)，接着我们根据倾向性得分对于用户进行匹配。

但是这种方式也是有其局限性的，就是在于treatment group和control group的用户需要足够齐全且维度较高，  
实际应用的时候我们会发现本身两个group的用户有较明显的差异，想要找到倾向性得分相近程度达到要求的用户非常难。  
相关的方法非常多：

PSS：用倾向性得分来对用户进行分组，称为 subclassification  
PSW：用倾向性得分来对用户进行加权抽样，称为 Inverse Propensity Score Weighting (IPSW) ，权重为  
![在这里插入图片描述][ea30c3795175407daafa38cf628438f2.png]

## 1.3 精准匹配 -> 倾向评分匹配 ##

之前提到，  
![在这里插入图片描述][3833e2ef924a4da6b927e8a6e32f0a49.png]  
要想估计出ATE且符合CIA假定，就需要严苛的筛选，譬如，想要了解大学生绩点高低与奖学金的比较，需要：`绩点3.7的，发了奖金的人 - 没法奖金的人`  
条件比较严苛的话，可能找不到，那么就放宽成为区间：`绩点 3.5-3.7，发了奖金的人 - 没法奖金的人`

具体可参考：[DID, PSM 及 DID+PSM 有何差异？DID 要假定不可观测效应随时间变化趋势相同？][DID_ PSM _ DID_PSM _DID]

# 2 基于倾向性评分法的因果推断 #

倾向性评分法由Rosenbaum和Rubin于1983年首次提出，是控制混淆变量的常用方法，其基本原理是将多个混淆变量的影响用一个综合的倾向性评分来表示，从而降低了混淆变量的维度。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70]

倾向性评分是给定混淆变量W的条件下，个体接受Treatment的概率估计，即 P(T=1|W)。  
需要以Treatment为因变量，混淆变量Confounders为自变量，建立回归模型(如Logistic回归)来估计每个研究对象接受Treatment的可能性。

> 回归：T~W  
> ![在这里插入图片描述][eceb8fa2e5f847a2a5b7b99ec1588752.png]

## 2.1 因果效应估计三种方法 ##

需要择一选择：

*  **倾向性评分匹配法(Propensity Score Matching，PSM)**  
    PSM将处理组和对照组中倾向性评分接近的样本进行匹配后得到匹配群体，再在匹配群体中计算因果效应。最常用的匹配方法是最近邻匹配法(nearest neighbor matching)，对于每一个处理组的样本，从对照组选取与其倾向评分最接近的所有样本，并从中随机抽取一个或多个作为匹配对象，未匹配上的样本则舍去。
 *  **倾向性评分分层法(Propensity Score Stratification，PSS)**  
    PSS将所有样本按照倾向性评分大小分为若干层(通常分为5-10层)，此时层内组间混淆变量的分布可以认为是均衡的，当层内有足够样本量时，可以直接对单个层进行分析，也可以对各层效应进行加权平均。当两组的倾向性评分分布偏离较大时，可能有的层中只有对照组个体，而有的层只有试验组的个体，这些层不参与评估因果效应。PSS的关键问题是分层数和权重的设定。可通过比较层内组间倾向性评分的均衡性来检验所选定的层数是否合理，权重一般由各层样本占总样本量的比例来确定。
 *  **倾向性评分加权法(Propensity Score Weighting，PSW)**  
    PSW在计算得出倾向性评分的基础上，通过倾向性评分值赋予每个样本一个相应的权重进行加权，使得处理组和对照组中倾向性评分分布一致，从而达到消除混淆变量影响的目的。

## 2.2 倾向性评分法的均衡性检验 ##

倾向性评分法要求匹配后样本的所有混淆变量在处理组和对照组达到均衡，否则后续分析会有偏差，因此需要对匹配之后的样本进行均衡性检验。  
![在这里插入图片描述][4548c94b9b8b4c738725ebf2fc1da769.png]

## 2.3 反驳 ##

反驳(Refute)使用不同的数据干预方式进行检验，以验证倾向性评分法得出的因果效应的有效性。反驳的基本原理是，对原数据进行某种干预之后，对新的数据重新进行因果效应的估计。  
理论上，如果处理变量(Treatment)和结果变量(Outcome)之间确实存在因果效应，那么这种因果关系是不会随着环境或者数据的变化而变化的，即新的因果效应估计值与原估计值相差不大。  
反驳中进行数据干预的方式有：

*  安慰剂数据法：用安慰剂数据(Placebo)代替真实的处理变量，其中Placebo为随机生成的变量或者对原处理变量进行不放回随机抽样产生的变量。
 *  添加随机混淆变量法：增加一个随机生成的混淆变量。
 *  子集数据法：随机删除一部分数据，新的数据为原数据的一个随机子集。

# 3 倾向性得分案例解读一（无代码）：就业与收入的影响 #

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 1]

Lalonde数据集是因果推断领域的经典数据集  
数据集共包含445个观测对象，一个典型的因果推断案例是研究个人是否参加就业培训对1978年实际收入的影响。  
treatment变量，就业培训与否0/1  
混淆变量为age、educ、black、hisp、married、nodeg。

## 3.1 第一步:使用倾向性评分法估计因果效应 ##

![在这里插入图片描述][b7ab339679624bbcb8fc82cc65e551f5.png]

各种倾向性评分法的因果效应估计值在图表7中，由于不同方法的原理不同，估计的因果效应值也不同。其中倾向性评分匹配法(PSM)因果效应估计值为2196.61，即参加职业培训可以使得一个人的收入增加约2196.61美元  
我们计算ATE(Average Treatment Effect)，即在不考虑任何混淆变量的情况下，参加职业培训(treat=1)和不参加职业培训(treat=0)两个群体收入(re78)的平均差异，在不考虑混淆变量下，参加职业培训可以使得一个人的收入增加约1794.34美元。

另外从ATE和几个估计方法的差异来看，ATE 与PSS/PSW差异不大（说明混淆变量影响不大），PSM差异较大，所以可能PSM不太稳定。

## 3.2 第二步：评估各倾向性评分方法的均衡性 ##

图表8展示了各倾向性评分方法中，每个混淆变量的标准化差值stddiff。总体来看，倾向性评分加权法(PSW)中各混淆变量的标准化差值最小(除了hisp)，说明PSW中混淆变量在处理组和对照组间较均衡，其因果效应估计值可能更可靠。

![在这里插入图片描述][cfc8628ef8b04b9bbc7b07f6ae0c8298.png]

## 3.3 第三步：反驳 ##

图表8展示了100次反驳测试中，三种倾向性评分法的每类反驳测试结果的均值。我们将三种倾向性评分法在真实数据下的因果效应估计值放在图表9最右侧进行对比。  
在安慰剂数据法中，由于生成的安慰剂数据(Placebo)替代了真实的处理变量，每个个体接收培训的事实已不存在，因此反驳测试中的因果估计效应大幅下降，接近0，这反过来说明了处理变量对结果变量具有一定因果效应。  
在添加随机混淆变量法和子集数据法中，反驳测试结果的均值在1585.19~1681.75之间。

对比真实数据的因果估计效应值，PSM的反驳测试结果大符下降，说明其估计的因果效应不太可靠；  
PSW的反驳测试结果与真实数据因果效应估计值最接近，说明其因果效应估计值可能更可靠。

![在这里插入图片描述][28b64b19bcc04a799a8f62f8f2081801.png]  
所以，可以需要得到的结论：

*  需要挑选PSM/PSS/PSW中一个合适的方法
 *  然后来看因果是正向还是负向（因果效应估计值的正负）

--------------------

# 4 倾向性匹配案例二（无代码）：是否有诊所与死亡率的因果关系 #

[Propensity Score Matching（PSM）倾向性匹配评分是什么？盘她！][Propensity Score Matching_PSM]

> 比较清晰了给出了PSM匹配步骤的做法，不过可惜没有代码

对比对照组和实验组项目实施前和项目实施后的新生儿死亡率，进行Difference in Difference 的研究，目前我们假设，并没有历史的新生儿死亡率数据  
数据样式：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 2]  
这里把treatment-T（有没有诊所），Outcome - Y（infant mortality 死亡率），都标明了；两个混杂变量，poverty rate 贫穷率 和per capita doctors 人均医生数。

> 目标：为每一个实验组的村庄创建/找到新的对照组： 针对每一个实验组的村庄，找到其特征类似的对照组。

## 4.1 PSM Step1：计算Propensity Score ##

import statsmodels.api as sm
    formula = "T ~ poverty_rate + per_capita_doctors"
    model = sm.Logit.from_formula(formula, data = data)
    re = model.fit()
    X = data[['poverty_rate', 'per_capita_doctors']]
    data['ps'] = re.predict(X)

这里注意，不是`Y ~ f(T,X)`，而是`T~f(Y,X)`

那么最终出来的结果就是，`每个村庄有诊所的可能性。`

## 4.2 PSM Step2：Matching ##

计算出Propensity Score后，在对照组中需要寻找到与实验组行为（贫穷率、人均医生数）相似的村庄，此过程被称为Matching。  
在这里我们采取最简单的临近匹配法，对每一个实验组村庄进行遍历，找到ps值最接近的对照组村庄作为新对照组集合中的元素，即为new\_control\_index。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 3]  
因为我们需要为有诊所的村庄（T = 1） 一一对应找无诊所（T=0）村庄  
所以，以index=0的实验组村庄为例（ps=0.416571），在健康诊所项目启动前，与其贫穷率、人均医生数最为接近的对照组成员为index=5村庄（ps=0.395162）。

> 这里的计算方式就是非常简单的：ps(index=5) - ps(index = 1)是最小

到此为止，每个实验组村庄都找到了其新的对照组归宿~

## 4.3 PSM Step3：实验组 VS 新对照组 评估建立健康诊所对新生儿死亡率的影响 ##

![在这里插入图片描述][b92d3f6b6e9e445493f1c3adacd06971.png]  
新对照组村庄（未建立健康诊所）新生儿死亡率比实验组村庄（建立健康诊所）足足高出7%，从而证明这个NGO组织的健康诊所项目对新生儿死亡率的降低有显著作用。

--------------------

# 5 倾向性匹配案例（含代码）案例三 #

整个代码非常简单，抽取 + 参考  
文档：[psmatching][]  
github:[psmatching][psmatching 1]

笔者练习代码：[causal\_inference\_demo/psm][causal_inference_demo_psm]

基本实现与案例二无差别，那我们就来看一下。

import psmatching.match as psm
    import pytest
    
    
    path = "simMATCH.csv"
    model = "CASE ~ AGE + TOTAL_YRS"
    k = "3"
    
    # m = psm.PSMatch(path, model, k)
    
    # Instantiate PSMatch object
    m = PSMatch(path, model, k)
    
    # Calculate propensity scores and prepare data for matching
    m.prepare_data()
    
    # Perform matching
    m.match(caliper = None, replace = False)
    
    # Evaluate matches via chi-square test
    m.evaluate()

这里的K=3,代表会找出三个候选集，之前案例二中是一个。

数据集的样子：

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 4]

那么此时：case是干预treatment；  
这个公式，`"CASE ~ AGE + TOTAL_YRS"`就是计算倾向性得分的时候会使用到的

## 5.1 倾向性得分计算 ##

此时在代码内部计算倾向性得分的公式为：

glm_binom = sm.formula.glm(formula = model, data = data, family = sm.families.Binomial())
    
    data = m.df

使用的是广义线性估计 + 二项式Binomial，最后得到了propensity

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 5]

## 5.2 Matching ##

m.matches

这里计算距离的方式，也是最粗暴的：

dist = abs(g1[m]-g2)

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 6]  
所以，matching的过程其实是筛选样本的过程，  
这里最终挑选出来的只有1700+样本（总共有15000个左右）。  
比如这里就是，`case_id = 10384`的，匹配到的三个样本编码分别为：`9504 / 11978 / 595`

当然该案例后续没有给出因果估计量，不过大可以根据`m.matched_data`后续自己计算。  
整个案例代码开源的比较简陋。。将就着看了。。

--------------------

# 6 dowhy实现PSM（含代码）案例四 #

未完待续

[_Rubin Potential_Pearl]: https://mattzheng.blog.csdn.net/article/details/119855174
[dowhy]: https://blog.csdn.net/sinat_26917383/article/details/119851996
[b91edadb726942eba5d60c93401a0a21.png]: /images/20220829/faa224ed10d242a6b98fca0aec846a46.png
[3833e2ef924a4da6b927e8a6e32f0a49.png]: /images/20220829/ff3d834addc24f7aae41587a22a143cb.png
[ea30c3795175407daafa38cf628438f2.png]: /images/20220829/12a12f0a1e994508a3a298a80f601238.png
[DID_ PSM _ DID_PSM _DID]: https://www.zhihu.com/question/35528132
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70]: /images/20220829/2bbab14a37ef4659abfbe4c68ba90767.png
[eceb8fa2e5f847a2a5b7b99ec1588752.png]: /images/20220829/264ebb5443be46b891c1af2cc7073d11.png
[4548c94b9b8b4c738725ebf2fc1da769.png]: /images/20220829/4b9b48535b7d44779b517b219d30f9d1.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 1]: /images/20220829/1913f54620e54d0cb55baad5d88c3990.png
[b7ab339679624bbcb8fc82cc65e551f5.png]: /images/20220829/a8346202fa0c43209ec011743455fb7b.png
[cfc8628ef8b04b9bbc7b07f6ae0c8298.png]: /images/20220829/769b5a79291b4731add19b4d65d5b898.png
[28b64b19bcc04a799a8f62f8f2081801.png]: /images/20220829/0a7ce3eca7aa4d36bd121a6bf20163c2.png
[Propensity Score Matching_PSM]: https://zhuanlan.zhihu.com/p/184958085
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 2]: /images/20220829/511fbb2439d84d12b05cf0ae6f7aee2b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 3]: /images/20220829/cfab5ae93386439096554ed9cc703e4d.png
[b92d3f6b6e9e445493f1c3adacd06971.png]: /images/20220829/34d6d68161364459bd9b9c620b22355f.png
[psmatching]: https://psmatching.readthedocs.io/en/latest/#
[psmatching 1]: https://github.com/rlirey/psmatching
[causal_inference_demo_psm]: https://github.com/mattzheng/causal_inference_demo/tree/main/psm
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 4]: /images/20220829/b3d17f188dcd458683cf4fbea68ee641.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 5]: /images/20220829/e427c2c310a543679c1900e8b005ae52.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3NpbmF0XzI2OTE3Mzgz_size_16_color_FFFFFF_t_70 6]: /images/20220829/49de6872e01248abb5fa27595b3f7e26.png