本篇是一个学习的总结：psm背后的原理及其应用，主要目的是梳理自己对知识的理解，也同时分享给和我遇到同样业务问题的同学。

问题背景：

互联网经常会有许多的活动，比如剁手的618，双11，如何衡量这些活动带来的效果呢？

一般我们会用活动带来的增量指标来衡量效果。它是和假设我不做这个活动进行对比的。说到这，自然会想到AB测试，但是也有很多情况下，我们无法使用AB测试来做衡量。

这种场景下，如果还需要衡量增量效果，该如何衡量呢？

解决方案：

使用因果推断的科学方法进行衡量。因果推断中有多种推断方法，其中常用的有DID、PSM等，先从简单的PSM介绍。

使用一个栗子来开始介绍：

假设一个平台，在暑假期间，通过对活跃用户发送红包（需要用户参与并领取）的活动，期望拉动用户付费和活跃。

平台期间有8个用户，4个参与了活动，4个未参与活动，活动前后的数据如下图。现在需要测算下活动带来的效果。

如果我们直接计算参与活动和未参与活动的用户的消费金额进行增量对比，得出的结论是，活动带来的效果是-300。这个数据给到业务，业务人员可能会被气死。。

为什么不能这样衡量？因为忽略了选择性偏差。

选择性偏差的概念：参与活动的用户和不参与活动的用户，他们本身是不同的。

看上图中的最后一列，可以看出来，用户活动之前就有了比较大的区别。也称之为混淆变量。

假如如果我们能够知道A用户在这个活动期间不参加活动的消费金额的话，我就能知道真实的增量了。比如下图中，如果我有第4列的数据，就能计算出来活动增量是200（650-450）

但是我们没有平行时空的数据，无法让我们知道一个用户（A），他参加活动的效果和他同时不参加活动的效果。所以这个假想被称之为反事实框架。

用数学语言来表示则为：

表示第i个人，他同时在实验组或者在控制组，他的结果变量是什么

PSM背后的原理：

虽然我们没有平行时空的数据，但是如果我们能够为实验组用户找到满足以下2个条件的用户，就可以近似得到活动的效果

• 不参与活动的活跃用户：如果我们称参加活动的是实验组，不参加活动的是对照组，则需要从对照组中去找用户

• 和ABCD非常相似的用户：在对照组中，去找和A相似、B相似、C相似、D相似的人群（活动前付费相关特征一致）

相似用户如何衡量和匹配：

如何去衡量两个用户是否相似？倾向性得分psm的叫法也是来源如此。

psm常使用逻辑回归来测算两个用户是否参与实验的概率，将对比用户的类似性从多维度降维到一维数据，再对倾向性得分进行匹配（常用的KNN等方法匹配）。

比如，下表中假设ABCDEFGH有x1,x2,x3 3个特征，然后根据逻辑回归来预测

P（T=1|x1,x2,x3）的概率值。等于原本需要匹配3个维度，通过逻辑回归，降维到匹配一个维度的数值。

        

处理效应（增量如何来衡量以及是否显著）：

上面的例子中实验组的处理效应（增量）为200，如果用数学公式来表示每个个体的平均处理效应，有以下几种表示方法

实际案例：

python现在是有包可以实现计算评分和匹配了，这个实践是了解底层实现原理

使用kaggle上面titanic的数据，了解是否有cabin对最终是否活下去的概率有影响，以及多大的影响。

df = pd.read_csv(r'./kaggle_titanic-master/train.csv')
df.head()

#数据处理
df['treatment'] = df.Cabin.apply(hasCabin)
df_data = df[['treatment','Sex','Age','SibSp','Parch','Embarked', 'Pclass', 'Fare']]
T = df_data.treatment
X = df_data.loc[:,df_data.columns !='treatment']
X_encoded = pd.get_dummies(X, columns = ['Sex','Embarked', 'Pclass'], \
                           prefix = {'Sex':'sex', 'Embarked' : 'embarked', 'Pclass' : 'class'}, drop_first=False)


# 逻辑回归，得到倾向性得分
pipe = Pipeline([
    ('scaler', StandardScaler()),
    ('logistic_classifier', lr())
])
pipe.fit(X_encoded, T)                
predictions = pipe.predict_proba(X_encoded)
 #得到p 
def logit(p):
    logit_value = math.log(p / (1-p))
return logit_value


predictions_logit = np.array([logit(xi) for xi in predictions[:,1]])
df_data.loc[:,'propensity_score'] = predictions[:,1]
df_data.loc[:,'propensity_score_logit'] = predictions_logit
df_data.loc[:,'outcome'] = y.Survived    
X_encoded.loc[:,'propensity_score'] = predictions[:,1]
X_encoded.loc[:,'propensity_score_logit'] = predictions_logit
X_encoded.loc[:,'outcome'] = y.Survived
X_encoded.loc[:,'treatment'] = df_data.treatment
sns.boxplot(x="Pclass", y="propensity_score", data=X_encoded)

从上面这个图中可以看出来，Pclass=1的用户更可能拥有cabin

#KNN matching
caliper = np.std(df_data.propensity_score) * 0.25
print('\nCaliper (radius) is: {:.4f}\n'.format(caliper))
df_data = X_encoded
knn = NearestNeighbors(n_neighbors=10 , p = 2, radius=caliper)
knn.fit(df_data[['propensity_score_logit']].to_numpy())


#得到距离每条记录最近的10条记录
distances , indexes = knn.kneighbors(
df_data[['propensity_score_logit']].to_numpy(), \
n_neighbors=10)


def perfom_matching_v2(row, indexes, df_data):
current_index = int(row['index']) # Obtain value from index-named column, not the actual DF index.
prop_score_logit = row['propensity_score_logit']
for idx in indexes[current_index,:]:
if (current_index != idx) and (row.treatment == 1) and (df_data.loc[idx].treatment == 0):
return int(idx)
        
#得到实验组匹配到在control组的记录
df_data['matched_element'] = df_data.reset_index().apply(perfom_matching_v2, axis = 1, args = (indexes, df_data))

对比匹配前后各个特征的分布情况，从下面可以看出来，匹配之后，实验组和控制组的Age分布相比之前一致很多。

#最终得到处理效应
overview = all_mached_data[['outcome','treatment']].groupby(by = ['treatment']).aggregate([np.mean, np.var, np.std, 'count'])
treated_outcome = overview['outcome']['mean'][1]
treated_counterfactual_outcome = overview['outcome']['mean'][0]
ATT = treated_outcome - treated_counterfactual_outcome
print('The Average Treatment Effect (ATT): {:.4f}'.format(ATT))
#The Average Treatment Effect (ATT): 0.1119

参考资料：

1. B站视频：https://www.bilibili.com/video/BV1gV41117Md?t=681

2. https://github.com/konosp/propensity-score-matching/blob/main/propensity_score_matching_v2.ipynb
   

●Graveyard分析模型是真的牛X！

●品牌知名度分析实例