因果推断与反事实预测

2024-07-17 15:05:20| 来源: 网络整理| 查看: 265

文章目录1 导言1.1 价格需求弹性介绍1.2 由盒马反事实预测论文开始1.3 DML - 价格弹性预测推理步骤2 案例详解2.1 数据清理2.2 [v1版]求解价格弹性：OLS回归2.3 [v2版]求解价格弹性：Poisson回归+多元岭回归2.4 [v3版]求解价格弹性：DML2.4.1 DML数据准备 + 建模 + 求残差2.4.2 三块模型对比2.4.3 稳健性评估1 导言1.1 价格需求弹性介绍

经济学课程里谈到价格需求弹性，描述需求数量随商品价格的变动而变化的弹性。价格一般不直接影响需求，而是被用户决策相关的中间变量所中介作用。假设 Q 为某个商品的需求的数量，P 为该商品的价格，则计算需求的价格弹性为，

通过上式可以简单知道，价格改变 1 元比价格改变 100 元，会导致更大的需求改变。比如以 5 元的价格每日可以卖 100 单位产品，如果价格需求弹性为 -3 ，那供应商将价格提升 5%（dp /P，从 5 元-> 5.25 元），需求将下降 15%（dQ/Q ，从 100->85）。那么收入将减少 100_5-5.25_85=53.75。

如果单价降低 5%，那么收入同理将提升 46.25。如果供应商知道了产品的价格弹性，那无须反复测试，即可清楚为提升收入到底应该是提价还是降价。

1.2 由盒马反事实预测论文开始

之前在因果推断笔记——DML ：Double Machine Learning案例学习（十六）这篇提到用DML求价格弹性，不过没有实操模块，本篇是在看过因果推断与反事实预测——盒马KDD2021的一篇论文（二十三）盒马论文之后，想实操一下价格弹性这块。

先来提一下盒马这篇，在反事实预测任务上（随着折扣改变销量如何改变）的尝试半参数模型、XGBtree模型、DeepIV：

第一种，半参数模型，不过这篇对动态折扣下销量的预估的半参数笔者还没深入了解，感觉用分层的价格弹性（平均折扣tree销量预测 + 价格弹性拟合动态折扣销量增量）来规避了核心因果推理的问题，后续要再理解一下该模型第二种，错误尝试，将折扣当作treatment，动态将treatment作为特征来预测销量第三种，deepIV，将三级品类的平均价格（treatment）作为工具变量

三者效果如图，还是semi-para好多了：

本篇是想放大价格弹性的因果计算模块，与盒马的不同：

推估弹性的方法不同（本篇是用DML预测）粒度不同，本篇案例可没顾得上商品分类，一股脑子全放一起了，盒马那篇弹性系数是By 每个商品1.3 DML - 价格弹性预测推理步骤

最好的方式，当然是直接进行 A/B 实验测试不同价格对用户的需求反应，但是价格这类的外生因素在同一产品同一阶段上，对不同用户展示不同的价格会直接损坏用户体验。因此从观察历史数据进行因果推断，但混杂因素（季节性、产品质量等）如何控制是因果推断的挑战。

这里采用 DML（Double Machine Learning）方法进行因果推断，该方法主要解决两个问题：

第一，通过正则化挑拣重要控制变量；第二，对比传统的线性回归模型，用非参数推断可以解决非线性问题。

DML 先应用机器学习算法去分别通过特征变量 X, W 拟合结果变量 Y 和处理变量 T，然后通过线性模型，使用处理变量的残差拟合出结果变量的残差。

目标是估计，这里的 Y 函数构成为 T 的因果作用和 X、W 的协同作用之和。

本篇整个价格弹性的推理过程：

将数据分为两部分，一部分样本选用随机森林等模型，用混杂变量预测处理变量（价格 P），得到 EP|X；另外的样本同样可选择随机森林模型，用混杂变量预测结果变量（需求量 Q），得到 EQ|X。计算残差，得到不受混杂变量影响的价格 P 和需求量 Q，即为 \widetilde P , \widetilde Q\widetilde P = P-E[P|X]\widetilde Q = Q-E[Q|X]因此直接将\widetilde P , \widetilde Q 进行 log-log 回归就能得到弹性系数 𝜃:

需要得到

\theta=\frac{d \widetilde Q / \widetilde Q}{d \widetilde P / \widetilde P}

倒推用log-log回归得到回归系数，即

