事件研究双重差分回归方程：如何将给定Stata代码转换为R代码？

你提供的Stata命令是带分析权重、聚类标准误的事件研究DID回归，以下是对应的R代码实现，分两种常用方式：

核心对应关系 #

先明确Stata命令各部分在R中的对应：

• 被解释变量：laccidentsvso2
• 核心解释变量：weakban、strongban
• 控制变量：lpop、lunemp、permale2、lrgastax、laccidentmv2、time
• 固定效应：st1-st50（州虚拟变量）、t1-t48（时间虚拟变量）、stt1-stt50（州-时间交互虚拟变量）
• 分析权重：[aweight=pop] → R中用weight=pop
• 聚类标准误：cluster(state) → R中指定聚类分组为state

方法1：基础lm + 稳健聚类标准误 #

使用lm拟合加权回归，再通过sandwich和lmtest包计算聚类标准误：

# 加载依赖包
library(sandwich)
library(lmtest)

# 拟合加权回归（若已生成st1-st50等虚拟变量，直接列出即可；否则用factor(state)生成州固定效应）
model_lm <- lm(laccidentsvso2 ~ weakban + strongban + lpop + lunemp + permale2 + 
                 lrgastax + laccidentmv2 + factor(state) + factor(time_var) + 
                 time + factor(state):factor(time_var), 
               data = your_data, 
               weight = pop)

# 计算按state聚类的稳健标准误
cluster_se <- vcovCL(model_lm, cluster = ~ state)

# 输出带聚类标准误的结果
coeftest(model_lm, vcov = cluster_se)

方法2：fixest包（高效处理固定效应） #

fixest包专为固定效应模型设计，代码更简洁高效：

# 加载fixest包
library(fixest)

# 拟合DID模型：用|分隔控制变量和固定效应
model_fe <- feols(laccidentsvso2 ~ weakban + strongban + lpop + lunemp + permale2 + 
                    lrgastax + laccidentmv2 + time | 
                    factor(state) + factor(time_var) + factor(state):factor(time_var),
                  data = your_data, 
                  weight = pop, 
                  cluster = ~ state)

# 查看回归结果
summary(model_fe)

注意事项 #

• 替换your_data为你的实际数据框名称，time_var为原始时间变量（若已手动生成t1-t48等虚拟变量，可直接替换factor(time_var)为这些变量）
• 若Stata中stt1-stt50是自定义的州-时间交互虚拟变量，而非全量交互，需手动列出这些变量，而非用factor(state):factor(time_var)

内容的提问来源于stack exchange，提问作者user20219888