中断时间序列分析#

本笔记本重点介绍如何进行简单的贝叶斯中断时间序列分析。这在准实验设置中非常有用，其中干预措施应用于所有治疗单位。

例如，如果对网站进行了更改，并且您想知道该网站更改的因果影响，那么如果该更改被有选择性地且随机地应用于网站用户的一个测试组，那么您可能能够使用A/B测试方法来做出因果声明。

然而，如果网站的更改已经对所有用户实施，那么您就没有对照组。在这种情况下，您无法直接测量反事实，即如果没有进行网站更改，会发生什么。在这种情况下，如果您有足够多的时间点的数据，那么您可能能够利用中断时间序列方法。

感兴趣的读者可以参考优秀的教材 The Effect [亨廷顿-克莱因，2021]。第17章涵盖了作者更倾向于使用的“事件研究”，而不是中断时间序列术语。

因果DAG#

下面给出了一个简单的因果DAG用于中断时间序列，但请参阅[亨廷顿-克莱因，2021]以获取更一般的DAG。简而言之，它表示：

  • 结果受到时间（例如，随时间变化的其他因素）和治疗的影响。

  • 治疗受时间因果影响。

直观地，我们可以将该方法的逻辑描述为：

  • 我们知道结果随时间变化。

  • 如果我们建立一个模型来描述结果在治疗之前随时间的变化，那么我们就可以预测如果治疗没有发生，我们预期会发生什么。

  • 我们可以将这个反事实与干预后时期的观察结果进行比较。如果存在显著的差异，那么我们可以将其归因于干预的因果影响。

如果我们已经排除了在干预的同时（或之后）发生的其他合理原因，这是合理的。随着干预后时间的推移，这种解释变得越来越难以证明，因为可能发生了其他相关事件，这些事件可能提供了替代的因果解释。

如果这不能立即理解，我建议查看下面的示例数据图，然后重新阅读本节。

import arviz as az
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import xarray as xr

from scipy.stats import norm

%config InlineBackend.figure_format = 'retina'
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")

现在让我们定义一些辅助函数

Hide code cell content 
def format_x_axis(ax, minor=False):
    # major ticks
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y %b"))
    ax.xaxis.set_major_locator(mdates.YearLocator())
    ax.grid(which="major", linestyle="-", axis="x")
    # minor ticks
    if minor:
        ax.xaxis.set_minor_formatter(mdates.DateFormatter("%Y %b"))
        ax.xaxis.set_minor_locator(mdates.MonthLocator())
        ax.grid(which="minor", linestyle=":", axis="x")
    # rotate labels
    for label in ax.get_xticklabels(which="both"):
        label.set(rotation=70, horizontalalignment="right")


def plot_xY(x, Y, ax):
    quantiles = Y.quantile((0.025, 0.25, 0.5, 0.75, 0.975), dim=("chain", "draw")).transpose()

    az.plot_hdi(
        x,
        hdi_data=quantiles.sel(quantile=[0.025, 0.975]),
        fill_kwargs={"alpha": 0.25},
        smooth=False,
        ax=ax,
    )
    az.plot_hdi(
        x,
        hdi_data=quantiles.sel(quantile=[0.25, 0.75]),
        fill_kwargs={"alpha": 0.5},
        smooth=False,
        ax=ax,
    )
    ax.plot(x, quantiles.sel(quantile=0.5), color="C1", lw=3)


# default figure sizes
figsize = (10, 5)

生成数据#

本示例的重点是使用中断时间序列方法进行因果推断。因此，我们将使用一些相对简单的合成数据，这些数据只需要一个非常简单的模型。

treatment_time = "2017-01-01"
β0 = 0
β1 = 0.1
dates = pd.date_range(
    start=pd.to_datetime("2010-01-01"), end=pd.to_datetime("2020-01-01"), freq="M"
)
N = len(dates)


def causal_effect(df):
    return (df.index > treatment_time) * 2


df = (
    pd.DataFrame()
    .assign(time=np.arange(N), date=dates)
    .set_index("date", drop=True)
    .assign(y=lambda x: β0 + β1 * x.time + causal_effect(x) + norm(0, 0.5).rvs(N))
)
df
time y
date
2010-01-31 0 0.302698
2010-02-28 1 0.396772
2010-03-31 2 -0.433908
2010-04-30 3 0.276239
2010-05-31 4 0.943868
... ... ...
2019-08-31 115 12.403634
2019-09-30 116 13.185151
2019-10-31 117 14.001674
2019-11-30 118 13.950388
2019-12-31 119 14.383951

120 行 × 2 列

# Split into pre and post intervention dataframes
pre = df[df.index < treatment_time]
post = df[df.index >= treatment_time]

fig, ax = plt.subplots()
ax = pre["y"].plot(label="pre")
post["y"].plot(ax=ax, label="post")
ax.axvline(treatment_time, c="k", ls=":")
plt.legend();
../_images/49d3e0b1e33ef95f1905ed0e648bc790fadb21183cf29156f75e742ff162e1b2.png

