【python因果推断库3】使用 CausalPy 进行贝叶斯geolift 分析

目录

导入数据

丹麦的销售额是否有地理提升(GeoLift)?

结果


本笔记本介绍如何使用 CausalPy 的贝叶斯{术语}合成控制功能来评估“地理提升”(GeoLift)。我们的假设情景如下:

你是一家在欧洲运营的公司的数据科学家。你得到了一份过去四年的历史销售量数据集,数据按千单位计算,按周频率收集,并按国家细分。数据涵盖了过去四年的时间。

从2022年初开始,市场部门决定翻新丹麦所有的门店。到了2022年底,公司希望你能评估这次翻新计划是否增加了销售额。如果你告诉他们门店翻新计划提高了销售额,那么他们会将这一计划推广到其他国家。虽然没有人明说,但在你的潜意识里你担心如果报告翻新能提高销售额但实际上未来并没有发生这种情况,那么公司的利润将会下降,你的股票价值也会减少,甚至你的工作也可能面临风险。

你的上司非常明白这些担忧。她也有同样的顾虑。她知道虽然很容易建立门店翻新与销售额变化之间的关联,但我们真正想知道的是门店翻新是否导致了销售额的增加。

我们知道提出因果主张的最佳方式是运行随机对照试验(有时称为A/B测试)。如果我们随机选择欧洲各地的门店(或选择一个国家),那么也许A/B测试就能完成任务。但我们并没有随机选择丹麦,所以我们担心可能存在混淆变量。

但我们听说过合成控制方法以及所谓的GeoLift。经过一番研究后,我们决定这就是我们要做的事情。但我们特别关心我们对检测到的任何提升程度的确信程度,所以我们希望使用贝叶斯方法并获得易于解释的贝叶斯可信区间。你找到了一个叫做 CausalPy 的库,正好支持这种应用场景,感到非常高兴。

import arviz as az
import matplotlib.dates as mdates
import pandas as pdimport causalpy as cp
%load_ext autoreload
%autoreload 2
pd.set_option("display.precision", 2)
seed = 42

导入数据

CausalPy 包含了一个适用于探索地理提升测试概念的示例(模拟)数据集。我们所需要做的就是加载这个数据集,适当地在 Pandas 数据框中设置观察日期,并定义处理时间。

df = (cp.load_data("geolift1").assign(time=lambda x: pd.to_datetime(x["time"])).set_index("time")
)treatment_time = pd.to_datetime("2022-01-01")
df.head()

在我们的数据集中,各列代表我们在欧洲运营的不同国家。我们还有一个索引,用于标记每一行的日期 —— 观测频率为每周一次。表格中的值代表销售量,并以千为单位。因此,一个2.4的值代表每周售出2,400单位。

那么,让我们来绘制一下这个图表。

untreated = list(set(df.columns).difference({"Denmark"}))
ax = df[untreated].plot(color=[0.7, 0.7, 0.7])
df["Denmark"].plot(color="r", ax=ax)
ax.axvline(treatment_time, color="k", linestyle="--")
ax.set(title="Observed data", ylabel="Sales volume (thousands)");

看起来挺不错的,但也有些令人失望。从视觉检查中并不明显看出,在门店翻新后(由虚线表示),丹麦的销售数据有明显的变化。这一点因为数据中存在的年度季节性(每个国家都不同)以及周销售数据本身的随机性而变得更糟。

为了更好地分析丹麦门店翻新后的销售变化,我们可以考虑对数据做一些处理,例如使用移动平均来平滑数据,或者建立一个时间序列模型来分离季节性和趋势成分。这样可以帮助我们更清晰地观察到翻新前后销售的具体变化。

丹麦的销售额是否有地理提升(GeoLift)?

为了计算门店翻新所带来的因果效应(如果有),我们需要比较干预后丹麦的实际销售额与如果没有进行干预时丹麦的反事实销售额。这就是为什么称之为反事实的原因——我们在丹麦确实进行了门店翻新,所以这是一个完全假设的情景,与事实相反。但如果我们可以测量(或者更现实地说是估算)这个反事实情形,那就会成为我们的对照组。

