【python因果推断库2】使用 PyMC 模型进行差分-in-差分（Difference in Differences, DID）分析
目录
使用 PyMC 模型进行差分-in-差分（Difference in Differences, DID）分析
导入数据
分析
使用 PyMC 模型建模银行业数据集
导入数据
分析 1 - 经典 2×2 差分-in-差分 (DiD)
分析 2 - 具有多个干预前后观测值的差分-in-差分 (DiD) 分析
使用 PyMC 模型进行差分-in-差分（Difference in Differences, DID）分析
import arviz as az
import causalpy as cp
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = 'retina'
seed = 42导入数据
df = cp.load_data("did")
df.head()
分析
`random_seed` 这个关键词参数对于 PyMC 采样器来说并不是必需的。我们在这里使用它是为了确保结果是可以重现的。
result = cp.pymc_experiments.DifferenceInDifferences(
df,
formula="y ~ 1 + group*post_treatment",
time_variable_name="t",
group_variable_name="group",
model=cp.pymc_models.LinearRegression(sample_kwargs={"random_seed": seed}),
)
fig, ax = result.plot()
result.summary()===========================Difference in Differences============================
Formula: y ~ 1 + group*post_treatment
Results:
Causal impact = 0.5, $CI_{94\%}$[0.4, 0.6]
Model coefficients:
Intercept 1.1, 94% HDI [1, 1.1]
post_treatment[T.True] 0.99, 94% HDI [0.92, 1.1]
group 0.16, 94% HDI [0.094, 0.23]
group:post_treatment[T.True] 0.5, 94% HDI [0.4, 0.6]
sigma 0.082, 94% HDI [0.066, 0.1]
ax = az.plot_posterior(result.causal_impact, ref_val=0)
ax.set(title="Posterior estimate of causal impact");
使用 PyMC 模型建模银行业数据集
本笔记本分析了来自 Richardson (2009) 的历史银行业关闭数据，并将其作为差分-in-差分分析的一个案例研究，该案例研究来源于优秀的书籍《Mastering Metrics》（Angrist 和 Pischke, 2014）。在这里，我们复制了这项分析，但是使用了贝叶斯推断的方法。
import arviz as az
import pandas as pd
import causalpy as cp
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = 'retina'
seed = 42导入数据
原始数据集包含一个日期列，这个列中的数字无法直接解读——对于我们分析而言只需要年份列。数据集中还有 `bib6`, `bio6`, `bib8`, `bio8` 这几列。我们知道数字 6 和 8 分别代表第 6 和第 8 联邦储备区。我假设 `bib` 表示“营业中的银行”，所以我们保留这些列。数据是以天为单位的，但我们将把它转换成年为单位。从 Angrist 和 Pischke (2014) 一书中的图 5.2 来看，他们似乎展示了每年营业银行数量的中位数。现在让我们加载数据并执行这些步骤。
df = (
cp.load_data("banks")
# just keep columns we want
.filter(items=["bib6", "bib8", "year"])
# rename to meaningful variables
.rename(columns={"bib6": "Sixth District", "bib8": "Eighth District"})
# reduce from daily resolution to examine median banks open by year
.groupby("year")
.median()
)
treatment_time = 1930.5
# set treatment time to zero
df.index = df.index - treatment_time
treatment_time = 0ax = df[["Sixth District", "Eighth District"]].plot(style="o-")
ax.set(ylabel="Number of banks in business")
ax.axvline(x=treatment_time, c="r", lw=3, label="intervention")
ax.legend();让我们可视化我们现在得到的数据。这与 Angrist 和 Pischke (2014) 一书中的图 5.2 完全匹配。
只需几个最后的数据处理步骤，就可以使数据适合于分析。我们将把数据从宽格式转换为长格式。然后我们将添加一个新的列 `treated`，用来指示已经实施处理的观测值。
df.reset_index(level=0, inplace=True)
df_long = pd.melt(
df,
id_vars=["year"],
value_vars=["Sixth District", "Eighth District"],
var_name="district",
value_name="bib",
).sort_values("year")
# We also need to create a column called `unit` which labels each distinct
# unit. In our case we just have two treatment units (each district). So
# we can build a `unit` column from `district`.
df_long["unit"] = df_long["district"]
# We also need to create a `post_treatment` column to define times after
# the intervention.
df_long["post_treatment"] = df_long.year >= treatment_time
df_long
# Dummy coding for district
df_long = df_long.replace({"district": {"Sixth District": 1, "Eighth District": 0}})
df_long
分析 1 - 经典 2×2 差分-in-差分 (DiD)
首先，我们只分析 1930 年和 1931 年的数据。这样我们就只有一个干预前的测量和一个干预后的测量。
我们将使用公式：`bib ~ 1 + district * post_treatment`，这等价于以下的期望值模型：
让我们逐一理解这些内容：
- 是第 i个观测值的结果（营业中的银行数量）的期望值。
- 是一个截距项，用来捕捉对照组在干预前营业中银行的基础数量。
- `district` 是一个虚拟变量，所以 将代表地区的主要效应，即相对于对照组，处理组的任何偏移量。
- `post_treatment` 也是一个虚拟变量，用来捕捉无论是否接受处理，在处理时间之后结果的任何变化。
- 两个虚拟变量的交互作用 `district:post_treatment` 只会在干预后处理组中取值为 1。因此， 将代表我们估计的因果效应。
result1 = cp.pymc_experiments.DifferenceInDifferences(
df_long[df_long.year.isin([-0.5, 0.5])],
formula="bib ~ 1 + district * post_treatment",
time_variable_name="post_treatment",
group_variable_name="district",
model=cp.pymc_models.LinearRegression(
sample_kwargs={"target_accept": 0.98, "random_seed": seed}
),
)这里我们遇到了一些发散的情况，这是不好的迹象。这很可能与我们只有4个数据点却有5个参数有关。这对于贝叶斯分析来说并不总是致命的（因为我们有先验分布），不过当我们遇到发散的样本时，这确实是一个警告信号。
使用下面的代码，我们可以看到我们遇到了经典的“漏斗问题”，因为当采样器探索测量误差（即 σ 参数）接近零的值时出现了发散。
az.plot_pair(result1.idata, var_names="~mu", divergences=True);为了进行“正规”的工作，我们需要解决这些问题以避免出现发散情况，例如，可以尝试探索不同的 σ 参数的先验分布。
fig, ax = result1.plot(round_to=3)
result1.summary()===========================Difference in Differences============================
Formula: bib ~ 1 + district * post_treatment
Results:
Causal impact = 19, $CI_{94\%}$[15, 22]
Model coefficients:
Intercept 165, 94% HDI [163, 167]
post_treatment[T.True] -33, 94% HDI [-36, -30]
district -30, 94% HDI [-32, -27]
district:post_treatment[T.True] 19, 94% HDI [15, 22]
sigma 0.84, 94% HDI [0.085, 2.2]
ax = az.plot_posterior(result1.causal_impact, ref_val=0, round_to=3)
ax.set(title="Posterior estimate of causal impact");
分析 2 - 具有多个干预前后观测值的差分-in-差分 (DiD) 分析
现在我们将对整个数据集进行差分-in-差分分析。这种方法与{术语}CITS（比较性中断时间序列）具有相似之处，其中涉及单个对照组随时间的变化。虽然这种区分稍微有些武断，但我们根据是否有足够的时间序列数据让CITS能够捕捉时间序列模式来区别这两种技术。
我们将使用公式：`bib ~ 1 + year + district*post_treatment`，这等价于以下的期望值模型：
与上面的经典 2×2 差分-in-差分模型相比，这里唯一的改变是增加了年份的主要效应。因为年份是数值编码的（而不是分类编码），它可以捕捉结果变量随时间发生的任何线性变化。
result2 = cp.pymc_experiments.DifferenceInDifferences(
df_long,
formula="bib ~ 1 + year + district*post_treatment",
time_variable_name="year",
group_variable_name="district",
model=cp.pymc_models.LinearRegression(
sample_kwargs={"target_accept": 0.95, "random_seed": seed}
),
)
fig, ax = result2.plot(round_to=3)
result2.summary()===========================Difference in Differences============================
Formula: bib ~ 1 + year + district*post_treatment
Results:
Causal impact = 20, $CI_{94\%}$[15, 26]
Model coefficients:
Intercept 160, 94% HDI [157, 164]
post_treatment[T.True] -28, 94% HDI [-33, -22]
year -7.1, 94% HDI [-8.5, -5.7]
district -29, 94% HDI [-34, -24]
district:post_treatment[T.True] 20, 94% HDI [15, 26]
sigma 2.4, 94% HDI [1.7, 3.2]
通过观察交互项，它可以捕捉干预措施的因果影响，我们可以看出干预似乎挽救了大约20家银行。尽管对此存在一定的不确定性，但我们可以在下方看到这一影响的完整后验估计。
ax = az.plot_posterior(result2.causal_impact, ref_val=0, round_to=3)
ax.set(title="Posterior estimate of causal impact");