在这个简单的数据集中，我们有一个向上的噪声线性趋势，并且由于这些数据是合成的，我们知道在干预时间点上结果有一个阶跃增加，并且这种效应随时间持续存在。

建模#

这里我们构建一个简单的线性模型。记住，我们正在构建一个干预前数据的模型，目标是它能够合理地预测如果未进行干预会发生什么。换句话说，我们不对干预后的观察结果进行任何方面的建模，例如截距的变化、斜率的变化或效果是暂时的还是永久的。

with pm.Model() as model:
    # observed predictors and outcome
    time = pm.MutableData("time", pre["time"].to_numpy(), dims="obs_id")
    # priors
    beta0 = pm.Normal("beta0", 0, 1)
    beta1 = pm.Normal("beta1", 0, 0.2)
    # the actual linear model
    mu = pm.Deterministic("mu", beta0 + (beta1 * time), dims="obs_id")
    sigma = pm.HalfNormal("sigma", 2)
    # likelihood
    pm.Normal("obs", mu=mu, sigma=sigma, observed=pre["y"].to_numpy(), dims="obs_id")

pm.model_to_graphviz(model)
../_images/07d35f8ce735c0bef9e204ad7e609371f84bfdcf6b938a252845dc6e65ff3004.svg

先验预测检查#

作为贝叶斯工作流程的一部分，我们将绘制先验预测图，以查看模型在观察到任何数据之前得出的结果。

with model:
    idata = pm.sample_prior_predictive(random_seed=RANDOM_SEED)

fig, ax = plt.subplots(figsize=figsize)

plot_xY(pre.index, idata.prior_predictive["obs"], ax)
format_x_axis(ax)
ax.plot(pre.index, pre["y"], label="observed")
ax.set(title="Prior predictive distribution in the pre intervention era")
plt.legend();

Sampling: [beta0, beta1, obs, sigma]
../_images/e484bdf704fd0f950e188ce71bbb35a8fa14ea627568f8d17aa82e12013dd09e.png

这在某种程度上是合理的，因为截距和斜率的先验足够广泛，以至于预测的观测值很容易包含实际数据。这意味着所选择的特定先验不会不当地限制后验参数估计。

推理#

绘制后验分布的样本，并记住我们仅为此干预前数据执行此操作。

with model:
    idata.extend(pm.sample(random_seed=RANDOM_SEED))

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [beta0, beta1, sigma]
100.00% [8000/8000 00:01<00:00 Sampling 4 chains, 0 divergences] 
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 1 seconds.

az.plot_trace(idata, var_names=["~mu"]);
../_images/75103b917521426142b814b8dfa564db306f82ce5bd9f1ad279ccbb33eeccdc6.png

后验预测检查#

贝叶斯工作流程的另一个重要方面是绘制模型的后验预测图，使我们能够看到模型对已观察数据的回溯预测效果如何。正是在这一点上，我们可以决定模型是否过于简单（那么我们会在模型中增加更多复杂性）或者是否已经足够好。

with model:
    idata.extend(pm.sample_posterior_predictive(idata, random_seed=RANDOM_SEED))

fig, ax = plt.subplots(figsize=figsize)

az.plot_hdi(pre.index, idata.posterior_predictive["obs"], hdi_prob=0.5, smooth=False)
az.plot_hdi(pre.index, idata.posterior_predictive["obs"], hdi_prob=0.95, smooth=False)
ax.plot(pre.index, pre["y"], label="observed")
format_x_axis(ax)
ax.set(title="Posterior predictive distribution in the pre intervention era")
plt.legend();

Sampling: [obs]
100.00% [4000/4000 00:00<00:00]
../_images/24082c83f5d3ecd6c1c0484718b05667f5ea47c1594abc40d2a6e3a60b223b81.png

下一步并不是严格必要的，但我们可以计算模型反向预测与数据之间的差异，以查看误差。这有助于识别任何意外的反向预测干预前数据的能力不足。

Hide code cell source 
# convert outcome into an XArray object with a labelled dimension to help in the next step
y = xr.DataArray(pre["y"].to_numpy(), dims=["obs_id"])

# do the calculation by taking the difference
excess = y - idata.posterior_predictive["obs"]

fig, ax = plt.subplots(figsize=figsize)
# the transpose is to keep arviz happy, ordering the dimensions as (chain, draw, time)
az.plot_hdi(pre.index, excess.transpose(..., "obs_id"), hdi_prob=0.5, smooth=False)
az.plot_hdi(pre.index, excess.transpose(..., "obs_id"), hdi_prob=0.95, smooth=False)
format_x_axis(ax)
ax.axhline(y=0, color="k")
ax.set(title="Residuals, pre intervention");
../_images/d18924b015c58c354c1ea453f673ec1f45db248a3eca1d3a2969463ecbed546a.png

反事实推理#

现在我们将使用我们的模型来预测在没有干预的“假设”情景中的观察结果。

