在双重差分法（DID）实证分析中，平行趋势检验是判断政策效应是否可信的重要环节。
当事件研究图显示政策实施前的系数显著时，很多研究会尝试增加控制变量、更换变量定义、缩短事件窗口，甚至删除部分年份，希望让平行趋势检验“通过”。
但这些操作未必能够真正解决问题。
平行趋势不成立，往往意味着处理组与对照组在政策实施前就处于不同的发展轨迹。此时，更规范的处理方式应当是重新检查研究设计，尤其是：
当前选择的对照组，能否合理代表处理组未接受政策时的反事实变化？
本文使用 Python 演示如何重新定义 DID 对照组，将“从未接受处理的单位”替换为“未来会接受处理、但当前尚未接受处理的单位”，并比较两种研究设计下的平行趋势结果。
案例背景
DID的基本逻辑，是比较处理组和对照组在政策前后的变化差异。
假设结果变量为 (Y_{it})，基础模型可以写为：
[
Y_{it}
\alpha_i
+
\lambda_t
+
\beta DID_{it}
+
\gamma X_{it}
+
\varepsilon_{it}
]
其中：
DID真正需要回答的问题是：
如果处理组没有接受政策，它的结果变量原本会怎样变化？
由于处理组未接受政策时的结果无法直接观察，因此需要通过对照组构造反事实。
这意味着，对照组不能只是“没有接受政策的单位”，还需要与处理组具有相近的潜在发展趋势。
以贫困县扶持政策为例。
部分县之所以较早被纳入政策，往往是因为这些地区具有以下特征：
这些特征不仅影响一个县是否会进入政策处理组，也可能直接影响当地收入、就业和产业结构的变化趋势。
如果将较早进入政策的贫困县作为处理组，将所有普通非贫困县作为对照组，那么两组在政策实施前就可能处于不同的发展轨迹。
此时，事件研究模型中的事前系数显著，未必是因为模型中少放了某个控制变量，更可能是因为：
原始对照组无法代表处理组的反事实变化。
分析思路
本文比较两种对照组定义。
第一种定义
使用从未接受政策的地区作为对照组：
这种设定操作简单，但从未处理县可能与政策县存在明显的经济基础和发展趋势差异。
第二种定义
使用晚处理地区作为对照组：
在2016—2019年期间，晚处理县尚未接受政策，因此可以暂时作为早处理县的对照组。
与从未处理县相比，晚处理县和早处理县通常更可能具有相近的：
因此，晚处理组更有可能代表早处理组未接受政策时的反事实结果。
这里需要注意：
晚处理组只能在正式接受政策之前作为对照组。
一旦晚处理县在2020年进入政策，就不能继续把它当作未处理单位。
事件研究模型
为检验平行趋势，可以估计以下动态DID模型：
[
Y_{it}
\alpha_i
+
\lambda_t
+
\sum_{k\neq-1}
\beta_k
\left(
Treat_i
\times
1{t-T_i=k}
\right)
+
\gamma X_{it}
+
\varepsilon_{it}
]
其中，(k=-1)通常作为基准期。
政策实施前的系数主要用于判断处理组与对照组是否存在提前趋势：
[
\beta_{-5},\beta_{-4},\beta_{-3},\beta_{-2}
]
如果这些系数整体接近零且不显著，说明平行趋势假设相对合理。
政策实施后的系数：
[
\beta_0,\beta_1,\beta_2,\beta_3
]
则反映政策效应的动态变化。
Python中可以使用PyFixest估计包含个体固定效应、时间固定效应和聚类标准误的模型。该工具支持类似fixest的公式写法，也支持通过i()构造事件时间与处理组的交互项。
核心代码
1. 安装环境
pip install pandas numpy pyfixest matplotlib openpyxl
导入需要使用的库：
import re
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyfixest as pf
数据至少需要包含以下字段：
| |
|---|
county_id | |
year | |
outcome | |
policy | |
control1 | |
control2 | |
control3 | |
读取数据：
df = pd.read_excel("did_data.xlsx")
required_columns = [
"county_id",
"year",
"outcome",
"policy",
"control1",
"control2",
"control3",
]
missing_columns = [
column for column in required_columns
if column not in df.columns
]
if missing_columns:
raise ValueError(
f"数据缺少以下字段：{missing_columns}"
)
df = df.sort_values(
["county_id", "year"]
).reset_index(drop=True)
2. 确定首次处理年份
首先计算每个县第一次接受政策的年份。
first_treat = (
df.loc[df["policy"].eq(1)]
.groupby("county_id")["year"]
.min()
.rename("first_treat")
)
df = df.merge(
first_treat,
on="county_id",
how="left",
)
检查首次处理年份：
treatment_timing = (
df[["county_id", "first_treat"]]
.drop_duplicates()
.sort_values("first_treat")
)
print(treatment_timing.head(20))
对于样本期内从未接受政策的地区，first_treat为空值。
构造是否曾经接受政策的变量：
df["ever_treated"] = (
df["first_treat"].notna().astype(int)
)
df["never_treated"] = (
df["first_treat"].isna().astype(int)
)
3. 设置处理批次
假设：
EARLY_YEAR = 2016
LATE_YEAR = 2020
PRE_WINDOW = 5
POST_WINDOW = 3
构造早处理组变量：
df["early_group"] = (
df["first_treat"].eq(EARLY_YEAR)
).astype(int)
事件时间以早处理组的政策年份为基准：
df["event_time"] = df["year"] - EARLY_YEAR
保留政策实施前5年至实施后3年：
window_condition = df["event_time"].between(
-PRE_WINDOW,
POST_WINDOW,
)
同时，为了保证晚处理组尚未受到政策影响，分析期必须早于晚处理年份：
pre_late_treatment = df["year"] < LATE_YEAR
4. 原始对照组
第一种设定使用从未接受政策的地区作为对照组。
original_sample = df.loc[
window_condition
& pre_late_treatment
& (
df["first_treat"].eq(EARLY_YEAR)
| df["first_treat"].isna()
)
].copy()
检查处理组和对照组数量：
original_group_count = (
original_sample[
["county_id", "early_group"]
]
.drop_duplicates()
["early_group"]
.value_counts()
)
print("原始对照组设定：")
print(original_group_count)
其中：
估计事件研究模型：
formula = """
outcome
~ i(event_time, early_group, ref=-1)
+ control1
+ control2
+ control3
| county_id
+ year
"""
fit_never = pf.feols(
formula,
data=original_sample,
vcov={"CRV1": "county_id"},
)
print(fit_never.summary())
模型控制了：
PyFixest的feols()可以通过公式中| county_id + year吸收县固定效应和年份固定效应，并通过vcov={"CRV1": "county_id"}计算县级聚类标准误。
绘制平行趋势图：
fit_never.iplot(
coord_flip=False,
title="从未处理组作为对照组",
yintercept=0,
)
5. 晚处理对照组
第二种设定仅保留早处理县和晚处理县。
late_control_sample = df.loc[
window_condition
& pre_late_treatment
& df["first_treat"].isin(
[EARLY_YEAR, LATE_YEAR]
)
].copy()
检查样本数量：
late_group_count = (
late_control_sample[
["county_id", "first_treat"]
]
.drop_duplicates()
["first_treat"]
.value_counts()
.sort_index()
)
print("晚处理对照组设定：")
print(late_group_count)
在这个样本中：
重新估计事件研究模型：
fit_late = pf.feols(
formula,
data=late_control_sample,
vcov={"CRV1": "county_id"},
)
print(fit_late.summary())
绘制重新定义对照组后的事件研究图：
fit_late.iplot(
coord_flip=False,
title="晚处理组作为对照组",
yintercept=0,
)
6. 提取事件研究结果
为了统一输出结果，可以编写一个系数提取函数。
def extract_event_results(
model,
model_name: str,
) -> pd.DataFrame:
"""
提取PyFixest事件研究系数。
"""
result = model.tidy().reset_index()
term_column = result.columns[0]
result["event_time"] = (
result[term_column]
.astype(str)
.str.extract(
r"event_time::(-?\d+(?:\.\d+)?)",
expand=False,
)
)
result = result.loc[
result["event_time"].notna()
].copy()
result["event_time"] = (
result["event_time"].astype(float)
)
result = result.rename(
columns={
"Estimate": "estimate",
"Std. Error": "std_error",
"Pr(>|t|)": "p_value",
"2.5%": "ci_low",
"97.5%": "ci_high",
}
)
keep_columns = [
"event_time",
"estimate",
"std_error",
"p_value",
"ci_low",
"ci_high",
]
result = result[keep_columns]
result["model"] = model_name
# 补回基准期
base_period = pd.DataFrame(
{
"event_time": [-1.0],
"estimate": [0.0],
"std_error": [0.0],
"p_value": [np.nan],
"ci_low": [0.0],
"ci_high": [0.0],
"model": [model_name],
}
)
result = pd.concat(
[result, base_period],
ignore_index=True,
)
return result.sort_values(
"event_time"
).reset_index(drop=True)
提取两组模型结果：
result_never = extract_event_results(
fit_never,
"从未处理组",
)
result_late = extract_event_results(
fit_late,
"晚处理组",
)
event_results = pd.concat(
[result_never, result_late],
ignore_index=True,
)
event_results.to_excel(
"event_study_results.xlsx",
index=False,
)
print(event_results)
7. 绘制左右对比图
下面将两种对照组的平行趋势图放在同一张图片中。
def plot_event_study(
ax,
result: pd.DataFrame,
title: str,
) -> None:
"""
绘制单个事件研究图。
"""
ax.errorbar(
result["event_time"],
result["estimate"],
yerr=[
result["estimate"] - result["ci_low"],
result["ci_high"] - result["estimate"],
],
fmt="o-",
capsize=4,
)
ax.axhline(
y=0,
linewidth=1,
)
ax.axvline(
x=-0.5,
linestyle="--",
linewidth=1,
)
ax.set_title(title)
ax.set_xlabel("相对政策实施时间")
ax.set_ylabel("动态政策效应")
ax.set_xticks(
sorted(result["event_time"].unique())
)
fig, axes = plt.subplots(
1,
2,
figsize=(13, 5),
sharey=True,
)
plot_event_study(
axes[0],
result_never,
"从未处理组作为对照组",
)
plot_event_study(
axes[1],
result_late,
"晚处理组作为对照组",
)
plt.tight_layout()
plt.savefig(
"parallel_trend_comparison.png",
dpi=300,
bbox_inches="tight",
)
plt.show()
左图展示使用从未处理地区作为对照组的结果，右图展示使用晚处理地区作为对照组的结果。
8. 计算事前系数指标
除了观察图形，还可以计算政策实施前系数的平均绝对值和最大绝对值。
def pretrend_metrics(
result: pd.DataFrame,
) -> pd.Series:
"""
计算事前趋势的简单诊断指标。
"""
pre_result = result.loc[
result["event_time"] <= -2
].copy()
return pd.Series(
{
"事前平均绝对系数": (
pre_result["estimate"]
.abs()
.mean()
),
"事前最大绝对系数": (
pre_result["estimate"]
.abs()
.max()
),
"事前显著系数数量": (
pre_result["p_value"] < 0.05
).sum(),
"事前系数数量": len(pre_result),
}
)
比较两种对照组：
pretrend_comparison = pd.DataFrame(
{
"从未处理组": pretrend_metrics(
result_never
),
"晚处理组": pretrend_metrics(
result_late
),
}
).T
print(pretrend_comparison)
pretrend_comparison.to_excel(
"pretrend_comparison.xlsx"
)
需要注意，平均绝对系数和显著系数数量只是辅助诊断。
正式研究中还应当：
9. 分期DID扩展
如果样本中存在多个政策处理批次，并且不同批次的政策效应可能不同，可以进一步使用适合分期处理的事件研究估计。
PyFixest的event_study()支持传统双向固定效应模型和完全交互的饱和事件研究模型；其中saturated设定可以估计批次与事件时间相互作用的动态效应。
首先将从未处理单位的处理年份编码为零：
df["treat_cohort"] = (
df["first_treat"]
.fillna(0)
.astype(int)
)
估计传统TWFE事件研究：
fit_twfe = pf.event_study(
data=df,
yname="outcome",
idname="county_id",
tname="year",
gname="treat_cohort",
estimator="twfe",
)
fit_twfe.iplot()
估计完全交互事件研究：
fit_saturated = pf.event_study(
data=df,
yname="outcome",
idname="county_id",
tname="year",
gname="treat_cohort",
estimator="saturated",
)
fit_saturated.iplot()
按事件时间汇总不同处理批次的估计结果：
dynamic_att = fit_saturated.aggregate(
weighting="shares"
)
print(dynamic_att)
绘制聚合后的动态政策效应：
fit_saturated.iplot_aggregate(
weighting="shares"
)
官方文档还提供了did2s()和lpdid()等分期DID估计接口，可用于处理分期实施与动态效应问题。
结果展示
以下结果用于展示推文中的输出形式，实际系数需要以具体数据运行结果为准。
原始对照组
首先使用从未接受政策的地区作为对照组。
政策实施前第5年至第2年的估计系数均显著为负。
这说明，在政策正式实施之前，早处理县与从未处理县之间已经存在系统性的趋势差异。
此时，政策后的正向系数可能同时包含：
因此，从未处理地区可能不是一个理想的反事实对照组。
晚处理对照组
将对照组替换为2020年才进入政策的县后，得到以下结果。
重新定义对照组后，政策实施前的系数明显向零靠近，而且均不再显著。
这说明早处理县与晚处理县在政策实施前具有更加接近的发展趋势。
结果对比
可以看到，更换对照组后，政策实施后的估计方向没有发生根本变化，但政策实施前的趋势差异明显缩小。
这表明原始平行趋势不成立，很可能不是因为控制变量不足，而是因为从未处理县与早处理县缺乏可比性。
结果解释
重新定义对照组的目的，并不是为了让平行趋势检验机械地“通过”。
真正需要回答的是：
哪些单位能够更加合理地代表处理组在没有接受政策时的发展轨迹？
如果政策进入具有明显的选择性，那么所有未处理地区并不一定都适合作为对照组。
晚处理县与早处理县都会在未来进入同一政策，说明它们通常满足相近的政策资格或筛选条件，只是接受政策的时间存在差异。
因此，在晚处理县正式接受政策之前，将其作为早处理县的对照组，往往能够提高处理组与对照组的可比性。
不过，这种方法仍然需要注意以下问题。
预期效应
如果晚处理地区提前知道自己即将进入政策，可能会在正式实施前调整投资、财政支出或就业安排。
此时，晚处理组在正式处理前就可能受到政策影响。
样本变化
从“所有从未处理县”改为“晚处理县”后，对照组规模通常会下降。
因此，需要同时报告：
处理效应异质性
当政策分批实施，并且不同批次、不同年份的政策效果存在差异时，传统双向固定效应模型可能无法给出容易解释的平均处理效应。
此时应进一步使用：
研究设计优先
合理的分析顺序应当是：
平行趋势不通过时，不应首先通过反复调整参数寻找一张“更好看”的图。
DID实证分析的核心不是让平行趋势检验通过，而是构造一个可信的反事实。