log \widetilde Q ～ \theta * log \widetilde P + 截距2 案例详解

与本节关联的文章：

因果推断与反事实预测——盒马KDD2021的一篇论文（二十三）因果推断笔记——DML ：Double Machine Learning案例学习（十六）

下面的案例的来源：

利用机器学习因果推理进行弹性定价数据分析36计(29)：价格需求弹性和因果推断简单版代码：DML.ipynb数据集来源：Association Rules and Market Basket Analysis2.1 数据清理

数据集是kaggle的比赛数据集，原文ipynb直接读入的时候会格式报错，这里贴一段kaggle原生读入的方式，不会有报错：

代码语言：javascript复制data = pd.read_csv('OnlineRetail.csv',encoding= 'cp1252',parse_dates=['InvoiceDate']) data = data.sort_values(by='InvoiceDate') data = data.set_index('InvoiceDate')

原数据是购物篮分析数据，这个数据集包含了一家英国在线零售公司在8个月期间的所有购买行为。

每个商品，在每个国家，每家店，每个时间出售的件数与对应的单价。

这里需要额外加工收入：

代码语言：javascript复制df['revenue'] = df.Quantity * df.UnitPrice

同时对P / Q进行对数化处理：

代码语言：javascript复制# 将单价和数量取log df_mdl = df_mdl.assign( LnP = np.log(df_mdl['UnitPrice']), LnQ = np.log(df_mdl['Quantity']), )

2.2 v1版求解价格弹性：OLS回归

v1版 = LnQ~LnP，没有协变量，用最简单的OLS回归

最简单的求解，也不管啥因果推断，有偏无偏，将上述数据的lnp和lnQ，一股脑子都分段，比如(-2.814,-0.868)就是这区间内lnp和lnQ的平均值，如下：

新生成的LnP和LnQ直接回归即得回归系数：

代码语言：javascript复制x='LnP' y='LnQ' df = df_mdl n_bins=15 x_bin = x + '_bin' df[x_bin] = pd.qcut(df[x], n_bins) tmp = df.groupby(x_bin).agg({ x: 'mean', y: 'mean' }) # 回归 mdl = sm.OLS(tmp[y], sm.add_constant(tmp[x])) res = mdl.fit()

得到结果：

弹性系数为-0.6064，价格越高，销量越少

v1的计算也可以使用另外一种方式，计算方差，

因为只有两个变量可以：

代码语言：javascript复制df_mdl[['LnP', 'LnQ']].cov()

这里就是:

\theta=\frac{-0.52}{0.9}=-0.602.3 v2版求解价格弹性：Poisson回归+多元岭回归

v2版 = LnQ~LnP+Country+StockCode+Date，有多元协变量，用岭回归+泊松回归

