【python因果推断库14】饮酒年龄 - 贝叶斯分析

目录

饮酒年龄 - 贝叶斯分析

主效应模型

交互模型

将连续变量以治疗阈值为中心

饮酒年龄 - 贝叶斯分析

这个例子使用了回归断点设计来探讨最低合法饮酒年龄（在美国为21岁）对全因死亡率的因果效应。数据集来自carpenter2009effect 的一项研究。

import arviz as az
import matplotlib.pyplot as plt

import causalpy as cp
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = 'retina'
seed = 42
Load and process data.

df = (
    cp.load_data("drinking")
    .rename(columns={"agecell": "age"})
    .assign(treated=lambda df_: df_.age > 21)
)

主效应模型

首先，我们将考察一个简单的“主效应”模型。在这里，对于给定年龄 a（以月为单位），预期死亡率（以每千人每年的死亡数为单位）通过截距项、处理效应  和年龄效应  来建模。

其中 t(a) 描述了是否已经施加了“处理”。在这个例子中，它仅仅描述了给定年龄 a 是否超过了最低合法饮酒年龄 21 岁：

为了明确， 描述了预期死亡率随年龄变化的线性趋势。系数 是这一线性趋势在 (a = 0) 时的截距。这剩下 ，我们可以理解为在年龄阈值周围的线性趋势的不连续性。

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

result = cp.pymc_experiments.RegressionDiscontinuity(
    df,
    formula="all ~ 1 + age + treated",
    running_variable_name="age",
    model=cp.pymc_models.LinearRegression(
        sample_kwargs={"target_accept": 0.95, "random_seed": seed}
    ),
    treatment_threshold=21,
)

fig, ax = result.plot()

result.summary()

============================Regression Discontinuity============================
Formula: all ~ 1 + age + treated
Running variable: age
Threshold on running variable: 21

Results:
Discontinuity at threshold = 7
Model coefficients:
  Intercept      	106, 94% HDI [83, 128]
  treated[T.True]	7, 94% HDI [4.4, 9.5]
  age            	-0.66, 94% HDI [-1.8, 0.51]
  sigma          	2.4, 94% HDI [2, 2.9]

在这个主效应模型中，因果效应的大小等于我们对 的后验估计。让我们绘制参数估计（左侧）并放大处理主效应的后验分布（右侧）。

fig, ax = plt.subplots(1, 2, figsize=(10, 3))

az.plot_forest(result.idata.posterior, var_names="beta", ax=ax[0])
az.plot_posterior(
    result.idata.posterior.beta.sel(coeffs="treated[T.True]"),
    round_to=3,
    ax=ax[1],
)

ax[1].set(title="treated[T.True]");

我们可以看到这几乎完全匹配了第一个结果图中的“阈值处的不连续性”图。

它并不是完全相同，因为我们实际上是手动计算阈值处的不连续性。原因在于，在最简单的模型之外，阈值处的不连续性并不简单地等于这个参数。为了理解这一点，让我们来看一个稍微复杂一点的模型，其中我们加入了一个交互项。

交互模型

我们通过改变公式为 `all ~ 1 + age + treated + age:treated` 来添加一个交互项。这样就改变了统计模型为：

这个模型现在更加复杂，因为它可以根据数据来允许在21岁之前与之后的死亡率趋势发生变化。如果还不清楚的话，从下一个结果图中我们会看到阈值处的不连续性不再等于 参数。让我们运行这个模型来看看结果。

result2 = cp.pymc_experiments.RegressionDiscontinuity(
    df,
    formula="all ~ 1 + age + treated + age:treated",
    running_variable_name="age",
    model=cp.pymc_models.LinearRegression(sample_kwargs={"random_seed": seed}),
    treatment_threshold=21,
)

fig, ax = result2.plot()

我们现在可以看到，通过添加交互项，参数估计发生了很大的变化，阈值处的估计不连续性不再简单地由 参数给出。为了确认这一点，我们可以检查  的估计值（这对应于 treated[T.True] 系数）。

az.plot_forest(result2.idata.posterior, var_names="beta", figsize=(10, 3));

我们可以看到这个估计值现在与所需的“阈值处的不连续性”值相差很大。正因为如此，CausalPy 通过计算治疗阈值稍上方与稍下方的预测结果值之差，手动计算“阈值处的不连续性”。

将连续变量以治疗阈值为中心

另一种处理方法是将连续变量以阈值为中心，使得阈值（最低合法饮酒年龄）现在位于零点。这也使得参数更加易于解释。

df["age"] = df["age"] - 21
result3 = cp.pymc_experiments.RegressionDiscontinuity(
    df,
    formula="all ~ 1 + age + treated + age:treated",
    running_variable_name="age",
    model=cp.pymc_models.LinearRegression(sample_kwargs={"random_seed": seed}),
    treatment_threshold=0,
)

fig, ax = result3.plot()

fig, ax = plt.subplots(1, 2, figsize=(10, 3))

az.plot_forest(result3.idata.posterior, var_names="beta", ax=ax[0])
az.plot_posterior(
    result3.idata.posterior.beta.sel(coeffs="treated[T.True]"),
    round_to=3,
    ax=ax[1],
)

ax[1].set(title="treated[T.True]");