因此，我们使用干预后数据框中的时间数据更新模型，并运行后验预测采样来预测在这种反事实情景下我们将会观察到的观测值。我们也可以称之为“预测”。

with model:
    pm.set_data(
        {
            "time": post["time"].to_numpy(),
        }
    )
    counterfactual = pm.sample_posterior_predictive(
        idata, var_names=["obs"], random_seed=RANDOM_SEED
    )

Sampling: [obs]
100.00% [4000/4000 00:00<00:00]
Hide code cell source 
fig, ax = plt.subplots(figsize=figsize)

plot_xY(post.index, counterfactual.posterior_predictive["obs"], ax)
format_x_axis(ax, minor=False)
ax.plot(post.index, post["y"], label="observed")
ax.set(
    title="Counterfactual: Posterior predictive forecast of outcome if intervention not taken place"
)
plt.legend();
../_images/6e686adda5100776fc9b67e22b3746c48e0cb515c985b9cb816321f08a84d845.png

我们现在有了计算因果影响所需的成分。这仅仅是贝叶斯反事实预测与观测值之间的差异。

因果影响：自干预以来#

现在我们将使用反事实情景下的预测结果，并将其与观察到的结果进行比较，以得出我们的反事实估计。

# convert outcome into an XArray object with a labelled dimension to help in the next step
outcome = xr.DataArray(post["y"].to_numpy(), dims=["obs_id"])

# do the calculation by taking the difference
excess = outcome - counterfactual.posterior_predictive["obs"]

我们可以轻松计算累积的因果影响

# calculate the cumulative causal impact
cumsum = excess.cumsum(dim="obs_id")
Hide code cell source 
fig, ax = plt.subplots(2, 1, figsize=(figsize[0], 9), sharex=True)

# Plot the excess
# The transpose is to keep arviz happy, ordering the dimensions as (chain, draw, t)
plot_xY(post.index, excess.transpose(..., "obs_id"), ax[0])
format_x_axis(ax[0], minor=True)
ax[0].axhline(y=0, color="k")
ax[0].set(title="Causal impact, since intervention")

# Plot the cumulative excess
plot_xY(post.index, cumsum.transpose(..., "obs_id"), ax[1])
format_x_axis(ax[1], minor=False)
ax[1].axhline(y=0, color="k")
ax[1].set(title="Cumulative causal impact, since intervention");
../_images/66a6c14b531ea70dc05928d110200706908ae8e29bf312a96cdeeca168db4e84.png

就这样，我们使用 PyMC 通过中断时间序列方法完成了贝叶斯反事实推理！在短短几个步骤中，我们：

  • 构建了一个简单的模型来预测时间序列。

  • 基于干预前数据推断模型参数，运行先验和后验预测检查。我们注意到模型相当好。

  • 使用模型创建了反事实预测，即如果在干预时间后干预未发生，会发生什么情况。

  • 通过将观察到的结果与在没有干预的情况下反事实预期结果进行比较，计算了因果影响（以及累积因果影响）。

当然，在现实世界中，中断时间序列方法可以有更多的应用方式。例如，可能存在更多的时序结构，如季节性。如果是这样，我们可能希望使用特定的时间序列模型，而不仅仅是线性回归模型。还可能有其他信息丰富的预测变量需要纳入模型中。此外，一些设计不仅仅包括干预前和干预后时期（也称为A/B设计），还可能涉及干预不活跃、活跃，然后再次不活跃的时期（也称为ABA设计）。

参考资料#

[1] (1,2)

尼克·亨廷顿-克莱因。效果：研究设计和因果关系导论。查普曼和霍尔/CRC，2021年。

作者#

  • 由Benjamin T. Vincent于2022年10月撰写。

  • 由Benjamin T. Vincent于2023年2月更新，以在PyMC v5上运行

水印#

%load_ext watermark
%watermark -n -u -v -iv -w -p pytensor,aeppl,xarray

Last updated: Wed Feb 01 2023

Python implementation: CPython
Python version       : 3.11.0
IPython version      : 8.9.0

pytensor: 2.8.11
aeppl   : not installed
xarray  : 2023.1.0

pymc      : 5.0.1
xarray    : 2023.1.0
numpy     : 1.24.1
matplotlib: 3.6.3
arviz     : 0.14.0
pandas    : 1.5.3

Watermark: 2.3.1

许可证声明#

本示例库中的所有笔记本均在MIT许可证下提供，该许可证允许修改和重新分发，前提是保留版权和许可证声明。

引用 PyMC 示例#

要引用此笔记本，请使用Zenodo为pymc-examples仓库提供的DOI。

重要

许多笔记本是从其他来源改编的：博客、书籍……在这种情况下，您应该引用原始来源。

同时记得引用代码中使用的相关库。

这是一个BibTeX的引用模板：

@incollection{citekey,
  author    = "<notebook authors, see above>",
  title     = "<notebook title>",
  editor    = "PyMC Team",
  booktitle = "PyMC examples",
  doi       = "10.5281/zenodo.5654871"
}

渲染后可能看起来像：