用Python完成PSM（倾向性匹配评分）算法详解
PSM是什么？什么时候会用到PSM？
PSM（Propensity Score Matching）倾向性评分匹配，是使用非实验数据或观测数据进行干预效应分析的一类统计方法。
为了研究某项措施或是某个行为对人群的影响（例如吸烟对健康的影响，读北大对收入的影响），或者互联网中某项措施对于用户的影响，最直接有效的评估方法是我们从大量的样本中随机选择对照组（control group）和实验组（treated group），保证这些用户在实验前的一致性（通过aa test），然后进行AB test，之后来评估效果。
而在现实中，由于一些因素，并不能有效的开展上述实验，例如不能强迫可以上北大的人不上北大吧，以及组织实验者来吸烟吧。在这种情况下，我们期望能够引入一种方法来找到有效的对照组，以此来评估该项措施的效果。（具体PSM的假设原理和合理性在此不作过多的讨论）
结合Python将具体介绍PSM的实现方式（本文介绍的Python实现方式基于Python 3.X)
数据预处理
引用包
import warnings warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
%matplotlib inline
import os
from scipy import stats
import matplotlib.pyplot as plt
import patsy
import sys
from statsmodels.genmod.generalized_linear_model import GLM
import statsmodels.api as sm import seaborn as sns
读数据
path=os.getcwd() # 获取当前的路径
file=path + "/data.txt" # 写具体的引用文件名，该文件数据包括所有的X变量和Y变量（0，1）
df=pd.read_csv(file,sep='\t') # 结合具体文件类型，若为txt文件则用read_csv, excel 则用read_excel
确定变量
X_field=['X1','X2','X3','X4','X5','X6']
Y_field=['Y']
field=X_field+Y_field
数据清洗，分辨出实验集和预选对照集
data=df[field]
data=data.dropna() # 去除任何包含无效数据的行
treated=data[data['Y']==1]
control=data[data['Y']==0]
treated[field].describe().round(2) #PSM前对照组和实验组的差异
control[field].describe().round(2)
根据数据量的大小，选择PSM的实验组数据和对照组预选数据
（在此以文件量较大的情况来进行进一步筛选）
data.groupby(Y_field).size() #原数据样本量大小
treated_sample=treated.sample(10000,axis=0)[field].reset_index(drop=True) # 实验组
control_sample=control.sample(100000,axis=0)[field].reset_index(drop=True) #对照组，数据量必须大于实验组数据量
PSM算法，step1，计算逻辑回归方程
data_p = control_sample.append(treated_sample).reset_index(drop=True) #两个筛选的数据集为回归基础数据
y_f,x_f=patsy.dmatrices('{} ~ {}'.format(Y_field[0], '+'.join(X_field)), data=data_p, return_type='dataframe')
formula = '{} ~ {}'.format(Y_field[0], '+'.join(X_field))
print('Formula:\n'+formula)
print('n majority:', len(control_sample))
print('n minority:', len(treated_sample)) #确定回归方程
i=0
nmodels=50 可指定归回模型个数
errors=0
model_accuracy = []
models = []
while i < nmodels and errors < 5:
sys.stdout.write('\r{}: {}\{}'.format("Fitting Models on Balanced Samples", i, nmodels)) #第几个模型
control_sample.sample(len(treated_sample)).append(treated_sample).dropna() #模型选择相同的对照组和控制组样本
y_samp, X_samp = patsy.dmatrices(formula, data=df, return_type='dataframe') #选出模型的自变量和因变量
glm = GLM(y_samp, X_samp, family=sm.families.Binomial()) #逻辑回归模型
try:
res = glm.fit()
preds = [1.0 if i >= .5 else 0.0 for i in res.predict(X_samp)]
preds=pd.DataFrame(preds)
preds.columns=y_samp.columns
b=y_samp.reset_index(drop=True)
a=preds.reset_index(drop=True)
ab_score=((a.sort_index().sort_index(axis=1) == b.sort_index().sort_index(axis=1)).sum() * 1.0 / len(y_samp)).values[0] # 模型预测准确性得分
model_accuracy.append(ab_score)
models.append(res) i += 1
except Exception as e:
errors += 1
print('Error: {}'.format(e))
print("\nAverage Accuracy:", "{}%". format(round(np.mean(model_accuracy) * 100, 2))) # 所有模型的平均准确性
PSM step2 为实验组，即treated_sample，在对照组池，即control_group中找与预测值相似的行作为对照组（采用临近匹配法）
threshold=0.001
method='min'
nmatches=1
test_scores = data_p[data_p[Y_field[0]]==True][['scores']]
ctrl_scores = data_p[data_p[Y_field[0]]==False][['scores']]
result, match_ids = [], []
for i in range(len(test_scores)):
match_id = i
score = test_scores.iloc[i]
matches = abs(ctrl_scores - score).sort_values('scores').head(nmatches)
chosen = np.random.choice(matches.index, nmatches, replace=False) result.extend([test_scores.index[i]] + list(chosen))
match_ids.extend([i] * (len(chosen)+1))
ctrl_scores=ctrl_scores.drop(chosen,axis=0)
matched_data =data_p.loc[result]
matched_data['match_id'] = match_ids
matched_data['record_id'] = matched_data.index
