前言

大家好，我是顾先生，PSM倾向性得分匹配法的Python代码实操终于来啦！
​
对于PSM原理不太熟悉的同学可以看看前一篇文章：PSM倾向得分匹配法【上篇：理论篇】

目前网上PSM实操的相关文章都是R语言、SPSS和STATA实现的，少数Python版本代码不全，可读性有限（有些甚至要钱。。。）

所以我想出一版可读性强、能迅速复用的Python版本PSM，让各位同学既能看懂又能快速上手。

本次Python代码实操主要参考了psmatching的源码，并做了一定的修改，地址在文末参考资料。

这次我把每段代码都做了保姆级的注释，相信每位同学都能理解到位，当然肯定有注释不对的地方，也欢迎后台私信我。

本文的代码和数据集可关注我的公众号“顾先生聊数据”，后台发送“psm”后领取~

数据集介绍

为了更好地进行演示，我现编了一个数据集，该数据集以电商场景为基础，判断给用户发送优惠券PUSH是否会影响到用户。

数据集的类别主要由事实层标签（年龄、性别等）和行为层标签（最近一次购买diff、之前使用优惠券情况等）。

重申一下！数据集是我randbetween现编的，不用太较真具体内容。

完整代码

下面的代码我做了尽可能详细的注释，复用时需要修改的地方我也做了标注，如有不合理的地方欢迎后台私聊我哦~

安装psmatching包。

!pip install psmatching

import psmatching.match as psm
import pytest
import pandas as pd
import numpy as np
from psmatching.utilities import *
import statsmodels.api as sm

path及model设置。

#地址
path = "./data/psm/psm_gxslsj_data.csv"
#model由干预项和其他类别标签组成，形式为"干预项~类别特征+列别特征。。。"
model = "PUSH ~ AGE + SEX + VIP_LEVEL + LASTDAY_BUY_DIFF + PREFER_TYPE + LOGTIME_PREFER + USE_COUPON_BEFORE + ACTIVE_LEVEL"
#想要几个匹配项，如k=3，那一个push=1的用户就会匹配三个push=0的近似用户
k = "3"
m = psm.PSMatch(path, model, k)

获得倾向性匹配得分。

df = pd.read_csv(path)
#将用户ID作为数据的新索引
df = df.set_index("ID")
print("\n计算倾向性匹配得分 ...", end = " ")
#利用逻辑回归框架计算倾向得分，即广义线性估计 + 二项式Binomial
glm_binom = sm.formula.glm(formula = model, data = df, family = sm.families.Binomial())
#拟合拟合给定family的广义线性模型
#https://www.w3cschool.cn/doc_statsmodels/statsmodels-generated-statsmodels-genmod-generalized_linear_model-glm-fit.html?lang=en
result = glm_binom.fit()
# 输出回归分析的摘要
# print(result.summary)
propensity_scores = result.fittedvalues
print("\n计算完成!")
#将倾向性匹配得分写入data
df["PROPENSITY"] = propensity_scores

计算倾向性匹配得分 …
计算完成!

groups是干预项，propensity是倾向性匹配得分，这里要分开干预与非干预，且确保n1<n2。
注意：将PUSH替换成自己的干预项。

groups = df.PUSH
propensity = df.PROPENSITY
#把干预项替换成True和False
groups = groups == groups.unique()[1]
n = len(groups)
#计算True和False的数量
n1 = groups[groups==1].sum()
n2 = n-n1
g1, g2 = propensity[groups==1], propensity[groups==0]
#确保n2>n1，,少的匹配多的，否则交换下
if n1 > n2:
    n1, n2, g1, g2 = n2, n1, g2, g1
随机排序干预组，减少原始排序的影响
m_order = list(np.random.permutation(groups[groups==1].index))

根据倾向评分差异将干预组与对照组进行匹配
注意：caliper = None可以替换成自己想要的精度

matches = {}
k = int(k)
print("\n给每个干预组匹配 [" + str(k) + "] 个对照组 ... ", end = " ")
for m in m_order:
    # 计算所有倾向得分差异,这里用了最粗暴的绝对值
    # 将propensity[groups==1]分别拿出来，每一个都与所有的propensity[groups==0]相减
    dist = abs(g1[m]-g2)
    array = np.array(dist)
    #如果无放回地匹配，最后会出现要选取3个匹配对象，但是只有一个候选对照组的错误，故进行判断
    if k < len(array):
        # 在array里面选择K个最小的数字，并转换成列表
        k_smallest = np.partition(array, k)[:k].tolist()
        # 用卡尺做判断
        caliper = None
        if caliper:
            caliper = float(caliper)
            # 判断k_smallest是否在定义的卡尺范围
            keep_diffs = [i for i in k_smallest if i <= caliper]
            keep_ids = np.array(dist[dist.isin(keep_diffs)].index)
        else:
            # 如果不用标尺判断，那就直接上k_smallest了
            keep_ids = np.array(dist[dist.isin(k_smallest)].index)
        #  如果keep_ids比要匹配的数量多，那随机选择下，如要少，通过补NA配平数量
        if len(keep_ids) > k:
            matches[m] = list(np.random.choice(keep_ids, k, replace=False))
        elif len(keep_ids) < k:
            while len(matches[m]) <= k:
                matches[m].append("NA")
        else:
            matches[m] = keep_ids.tolist()
        # 判断 replace 是否放回
        replace = False
        if not replace:
            g2 = g2.drop(matches[m])