在这种情况下,我们将生成一个合成控制组,这是一种用来估算如果门店没有翻新的情况下丹麦反事实销售额的技术。你可以阅读合成控制维基百科页面了解更多关于合成控制方法的信息,但基本思想如下。对于熟悉传统(非贝叶斯)建模方法的人来说,基本的合成控制算法可以这样描述:

import my_custom_scikit_learn_model as weighted_combination# fit the model to pre-intervention (training) data
weighted_combination.fit(X_train, y_train)
# estimate the counterfactual given post-intervention (test) data
counterfactual = weighted_combination.predict(X_test)
# estimate the causal impact
causal_impact = y_test - counterfactual
# cumulative causal impact
np.cumsum(causal_impact)

这里并没有什么魔法——我们只是将未接受处理的单位作为加权和来估计一个合成的丹麦。我们基于干预前观察到的“训练”数据来做这件事。然后,我们使用这个加权和模型根据干预后未处理国家的“测试”数据来预测我们的合成丹麦。之后,我们可以简单地将这个反事实的估计与丹麦观察到的销售数据进行比较,从而得到因果影响的估计值。我们还可以(可选地)计算累积因果影响来回答这个问题:“门店翻新在2022年带来了多少额外的销售额?”

我们可以使用CausalPy的API来运行这个过程,但是采用贝叶斯推断方法如下:

对于PyMC采样器的random_seed关键字参数并不是必需的。我们在这里使用它是为了使结果可以重现。

formula = """Denmark ~ 0 + Austria + Belgium + Bulgaria + Croatia + Cyprus + Czech_Republic
"""result = cp.pymc_experiments.SyntheticControl(df,treatment_time,formula=formula,model=cp.pymc_models.WeightedSumFitter(sample_kwargs={"target_accept": 0.95, "random_seed": seed}),
)

结果

让我们使用ArviZ来检查每个国家的beta权重以及测量标准差sigma的后验参数估计值。

az.plot_forest(result.idata, var_names=["~mu"], figsize=(8, 3), combined=True);

 现在我们可以使用从 CausalPy 返回的结果对象的绘图方法。这将为我们提供一个非常详细的可视化输出。

fig, ax = result.plot(plot_predictors=False)# formatting
ax[2].tick_params(axis="x", labelrotation=-90)
ax[2].xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
ax[2].xaxis.set_major_locator(mdates.YearLocator())
for i in [0, 1, 2]:ax[i].set(ylabel="Sales (thousands)")

通过创建一个简单的模型公式并调用一次 CausalPy,我们就能够评估出处理单元产生的提升效果。

在这个例子中,测量噪声相当多,但由于我们在这里使用的是贝叶斯推断方法,因此我们对自己的不确定性有了精确且合理的原则性量化。

我们可以看到,对于预处理数据的贝叶斯值大约是 0.5。这个值虽然不是特别好,但对于真实世界的数据来说已经很不错了。它表明线性加权组合模型(合成控制的核心)在构建虚拟(即合成)丹麦直到处理期之前的工作表现还算合理。

这个合成控制的丹麦就是我们估计的反事实情况——如果商店翻新项目没有实施的话,销售额会是多少。通过对比我们可以估算因果影响,也可以称之为“地理位置提升”。

在实施后的一年内,我们可以看到丹麦的销售额累计因果影响接近 10,000 单位。让我们更详细地看看这一点。下面我们将查看时间序列末尾时累计因果影响的后验分布,在该方案实施一年之后。

# get index of the final time point
index = result.post_impact_cumulative.obs_ind.max()
# grab the posterior distribution of the cumulative impact at this final time point
last_cumulative_estimate = result.post_impact_cumulative.sel({"obs_ind": index})
# get summary stats
ax = az.plot_posterior(last_cumulative_estimate, figsize=(8, 4))
ax.set(title="Estimated cumulative causal impact (at end of 2022)",xlabel="Sales (thousands)",
);

