需求背景：

策略不适用随机分流，在某部分人群全量上线，需要同通过构建相似人群的方式，对策略进行评估。

评估方案：

1、使用PSM构建相似人群，确保实验组与对照组在AA期的评估指标趋势能够保持一致

2、通过DID对实验效果进行评估，确认策略对实验组的影响。

一、构建相似人群

1.1环境包导入


import warnings
warnings.filterwarnings('ignore')
# from pymatch.Matcher import Matcher
import pandas as pd
import numpy as np
%matplotlib inline
from scipy import stats
import matplotlib.pyplot as plt
import patsy
from statsmodels.genmod.generalized_linear_model import GLM
import statsmodels.api as sm
import seaborn as sns
import sys
import os
from jinja2 import Template
import prestodb
from pyhive import hive
from scipy import stats
import patsy
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib
#指定默认字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['font.family']='sans-serif'
#解决负号'-'显示为方块的问题
matplotlib.rcParams['axes.unicode_minus'] = Falsedef presto_read_sql_df(sql):conn=prestodb.dbapi.connect(host=os.environ['PRESTO_HOST'],port=os.environ['PRESTO_PORT'],user=os.environ['JUPYTER_HADOOP_USER'],password=os.environ['HADOOP_USER_PASSWORD'],catalog='hive')cur = conn.cursor()cur.execute(sql)query_result = cur.fetchall() colnames = [part[0] for part in cur.description]raw = pd.DataFrame(query_result, columns=colnames,dtype=np.float64)return rawdef hive_read_sql_df(sql):conn = hive.connect(host=os.environ['PYHIVE_HOST'], port=os.environ['PYHIVE_PORT'],username=os.environ['JUPYTER_HADOOP_USER'],password=os.environ['HADOOP_USER_PASSWORD'],auth='LDAP',configuration={'mapreduce.job.queuename': os.environ['JUPYTER_HADOOP_QUEUE'],'hive.resultset.use.unique.column.names':'false'})raw = pd.read_sql_query(sql,conn)return raw

1.2数据获取


# 1.2.1通过SQL查询数据
# df = hive_read_sql_df("""
#     select
#         XXXX
#     from
#         XXXXX
#     where
#         XXXXX
# """)
# df.to_csv('file_name.csv',encoding='gbk',sep=',',index = False)#1.2.2通过文件读取获得数据
df = pd.read_csv('file_name.csv',encoding='gbk')
y_date = df y_date.dtypes #check 数据的数据类型

1.3模型构建

# 设定匹配因子和分组标签X_field=['XXXX' ,'XXXX','XXXX' #分组标签名称
]
Y_field=['is_group']
field=X_field+Y_field
field

