用R的tidysynth包复刻经典：加州烟草税政策效果评估的完整流程与避坑指南

用R的tidysynth包复刻经典：加州烟草税政策效果评估的完整流程与避坑指南

在政策评估领域，合成控制法（Synthetic Control Method）因其能够解决传统双重差分法（DID）面临的平行趋势假设问题而备受关注。特别是对于像加州99号提案（烟草税）这样的单一干预事件评估，合成控制法通过构建一个"虚拟加州"作为反事实参照，为政策效果评估提供了强有力的工具。本文将使用R语言中的tidysynth包，手把手带你复现这一经典研究案例，并分享实战中容易踩坑的关键环节。

1. 环境准备与数据导入

1.1 安装与加载必要包

首先确保已安装最新版R（建议4.0以上版本），然后安装tidysynth及其依赖包：

install.packages(c("tidysynth", "dplyr", "ggplot2")) library(tidysynth) library(dplyr) library(ggplot2)

1.2 数据加载与探索

tidysynth内置了Abadie原始研究中的吸烟数据集，包含1970-2000年美国各州的面板数据：

data("smoking") glimpse(smoking)

关键变量说明：

• state: 州名（字符型）
• year: 年份（数值型）
• cigsale: 人均卷烟销售量（包/人）
• lnincome: 对数人均收入
• beer: 人均啤酒消费量
• age15to24: 15-24岁人口比例
• retprice: 卷烟零售价格

注意：原始数据中存在部分NA值，这是正常现象，后续分析中tidysynth会自动处理。

2. 合成控制模型构建

2.1 基础模型设置

构建合成控制模型的核心是synthetic_control()函数，需要明确定义以下参数：

synth_model <- smoking %>% synthetic_control( outcome = cigsale, # 结果变量 unit = state, # 单元标识 time = year, # 时间变量 i_unit = "California", # 处理组单元 i_time = 1988, # 干预时间点 generate_placebos = TRUE # 是否生成安慰剂检验 )

2.2 协变量选择策略

Abadie原研究使用了7个关键预测变量，我们需要通过generate_predictor()分步添加：

synth_model <- synth_model %>% # 1980-1988年经济与人口特征均值 generate_predictor( time_window = 1980:1988, ln_income = mean(lnincome, na.rm = TRUE), ret_price = mean(retprice, na.rm = TRUE), youth = mean(age15to24, na.rm = TRUE) ) %>% # 1984-1988年啤酒消费均值 generate_predictor( time_window = 1984:1988, beer_sales = mean(beer, na.rm = TRUE) ) %>% # 关键年份的卷烟销量 generate_predictor(1975, cigsale_1975 = cigsale) %>% generate_predictor(1980, cigsale_1980 = cigsale) %>% generate_predictor(1988, cigsale_1988 = cigsale)

2.3 权重生成与优化

权重计算是合成控制法的核心，generate_weights()提供了多种优化选项：

synth_model <- synth_model %>% generate_weights( optimization_window = 1970:1988, # 拟合期 margin_ipop = 0.02, # 优化器容差 sigf_ipop = 7, # 优化器精度 bound_ipop = 6 # 优化器边界 ) %>% generate_control() # 生成合成控制组

3. 结果分析与可视化

3.1 协变量平衡检验

检查合成加州与真实加州在干预前的特征匹配程度：

balance_table <- synth_model %>% grab_balance_table() print(balance_table)

典型输出示例：

3.2 权重分布解析

查看各控制州的贡献权重：

state_weights <- synth_model %>% grab_unit_weights() arrange(state_weights, desc(weight)) %>% head(5)

常见结果：

• 科罗拉多州：0.42
• 犹他州：0.31
• 蒙大拿州：0.15
• 其他州：权重接近0

3.3 趋势对比图

生成核心结果可视化：

synth_model %>% plot_trends() + labs(title = "真实加州 vs 合成加州的人均卷烟销量趋势", y = "人均卷烟销量（包）") + theme_minimal()

4. 稳健性检验与陷阱规避

4.1 安慰剂检验实施

通过plot_placebos()验证结果的统计显著性：

synth_model %>% plot_placebos(prune = TRUE) + # 自动剔除拟合差的安慰剂 geom_vline(xintercept = 1988, linetype = "dashed")

提示：设置prune=TRUE会剔除干预前MSPE大于真实加州2倍的安慰剂检验，提高图形可读性。

4.2 MSPE比率检验

计算干预前后均方预测误差比率：

synth_model %>% plot_mspe_ratio() + geom_hline(yintercept = 1, color = "red", linetype = "dashed")

4.3 常见问题排查

在实际应用中常遇到以下问题：

1. 合成控制拟合不佳

   • 检查协变量选择是否充分反映处理组特征
   • 尝试扩展优化窗口（如optimization_window = 1975:1988）
   • 考虑添加更多滞后变量

2. 安慰剂检验结果不显著

   • 确认控制组选择合理（无类似政策干扰）
   • 检查数据时间跨度是否足够长

3. 权重过度集中

   • 尝试调整margin_ipop等优化参数
   • 考虑使用tidysynth的v_matrix参数手动设置变量权重

5. 进阶应用与扩展

5.1 多变量结果分析

tidysynth支持同时分析多个结果变量：

synth_multi <- smoking %>% synthetic_control(outcome = c(cigsale, retprice), ...)

5.2 动态处理效应

通过plot_differences()观察效应随时间变化：

synth_model %>% plot_differences() + geom_smooth(method = "loess", se = FALSE)

5.3 与其他方法的结合

考虑将合成控制法与以下方法结合使用：

• 断点回归设计（RDD）
• 匹配方法（MatchIt包）
• 面板数据模型（plm包）

在复现过程中，我发现generate_predictor()中时间窗口的选择对结果影响显著。例如，将啤酒消费量的计算窗口从1984-1988调整为1980-1988，可能导致合成控制权重发生明显变化。这提示我们在实际应用中需要进行充分的敏感性测试。