【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 pd

import 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。”

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