代码语言：javascript复制import sklearn.preprocessing from sklearn import linear_model from sklearn.pipeline import Pipeline from sklearn.compose import ColumnTransformer from sklearn.preprocessing import OneHotEncoder, StandardScaler, RobustScaler from sklearn.feature_extraction.text import CountVectorizer feature_generator_basic = ColumnTransformer( [ ('StockCode', OneHotEncoder(), ['StockCode']), ('Date', OneHotEncoder(), ['Date']), ('Country', OneHotEncoder(), ['Country']), ('LnP', 'passthrough', ['LnP']), ], remainder='drop' ) mdl_basic = Pipeline([ ('feat_proc', feature_generator_basic), ('reg', linear_model.PoissonRegressor( alpha=1e-6, # l2 penalty strength; manually selected value for minimum interference on LnP-coef (elasticity) fit_intercept=False, # no need, since we have OneHot encodings without drop max_iter=100_000, )), ], verbose=True) mdl_basic_ols = Pipeline([ ('feat_proc', feature_generator_basic), ('reg', linear_model.Ridge( alpha=1e-20, # l2 penalty strength, "very small" fit_intercept=False, max_iter=100_000, )), ], verbose=True) mdl_basic.fit( df_mdl[['LnP', 'StockCode', 'Date', 'Country']], df_mdl['Quantity'] # Poisson regression has log-link, so LnQ is implicit in loss function )

但是训练数据不跟之前v1一样，不需要分组，直接用原始数据：

柏松回归中LnP的回归系数为 -2.87559，

Ridge—OLS回归中LnP的回归系数为 -1.79945，

尝试下来各个方法得到的结果差异很大。

2.4 v3版求解价格弹性：DML2.4.1 DML数据准备 + 建模 + 求残差

因为不同产品的单价差异很大，所以对于同一维度的单价需要减去该维度的单价均值：

dLnP_{i,t}=log(p_{i,t})-log(\overline p_{i})

这里，消除数据差异的方法，在盒马论文里面是：

Y_i^\text{o} 是常规渠道产品i 近期的平均销量Y_i/Y_i^\text{nor} 代表了折扣价格使得销量增加的百分比，因为不同商品销量差异很大，所以比率会比绝对值更有用代码语言：javascript复制df_mdl['dLnP'] = np.log(df_mdl.UnitPrice) - np.log(df_mdl.groupby('StockCode').UnitPrice.transform('mean')) df_mdl['dLnQ'] = np.log(df_mdl.Quantity) - np.log(df_mdl.groupby('StockCode').Quantity.transform('mean'))

混杂因子也做了一些处理：

季节性变量：该价格处于第几月、处于月里第几天和周里第几天产品上线的时长：用当期时间减去该产品的最小时间sku 的价格水平：单个sku内的价格中位数代码语言：javascript复制df_mdl = df_mdl.assign( month = lambda d: d.Date.dt.month, DoM = lambda d: d.Date.dt.day, DoW = lambda d: d.Date.dt.weekday, stock_age_days = lambda d: (d.Date - d.groupby('StockCode').Date.transform('min')).dt.days, sku_avg_p = lambda d: d.groupby('StockCode').UnitPrice.transform('median') )

有了混杂因子，lnp,lnq，来看看DML过程：

构建lnp-W ，lnq-W 计算残差代码语言：javascript复制# 混杂因子针对Q \ P 分别建模 model_y = Pipeline([ ('feat_proc', feature_generator_full), ('model_y', RandomForestRegressor(n_estimators=50, min_samples_leaf=3, n_jobs=-1, verbose=0)) # n_samples_leaf/n_estimators is set to reduce model (file) size and runtime # larger models yield prettier plots. ]) model_t = Pipeline([ ('feat_proc', feature_generator_full), ('model_t', RandomForestRegressor(n_estimators=50, min_samples_leaf=3, n_jobs=-1, verbose=0)) ]) # 上述模型得到预估值 # Get first-step, predictions to residualize ("orthogonalize") with (in-sample for now) q_hat = model_y.predict(df_mdl) p_hat = model_t.predict(df_mdl) # 用观测值减去预测得到的值求解残差 df_mdl = df_mdl.assign( dLnP_res = df_mdl['dLnP'] - p_hat, dLnQ_res = df_mdl['dLnQ'] - q_hat, )2.4.2 三块模型对比

此时经过数据处理，数据集中就有三种数据类型：

对数对数+去均值化对数+去均值化+求残差

然后三组数据，按照v1版的处理方式，先分段，后利用OLS求价格弹性：

代码语言：javascript复制# 初始ols模型 old_fit = binned_ols( df_mdl, x='LnP', y='LnQ', n_bins=15, ) # 初始去均值化后的ols模型 old_fit = binned_ols( df_mdl, x='dLnP', y='dLnQ', n_bins=15, plot_ax=plt.gca(), ) # 残差拟合的ols模型 old_fit = binned_ols( df_mdl, x='dLnP_res', y='dLnQ_res', n_bins=15, plot_title='Causal regression naively, with item controls, and after DML.', plot_ax=plt.gca() )

此时经过数据处理，数据集中就有三种数据类型，三者的价格弹性对比：

对数：\theta=-1.7 对数+去均值化：\theta=-1.7 对数+去均值化+求残差：\theta=-1.819

当然OLS还有截距项，绘图可得：

这里原文也给出了，DML求解过程中，两个随机森林模型的特征重要性：

代码语言：javascript复制feat_imp = pd.DataFrame({ 'feat': get_feat_generator_names(model_y['feat_proc']), 'importance_q': model_y['model_y'].feature_importances_, 'importance_p': model_t['model_t'].feature_importances_, }).set_index('feat') feat_imp.sort_values(by='importance_p').iloc[-15:].plot.barh( figsize=(5, 8), title='feature importances for DML estimators of treatment(p) and outcome(q)' )

2.4.3 稳健性评估

这章主要学习到的：

一种数据筛选的原则，残差正交化后，dLnP_{res}

总是很小，因此为了减少噪音，我们将丢弃所有非常小的价格变化观察值，它们不包含太多信息

训练数据分成多k-fold来检验弹性系数的稳定性

那么在盒马那篇文章里面来看一下这个图，使用training data的比例往上几个模型的稳定性分布情况

模型的预测推断结果是\hat \theta=\frac{dLnQ_{res}}{dLnP_{res}} 但是残差正交化后，dLnP_{res} 总是很小，因此为了减少噪音，我们将丢弃所有非常小的价格变化观察值，它们不包含太多信息。

Chernozhukov 提出了一个改进的 DML，传统的标准 OLS 方法估计\hat \theta=(\widetilde P^T \widetilde P)^{-1}\widetilde P^T \widetilde Q 但改进的\hat \theta=(\widetilde P^T P)^{-1}\widetilde P^T \widetilde Q^T

即第二个 P 矩阵用未残差化的。最后采取 2-fold 得到平均值使得结果更稳健，最终弹性系数结果为 -1.89

代码语言：javascript复制old_fit = binned_ols( df_mdl, x='dLnP', y='dLnQ', n_bins=15, plot_ax=plt.gca(), ) plt.gca().set( xlabel='log(price)', ylabel='log(quantity)', ) plt.gca().axvline(0, color='k', linestyle=':') plt.gca().axhline(0, color='k', linestyle=':') elast_estimates = list() for idx_aux, idx_inf in KFold(n_splits=2, shuffle=True).split(df_mdl): df_aux = df_mdl.iloc[idx_aux] df_inf = df_mdl.iloc[idx_inf].copy() # step 1: aux models and residualize in inferential set print('fitting model_y') model_y.fit(df_aux, df_aux.dLnQ) print('fitting model_t') model_t.fit(df_aux, df_aux.dLnP) df_inf = df_inf.assign( dLnP_res = df_inf['dLnP'] - model_t.predict(df_inf), dLnQ_res = df_inf['dLnQ'] - model_y.predict(df_inf), ) binned_ols( df_inf, x='dLnP_res', y='dLnQ_res', n_bins=15, plot_ax=plt.gca(), label='fold' ) # ignore observations where we residualized away all variation in price mask = (~(df_inf.dLnP_res.abs() < 0.01)) df_inf_censored = df_inf[mask] # step 2.1: Chernozhukov DML inference elast = ( df_inf_censored['dLnP_res'].dot(df_inf_censored['dLnQ_res']) / df_inf_censored['dLnP_res'].dot(df_inf_censored['dLnP']) # the last part here deviates from standard OLS solution ) print('DML elast: ', elast) elast_estimates.append(elast) print('OLS elasticity for comparison:', df_inf_censored['dLnP_res'].dot(df_inf_censored['dLnQ_res']) / df_inf_censored['dLnP_res'].dot(df_inf_censored['dLnP_res']) ) print("DML efficient estimate of elasticity:", np.mean(elast_estimates))

这里远程再对比一下盒马的那篇，貌似媲美他们的半参数模型？

【本文地址】

公司简介

联系我们

今日新闻

点击排行

实验室常用的仪器、试剂和: 说到实验室常用到的东西，主要就分为仪器、试剂和耗

不用再找了，全球10大实验: 01、赛默飞世尔科技（热电）Thermo Fisher Scientif

三代水柜的量产巅峰T-72坦: 作者：寞寒最近，西边闹腾挺大，本来小寞以为忙完这

通风柜跟实验室通风系统有: 说到通风柜跟实验室通风，不少人都纠结二者到底是不

集消毒杀菌、烘干收纳为一: 厨房是家里细菌较多的地方，潮湿的环境、没有完全密

实验室设备之全钢实验台如: 全钢实验台是实验室家具中较为重要的家具之一，很多

图片新闻

实验室药品柜的特性有哪些: 实验室药品柜是实验室家具的重要组成部分之一，主要

小学科学实验中有哪些教学: 计算机计算器一般打孔器打气筒仪器车显微镜

实验室各种仪器原理动图讲: 1.紫外分光光谱UV分析原理：吸收紫外光能量，引起分

高中化学常见仪器及实验装: 1、可加热仪器：2、计量仪器：（1）仪器A的名称：量

微生物操作主要设备和器具: 今天盘点一下微生物操作主要设备和器具，别嫌我啰嗦

浅谈通风柜使用基本常识: 　众所周知，通风柜功能中最主要的就是排气功能。在

因果推断与反事实预测

因果推断与反事实预测

今日新闻

点击排行

推荐新闻

图片新闻

专题文章