如果我们愿意,我们也可以以数值形式提取这些统计信息:

如此,在我们因果建模工作的最后,我们可以向我们的上司汇报如下:“我们认为商店翻新计划对推动总共约 9,150 个额外销售起到了因果作用。但我们对确切的额外销售数量有所不确定——我们的 94% 可信区间范围是从 7,420 到 10,790。”

当然,还有一些需要注意的限制条件。我们所进行的分析假设唯一可能选择性地影响丹麦销售额的重大事件就是商店翻新项目。如果这是一个合理的假设,那么我们在提出因果主张时可能会相对稳固。但如果还有其他事件选择性地影响了一些单位(国家)而不是其他的单位,那么我们在提出主张时就需要更加谨慎,并可能需要采取更深入的建模方法。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/415113.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

图文解析保姆级教程: IDEA里面创建SpringBoot工程、SpringBoot项目的运行和测试、实现浏览器返回字符串

文章目录 一、创建SpringBoot工程(需要联网)二、 定义请求处理类三、运行测试 此教程摘选自我的笔记:黑马JavaWeb开发笔记13——Springboot入门(创建、运行&测试项目)、Http协议(请求&响应协议&…

Json数据解析报错 -TAB

表现: n8n 解析服务器的数据 报错 json 解析错误 原理: tab键 在代码中为 string tab \t解决办法:tab键替换4个空格 string tab "\t" tab.replaceAll("\t", " ")问题: tab 键 和 空格 在普…

卷积公式的几何学理解

1、Required Knowledge 1.1、概率密度函数 用于描述连续型随机变量在不同取值上的概率密度,记作 f ( x ) f(x) f(x)。 如随机变量 X X X的分布为正态分布,则其概率密度函数为: f ( x ) 1 σ 2 π e − ( x − μ ) 2 2 σ 2 f(x)\frac{1}…

记忆化搜索【上】

509. 斐波那契数 题目链接:斐波那契数 递归(暴搜) 斐波那契数列,最传统的解法,采用递归: class Solution { public:int fib(int n){return dfs(n);}int dfs(int n){if(n 0 || n 1)return n;return d…

大数据-114 Flink DataStreamAPI 程序输入源 自定义输入源 Rich并行源 RichParallelSourceFunction

点一下关注吧!!!非常感谢!!持续更新!!! 目前已经更新到了: Hadoop(已更完)HDFS(已更完)MapReduce(已更完&am…

linux 高级IO

IO等(要进行io是要有条件的,要有数据或者有空间)拷贝。高效体现在等待的时间所占比重越低越高效。 阻塞IO:数据没有就绪,read不返回。在内核将数据准备好之前, 系统调用会一直等待。所有的套接字, 默认都是阻塞方式。…

nginx容器映射配置文件后,启动一直报错提示:failed (13: Permission denied)的排查

问题现象: 使用harbor 的install.sh 创建docker-compose之后,出现nginx容器一直重启。 查看日志发现是:配置文件无权限。报错信息如下: Sep 2 16:43:13 172.28.0.1 nginx[1344]: 2024/09/02 08:43:13 [emerg] 1#0: open() “/e…

百度地图绘制电子围栏(包括移动端绘制操作)以及检测坐标是否在电子围栏内

由于本人在PC端仅使用了多边形绘制,但矩形跟多边形用法基本一样,圆形并未使用,如不符合读者需求也可以参考一下。 绘制后得到的数据可能不同,但绘制方法仅仅是传递的参数不同。 关于给坐标数组在地图绘制图形的效果在移动端部分包…

深度学习系列74:语音中的mel谱

1 mel谱介绍 一个人说一句话,其 waveform 可以很不一样,但是 spectrogram 基本上会相似,甚至有人可以通过 spectrogram 来判断说话的内容。语谱图的横坐标是时间,纵坐标是频率,坐标点值为语音数据能量。由于是采用二维…

# 利刃出鞘_Tomcat 核心原理解析(十一)-- WebSocket -- 1

利刃出鞘_Tomcat 核心原理解析(十一)-- Tomcat 附加功能 WebSocket – 1 一、Tomcat专题 - WebSocket - 介绍 1、Tomcat 附加功能:websocket 介绍 1)websocket :是 HTML5 新增的协议,它的目的是在浏览器…

动态规划法-资源分配问题

动态规划法 - 资源分配问题 问题描述 把4个份额的资源分配给3个工程,给定利润表如下表所示,写出资源的最优分配方案的求解过程。 4份资源分配给3个工程的利润表 步骤一:求各个阶段不同分配份额时的最大利润及分配份额 目标 我们的目标是…

53 mysql pid 文件的创建

前言 接上一篇文章 mysql 启动过程中常见的相关报错信息 在 mysql 中文我们在 “service mysql start”, “service mysql stop” 经常会碰到 mysql.pid 相关的错误信息 比如 “The server quit without updating PID file” 我们这里来看一下 mysql 中 mysql.pid 文件的…

微积分复习笔记 Calculus Volume 1 - 1.3Trigonometric Functions

1.3 Trigonometric Functions - Calculus Volume 1 | OpenStax

自己开发完整项目一、登录功能-05(动态权限控制)

一、上节回顾 在上一节中,我们介绍了如何通过数据库查询用户的权限,并对方法级别的接口使用注解的方式进行权限控制,之后通过用户携带的tocken进行解析权限,判断是否可以访问。 具体步骤: 1.在查询用户信息的时候将用户…

map和set的封装

目录 一、红黑树的改造 1.1节点的定义 二、红黑树的迭代器 2.1用节点封装迭代器 2.2迭代器的实现 2.3map和set的迭代器 三、insert的返回值 3.1insert返回值的用处 3.2operator[ ]的实现 四、key不能修改的问题 封装map和set一般分为六步: 封装map和set …

MFC生成dll的区别

主要分三种: A. 动态链接库(dll) B.具有导出项的(dll)动态链接库 C.MFC动态链接库 对比项目:可以根据需要选择哪种dll方便 添加自定义导出功能Demo 1. 添加导出实现接口: A. 导出需要具有:__declspec(dllexport) B. 按照C语言…

Javascript LeeCode选题(汉诺塔求解)

LeeCode选题 汉诺塔递归求解move移动函数hanoi函数main方法测试代码:代码实现 汉诺塔递归求解 汉诺塔介绍: 汉诺塔的的图形(从上到下1,2,3个)实现: 这里我们可以看到因为必须要将第n个移动到…

Spring中基于redis stream 的消息队列实现方法

本文主要介绍了消息队列的概念性质和应用场景,介绍了kafka、rabbitMq常用消息队列中间件的应用模型及消息队列的实现方式,并实战了在Spring中基于redis stream 的消息队列实现方法。 一、消息队列 消息队列是一种进程间通信或者同一个进程中不同线程间的…

uni-app 获取当前位置的经纬度以及地址信息

文章目录 uni.getLocation(objc)获取经纬度和地址调试结果问题 uni-app 获取当前位置的经纬度以及地址信息 uni.getLocation(objc) uni-app官方文档定位API: uni.getLocation(OBJECT) uni.getLocation({type: wgs84,success: function (res) {console.log(当前位置的经度&…

【系统架构设计】嵌入式系统设计(1)

【系统架构设计】嵌入式系统设计(1) 嵌入式系统概论嵌入式系统的组成硬件嵌入式处理器总线存储器I/O 设备与接口 软件 嵌入式开发平台与调试环境交叉平台开发环境交叉编译环境调试 嵌入式系统概论 嵌入性、专用性、计算机系统是嵌入式系统的三个基本的核…