print("\n匹配完成!")

给每个干预组匹配 [3] 个对照组 …
匹配完成!

这里提一下抽取规则：
抽取规则分无放回匹配和有放回匹配两种。如对照组个体不多，可选择有放回匹配，但重复选择对照组的样本进行匹配，会降低最终匹配样本的样本量，估计精度下降。大家可依据自身样本情况修改代码。

再提一下对应关系：
对应关系分一对一和一对多两种。一对一匹配样本少，估计方差较大，且每个匹配都是最近的，故偏差小； 一对多匹配样本多，估计精度高，但由于干预组个体匹配多个，故相似度降低，偏差会增加。本文使用的是一对多，大家可依据自身样本情况修改代码。

将匹配完成的结果合并起来

matches = pd.DataFrame.from_dict(matches, orient="index")
matches = matches.reset_index()
column_names = {}
column_names["index"] = "干预组"
for i in range(k):
    column_names[i] = str("匹配对照组_" + str(i+1))
matches = matches.rename(columns = column_names)

从全部data中提取干预组和匹配对照的数据。
这里直接调用get_matched_data，注意输入的matches是匹配结果，df是全部数据。

matched_data = get_matched_data(matches, df)

将结果数据写入到文档
注意：可以自己更改地址

print("将倾向性匹配得分写入到文档 ...", end = " ")
save_file = path.split(".")[0] + "_倾向性匹配得分.csv"
df.to_csv(save_file, index = False)
print("完成!")
print("将匹配结果写入到文档 ...", end = " ")
save_file = path.split(".")[0] + "_匹配结果.csv"
matches.to_csv(save_file, index = False)
print("完成!")

将倾向性匹配得分写入到文档 … 完成!
将匹配结果写入到文档 … 完成!

对匹配结果进行变量检验
#注意：将PUSH替换成自己的干预项

variables = df.columns.tolist()[0:-2]
results = {}
print("开始评估匹配 ...")
#注意：将PUSH替换成自己的干预项
for var in variables:
        # 制作用于卡方检验的频率计数交叉表
    crosstable = pd.crosstab(df['PUSH'],df[var])
    if len(df[var].unique().tolist()) <= 2:
        # 计算 2x2 表的卡方统计量、df 和 p 值
        p_val = calc_chi2_2x2(crosstable)[1]
    else:
        # 计算 2x2 表的卡方统计量、df 和 p 值
        p_val = calc_chi2_2xC(crosstable)[1]
    results[var] = p_val
    print("\t" + var + '(' + str(p_val) + ')', end = "")
    if p_val < 0.05:
        print(": 未通过")
    else:
        print(": 通过")
if True in [i < 0.05 for i in results.values()]:
    print("\n变量未全部通过匹配")
else:
    print("\n变量全部通过匹配")

开始评估匹配 … AGE(0.9904): 通过
SEX(0.6688): 通过
VIP_LEVEL(0.0089): 未通过
LASTDAY_BUY_DIFF(0.5134): 通过
PREFER_TYPE(0.7107): 通过
LOGTIME_PREFER(0.2442): 通过
USE_COUPON_BEFORE(0.2961): 通过
ACTIVE_LEVEL(0.7934): 通过
变量未全部通过匹配

这里提一句，为啥P值<0.05未通过，P值>0.05反而通过呢？

这是因为这里的"通过"与假设检验里面的"通过"不太一样。
我们做PSM，是为了在对照组中找到与干预组类似的个体。
那单个个体中的各个变量，应该和干预项"PUSH"是相互独立的。

所以上面的检验中，VIP_LEVEL(0.0089)显示未通过。

后记

ok，PSM的python实操就介绍到这里了，相信看到这里的同学肯定已经学会了如何使用PSM了吧~

本文的代码和数据集可关注我的公众号“顾先生聊数据”，后台发送“psm”后领取~
如有不清楚的地方，也欢迎后台私信我一起讨论。

参考资料：
⼩洛学习群——因果推断（第⼀期）
陈强 高级计量经济学及stata应用（第二版）
https://github.com/rickydangc/psmatching
https://mattzheng.blog.csdn.net/
https://www.jianshu.com/p/34dd19ebe475