# 循环字段 ：根据需要选择，仅匹配一组，则不需要循环
city_list = {
#                 'XXXX','XXXX','XXXX'}

city_psm_dri = None for x in city_list :city_name = x print('='*20+'城市：'+city_name +'='*20)  query_str = "city_name=='"+city_name+"'"city_info = y_date.query(query_str)  #数据筛选，删除一些非必要的数据city_info = city_info.drop(columns=["city_name"]) #删除城市名称列 city_info = city_info.query("is_strategy_time=='0'") city_info = city_info.query("XXX > 0")#归一化 非必须 city_info.online_days = (city_info.online_days-city_info.online_days.min())/(city_info.online_days.max()-city_info.online_days.min())data=city_infodata=data.dropna()# 区分实验/对照组样本treated=data[(data['is_group']==1)]control=data[(data['is_group']==0)] ttest_result_t=pd.DataFrame(stats.ttest_ind(treated,control))data_p = data# 构建一个线性模型print("-"*20+'模型定义'+"-"*20)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:    '+formula)i=0nmodels=100errors=0model_accuracy = []models = []while i < nmodels and errors < 5:sys.stdout.write('\r{}: {}\{}'.format("Fitting Models on Balanced Samples", i, nmodels)) #第几个模型# sample from majority to create balance datasetdf = control.sample(len(treated)).append(treated).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.columnsb=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] # 模型预测准确性得分
#             print('   model_accuracy:{}'.format(ab_score))model_accuracy.append(ab_score)models.append(res)i += 1except Exception as e:errors += 1 # to avoid infinite loop for misspecified matrixprint('\nError: {}'.format(e))print("\nAverage Accuracy:", "{}%".format(round(np.mean(model_accuracy) * 100, 2)))# 所有模型的平均准确性print('Fitting 1 (Unbalanced) Model...')if errors >= 5: ##异常城市print("【异常警告：】"+city_name+'数据出现异常，该城市数据未记录\n')city_info.to_csv(city_name+'.csv',encoding='gbk',sep=',',index = False)continue glm = GLM(y_f, x_f, family=sm.families.Binomial())res = glm.fit()# model_accuracy.append(self._scores_to_accuracy(res, x_f, y_f))preds = [1.0 if i >= .5 else 0.0 for i in res.predict(x_f)]preds=pd.DataFrame(preds)preds.columns=y_f.columnsb=y_f.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_f)).values[0]model_accuracy.append(ab_score)models.append(res)print("\nAccuracy", round(np.mean(model_accuracy[0]) * 100, 2))scores = np.zeros(len(x_f))for i in range(nmodels):m = models[i]scores += m.predict(x_f)data_p['scores'] = scores/nmodels# 绘图 plt.figure(figsize=(10,5))sns.distplot(data_p[data_p[Y_field[0]]==0].scores, label='Control')sns.distplot(data_p[data_p[Y_field[0]]==1].scores, label='Test')plt.legend(loc='upper right')plt.xlim((0, 1))plt.title("Propensity Scores Before Matching")plt.ylabel("Percentage (%)")plt.xlabel("Scores")threshold=0.00001method='min' # 使用最近样本匹配方法nmatches=1test_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)):#print(i)if i % 10 == 0:sys.stdout.write('\r{}: {}'.format("Fitting Samples", i ))# uf.progress(i+1, len(test_scores), 'Matching Control to Test...')match_id = iscore = test_scores.iloc[i]if method == 'random':bool_match = abs(ctrl_scores - score) <= thresholdmatches = ctrl_scores.loc[bool_match[bool_match.scores].index]elif method == 'min':matches = abs(ctrl_scores - score).sort_values('scores').head(nmatches)else:raise(AssertionError, "Invalid method parameter, use ('random', 'min')")if len(matches) == 0:continue# randomly choose nmatches indices, if len(matches) > nmatchesselect = nmatches if method != 'random' else np.random.choice(range(1, max_rand+1), 1)chosen = np.random.choice(matches.index, min(select, nmatches), replace=False)#     print(chosen)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_idsmatched_data['record_id'] = matched_data.indexprint("\n匹配结果：")print(len(matched_data[matched_data['is_group']==0]['record_id'].unique()))print(len(matched_data[matched_data['is_group']==1]['record_id'].unique()))# m.plot_scores()plt.figure(figsize=(10,5))sns.distplot(matched_data[matched_data[Y_field[0]]==0].scores, label='Control')sns.distplot(matched_data[matched_data[Y_field[0]]==1].scores, label='Test')plt.legend(loc='upper right')plt.xlim((0, 1))plt.title("Propensity Scores After Matching")plt.ylabel("Percentage (%)")plt.xlabel("Scores")matched_data['city_name'] = city_namefile_name = 'result'+city_name+'.csv'matched_data.to_csv(file_name,encoding='gbk',sep=',',index = False)if type(city_psm_dri) == type(None) :city_psm_dri = matched_dataelse :city_psm_dri = pd.concat([city_psm_dri,matched_data])
city_psm_dri.to_csv('result.csv',encoding='gbk',sep=',',index = False)

二、DID效果评估

2.1环境包导入

import statsmodels.formula.api as smf
import pandas as pdimport matplotlib.pyplot as plt
import matplotlib
#指定默认字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['font.family']='sans-serif'
#解决负号'-'显示为方块的问题
matplotlib.rcParams['axes.unicode_minus'] = Falseimport prestodb
import os
from pyhive import hive
def presto_read_sql_df(sql):conn=prestodb.dbapi.connect(host=os.environ['PRESTO_HOST'],port=os.environ['PRESTO_PORT'],user=os.environ['JUPYTER_HADOOP_USER'],password=os.environ['HADOOP_USER_PASSWORD'],catalog='hive')cur = conn.cursor()cur.execute(sql)query_result = cur.fetchall() colnames = [part[0] for part in cur.description]raw = pd.DataFrame(query_result, columns=colnames,dtype=np.float64)return rawdef hive_read_sql_df(sql):conn = hive.connect(host=os.environ['PYHIVE_HOST'], port=os.environ['PYHIVE_PORT'],username=os.environ['JUPYTER_HADOOP_USER'],password=os.environ['HADOOP_USER_PASSWORD'],auth='LDAP',configuration={'mapreduce.job.queuename': os.environ['JUPYTER_HADOOP_QUEUE'],'hive.resultset.use.unique.column.names':'false'})raw = pd.read_sql_query(sql,conn)return raw

2.2数据读取与验证

df = pd.read_csv('20220617171533.csv',encoding='gbk')
print('shape:'+ str(df.shape))
df.dtypes
df.head().T
df.describe()

2.3增量评估

city_list = { 'XXXX'}data = df
city_num=0
for city_name in city_list :city_num=city_num+1"city_name=='"+city_name+"'"city_info = data.query("city_name=='"+city_name+"'")print("="*20+str(city_num)+"="*20)
#     print(city_info.shape)value     = list([city_info.tsh][0])group     = list([city_info.is_group][0])strategy  = list([city_info.is_strategy_time][0])tg        = list([city_info.tg][0])aa = pd.DataFrame({'strategy':strategy,'group':group,'tg':tg,'value':value })X = aa[['strategy', 'group','tg']]y = aa['value']
#     print(X.shape)
#     print(y.shape)est = smf.ols(formula='value ~ strategy + group + tg ', data=aa).fit() y_pred = est.predict(X)aa['value_pred'] = y_pred#策略增量 dat_tsh = est.params.tg / (est.params.strategy + est.params.group  + est.params.Intercept ) p_value = est.pvalues.tg if p_value < 0.05 :print(city_name+":△TSH = %.2f %%，效果显著 " % (dat_tsh*100))else :print(city_name+":△TSH = %.2f %%，效果不显著 " % (dat_tsh*100))#     print(est.params)
#     print(est.pvalues)

2.4趋势绘图

def get_trend_plot(matched_long_new,city_name):yx_trend = matched_long_new.groupby(['is_group','mon']).tsh.mean().reset_index()fig, axes = plt.subplots(1, 1,figsize=(10,4), sharex = False, subplot_kw=dict(frameon=False))x = yx_trend.loc[yx_trend.is_group==1,'mon'].astype(str)y1 = yx_trend.loc[yx_trend.is_group==1,'tsh']y2 = yx_trend.loc[yx_trend.is_group==0,'tsh']axes.plot(x, y1,color='coral', linestyle='-', marker='o', label='实验组')axes.plot(x, y2,color='cornflowerblue', linestyle='--', marker='o', label='对照组')# 取time变量为0，即匹配期的最大日期，作为分界线axes.axvline(x=matched_long_new.loc[matched_long_new.is_strategy_time==0].mon.max(), color='red', linestyle='--')axes.set_xlabel(None)axes.set_ylabel('TSH',rotation=0,labelpad=20)
#     axes.set_xticklabels(labels=[i.replace('2021-','') for i in x],rotation = 40)axes.legend()axes.set_title('%s:XXXXX效果by月'%city_name)plt.show()city_num=0
for city_name in city_list :city_num=city_num+1"city_name=='"+city_name+"'"city_info = data.query("city_name=='"+city_name+"'")get_trend_plot(city_info,city_name)
#     if city_num == 2 :
#         break

若在策略上线前，实验组对照组趋势非常相似，则认为psm匹配效果很好，策略增量置信度高

策略上线前，实验组对照组趋势GAP不稳定，则认为PSM人群不够相似，策略的增量置信度低，建议重新进行人群匹配。

结语：

1、psm选择指标时，需要确保指标对AB两组人群是根据样本自身行为有关的，与AB两组样本是无关的。

2、在构建相似人群时，需要注意人群是否本身存在一定的主观意愿的差异，即用户成为A人群和B人群是否是随机事件。如果存在主观的差异，尽管很多行为指标被psm拉齐，但人群仍然是不同质的，潜在的对策略的影响可能不同，从而导致结果虚高或者虚低。

例如A组是全量天猫商家，B组是全量淘宝商家，在全量B中选择A相似人群。
首先，用户成为A或着B，两组样本本身存在主管意愿的差异，已经具备了样本不同质，对策略的敏感度可能不同。

其次，如果相似指标中包含店铺评分指标，A组和B组虽然都有店铺评分，但A、B的评分体系并不相同，所以这类指标不适宜作为相似指标进行筛选。

（例子不是很合适，我也不能讲我自身的业务，但相信大家可以get到我想说的是什么

本人的相关文章都是相对入门的一些数据分析方法，不会过多的讲解原理，更多的讲解怎么帮助业务解决问题。先在工作中用起来。

如果对你的工作有帮助，帮你节约了一定的时间成本，慷慨的打赏起来吧，一元两元也是爱~

PSM+DID 效果评估python demo 、线性分类模型+双重差分法相关推荐

Bishop 模式识别与机器学习读书笔记 || 线性分类模型之判别函数的几何建模
线性分类模型之判别函数的几何建模文章目录线性分类模型之判别函数的几何建模 1. 判别函数 1.1 两类问题 1.2 多类问题 1.3 Fisher 线性判别 LDA 算法 1.3 代码实现 1.4 ...
逻辑回归和线性回归的区别_[PRML]线性分类模型贝叶斯逻辑回归
线性分类相关文章:1.Fisher线性判别分析(LDA)[1]2.广义模型与线性模型& 判别分析 [2]3.逻辑回归[3]4.线性分类模型简介5.感知机原理及代码复现6.概率生成模型7.概率判 ...
Py之rgf_python：使用Python实现高性能分类模型的完整指南
Py之rgf_python:使用Python实现高性能分类模型的完整指南在机器学习领域,分类问题一直是热门的研究方向.而近年来,随着机器学习技术的普及和算力的提升,越来越多的人开始关注如何构建高性能 ...
【机器学习】SVM支持向量机在手写体数据集上进行二分类、采⽤ hinge loss 和 cross-entropy loss 的线性分类模型分析和对比、网格搜索
2022Fall 机器学习 1. 实验要求考虑两种不同的核函数:i) 线性核函数; ii) ⾼斯核函数可以直接调⽤现成 SVM 软件包来实现⼿动实现采⽤ hinge loss 和 cross-e ...
Python scikit-learn，分类模型的评估，精确率和召回率，classification_report
分类模型的评估标准一般最常见使用的是准确率(estimator.score()),即预测结果正确的百分比. 混淆矩阵: 准确率是相对所有分类结果:精确率.召回率.F1-score是相对于某一个分类的预 ...
基于python的分类模型_python SVM 线性分类模型的实现
运行环境:win10 64位 py 3.6 pycharm 2018.1.1 导入对应的包和数据 import matplotlib.pyplot as plt import numpy as np ...
Python实现Catboost分类模型(CatBoostClassifier算法)项目实战
说明:这是一个机器学习实战项目(附带数据+代码+文档+视频讲解),如需数据+代码+文档+视频讲解可以直接到文章最后获取. 1.项目背景 CatBoost提供最先进的结果,在性能方面与任何领先的机器学习 ...
线性分类模型python_python SVM 线性分类模型的实现
运行环境:win10 64位 py 3.6 pycharm 2018.1.1 导入对应的包和数据 import matplotlib.pyplot as plt import numpy as np ...
python实现决策树分类模型（小白入门超简单实战）
注:由于我不喜欢研究机器学习的原理而更关注于实战,所以本文只讲解python实现决策树模型的代码. 数据集:Iris(鸢尾花卉数据集),是一类多重变量分析的数据集.数据集包含150个数据样本,分为3类 ...
神经网络与深度学习（四）线性分类
目录 3.1 基于Logistic回归的二分类任务 3.1.1 数据集构建 3.1.2 模型构建 3.1.3 损失函数 3.1.4 模型优化 3.1.4.1 梯度计算 3.1.4.2 参数更新 3.1 ...

PSM+DID 效果评估python demo 、线性分类模型+双重差分法