使用 R 的事件研究设计中具有交互项的公式
Formula with interaction terms in event-study designs using R
我正在估计 R 中通常称为“事件研究”的差异模型的规范。基本上，我们随着时间的推移观察处理和控制单元，并估计带有参数的双向固定效应模型对于每个时间段接受治疗的“效果”（省略一个时期，通常是治疗前的一个时期，作为参考时期）。我正在努力研究如何使用 R 公式紧凑地指定此模型。
例如，这里是模型...
library(lfe)
library(tidyverse)
library(dummies)
N <- 100
df <- tibble(
id = rep(1:N, 5),
treat = id >= ceiling(N / 2),
time = rep(1:5, each=N),
x = rnorm(5 * N)
)
# produce an outcome variable
df <- df %>% mutate(
y = x - treat * (time == 5) + time + rnorm(5*N)
)
head(df)
# easily recover the parameters with the true model...
summary(felm(
y ~ x + I(treat * (time == 5)) | id + time, data = df
))
现在，我想使用第 4 期作为基线进行事件研究设计，因为治疗发生在第 5 期。我们预计前期 (1–4) 的系数接近零，并且治疗效果为负对于治疗期间的治疗 (time == 5)
df$timefac <- factor(df$time, levels = c(4, 1, 2, 3, 5))
summary(felm(
y ~ x + treat * timefac | id + time, data = df
))
这看起来不错，但会产生很多 NA，因为一些系数被单位和时间效应吸收了。理想情况下，我可以指定没有这些系数的模型...
# create dummy for each time period for treated units
tdum <- dummy(df$time)
df <- bind_cols(df, as.data.frame(tdum))
df <- df %>% mutate_at(vars(time1:time5), ~ . * treat)
# estimate model, manually omitting one dummy
summary(felm(
y ~ x + time1 + time2 + time3 + time5 | id + time, data = df
))
现在，问题是如何以紧凑的方式指定此模型。我认为下面的方法会起作用，但它会产生非常不可预测的输出...
summary(felm(
y ~ x + treat:timefac | id + time, data = df
))
综上所述，R 不使用周期 4 作为参考周期，有时选择包括与未处理而不是处理的交互。输出是...
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x 0.97198 0.05113 19.009 < 2e-16 ***
treatFALSE:timefac4 NA NA NA NA
treatTRUE:timefac4 -0.19607 0.28410 -0.690 0.49051
treatFALSE:timefac1 NA NA NA NA
treatTRUE:timefac1 -0.07690 0.28572 -0.269 0.78796
treatFALSE:timefac2 NA NA NA NA
treatTRUE:timefac2 NA NA NA NA
treatFALSE:timefac3 0.15525 0.28482 0.545 0.58601
treatTRUE:timefac3 NA NA NA NA
treatFALSE:timefac5 0.97340 0.28420 3.425 0.00068 ***
treatTRUE:timefac5 NA NA NA NA
有没有一种方法可以指定此模型，而无需为每个时间段的处理单元手动生成虚拟变量和交互项？
如果你了解 Stata，我基本上是在寻找像这样简单的东西：
areg y x i.treat##ib4.time, absorb(id)
（请注意告诉 Stata 将变量视为分类变量是多么简单——i 前缀——无需为时间制作虚拟变量，并且还指出周期 4 应该是基期——b4前缀。）
您可以重新定义 timefac，以便将未处理的观察值编码为省略的时间类别。
df %>%
mutate(time = ifelse(treat == 0, 4, time),
timefac = factor(time, levels = c(4, 1, 2, 3, 5)))
然后，您可以在没有交互的情况下使用 timefac 并获得没有 NA 的回归 table。
summary(felm(
y ~ x + timefac | id + time, data = df
))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x 0.98548 0.05028 19.599 < 2e-16 ***
time_fac1 -0.01335 0.27553 -0.048 0.961
time_fac2 -0.10332 0.27661 -0.374 0.709
time_fac3 0.24169 0.27575 0.876 0.381
time_fac5 -1.16305 0.27557 -4.221 3.03e-05 ***
这个想法来自：https://blogs.worldbank.org/impactevaluations/econometrics-sandbox-event-study-designs-co
包 fixest 执行固定效应估计（如 lfe）并包含处理交互的实用程序。函数 i（或 interact）就是您要找的。
这里是一个例子，其中治疗与第 5 年相互作用，第 5 年被淘汰：
library(fixest)
data(base_did)
est_did = feols(y ~ x1 + i(treat, period, 5) | id + period, base_did)
est_did
#> OLS estimation, Dep. Var.: y
#> Observations: 1,080
#> Fixed-effects: id: 108, period: 10
#> Standard-errors: Clustered (id)
#> Estimate Std. Error t value Pr(>|t|)
#> x1 0.973490 0.045678 21.312000 < 2.2e-16 ***
#> treat:period::1 -1.403000 1.110300 -1.263700 0.206646
#> treat:period::2 -1.247500 1.093100 -1.141200 0.254068
#> treat:period::3 -0.273206 1.106900 -0.246813 0.805106
#> treat:period::4 -1.795700 1.088000 -1.650500 0.099166 .
#> treat:period::6 0.784452 1.028400 0.762798 0.445773
#> treat:period::7 3.598900 1.101600 3.267100 0.001125 **
#> treat:period::8 3.811800 1.247500 3.055500 0.002309 **
#> treat:period::9 4.731400 1.097100 4.312600 1.8e-05 ***
#> treat:period::10 6.606200 1.120500 5.895800 5.17e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -2,984.58 Adj. R2: 0.48783
好的一点是，您可以从估计中绘制出交互系数，以快速直观地表示结果（如果您觉得图表过于简洁，不用担心，您可以自定义其中的几乎所有内容）：
coefplot(est_did)
如果您不想使用 fixest 进行估算，您仍然可以使用函数 i 来创建交互。它的语法是 i(var, f, ref, drop, keep)：它将变量 var 与 f 中每个值的虚拟变量交互。您可以 select 使用参数 ref、drop 和 keep 保留 f 的哪些值。 drop 嗯...从 f 和 ref 中删除值与 drop 相同，但引用显示在 coefplot 中（而值在drop 没有出现在图表中。
下面是 i 的一个例子：
head(with(base_did, i(treat, period, keep = 3:7)))
#> treat:period::3 treat:period::4 treat:period::5 treat:period::6 treat:period::7
#> 1 0 0 0 0 0
#> 2 0 0 0 0 0
#> 3 1 0 0 0 0
#> 4 0 1 0 0 0
#> 5 0 0 1 0 0
#> 6 0 0 0 1 0
head(with(base_did, i(treat, period, drop = 3:7)))
#> treat:period::1 treat:period::2 treat:period::8 treat:period::9 treat:period::10
#> 1 1 0 0 0 0
#> 2 0 1 0 0 0
#> 3 0 0 0 0 0
#> 4 0 0 0 0 0
#> 5 0 0 0 0 0
#> 6 0 0 0 0 0
您可以在 fixest here.
上找到更多信息
我正在估计 R 中通常称为“事件研究”的差异模型的规范。基本上，我们随着时间的推移观察处理和控制单元，并估计带有参数的双向固定效应模型对于每个时间段接受治疗的“效果”（省略一个时期，通常是治疗前的一个时期，作为参考时期）。我正在努力研究如何使用 R 公式紧凑地指定此模型。
例如，这里是模型...
library(lfe)
library(tidyverse)
library(dummies)
N <- 100
df <- tibble(
id = rep(1:N, 5),
treat = id >= ceiling(N / 2),
time = rep(1:5, each=N),
x = rnorm(5 * N)
)
# produce an outcome variable
df <- df %>% mutate(
y = x - treat * (time == 5) + time + rnorm(5*N)
)
head(df)
# easily recover the parameters with the true model...
summary(felm(
y ~ x + I(treat * (time == 5)) | id + time, data = df
))
现在，我想使用第 4 期作为基线进行事件研究设计，因为治疗发生在第 5 期。我们预计前期 (1–4) 的系数接近零，并且治疗效果为负对于治疗期间的治疗 (time == 5)
df$timefac <- factor(df$time, levels = c(4, 1, 2, 3, 5))
summary(felm(
y ~ x + treat * timefac | id + time, data = df
))
这看起来不错，但会产生很多 NA，因为一些系数被单位和时间效应吸收了。理想情况下，我可以指定没有这些系数的模型...
# create dummy for each time period for treated units
tdum <- dummy(df$time)
df <- bind_cols(df, as.data.frame(tdum))
df <- df %>% mutate_at(vars(time1:time5), ~ . * treat)
# estimate model, manually omitting one dummy
summary(felm(
y ~ x + time1 + time2 + time3 + time5 | id + time, data = df
))
现在，问题是如何以紧凑的方式指定此模型。我认为下面的方法会起作用，但它会产生非常不可预测的输出...
summary(felm(
y ~ x + treat:timefac | id + time, data = df
))
综上所述，R 不使用周期 4 作为参考周期，有时选择包括与未处理而不是处理的交互。输出是...
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x 0.97198 0.05113 19.009 < 2e-16 ***
treatFALSE:timefac4 NA NA NA NA
treatTRUE:timefac4 -0.19607 0.28410 -0.690 0.49051
treatFALSE:timefac1 NA NA NA NA
treatTRUE:timefac1 -0.07690 0.28572 -0.269 0.78796
treatFALSE:timefac2 NA NA NA NA
treatTRUE:timefac2 NA NA NA NA
treatFALSE:timefac3 0.15525 0.28482 0.545 0.58601
treatTRUE:timefac3 NA NA NA NA
treatFALSE:timefac5 0.97340 0.28420 3.425 0.00068 ***
treatTRUE:timefac5 NA NA NA NA
有没有一种方法可以指定此模型，而无需为每个时间段的处理单元手动生成虚拟变量和交互项？
如果你了解 Stata，我基本上是在寻找像这样简单的东西：
areg y x i.treat##ib4.time, absorb(id)
（请注意告诉 Stata 将变量视为分类变量是多么简单——i 前缀——无需为时间制作虚拟变量，并且还指出周期 4 应该是基期——b4前缀。）
您可以重新定义 timefac，以便将未处理的观察值编码为省略的时间类别。
df %>%
mutate(time = ifelse(treat == 0, 4, time),
timefac = factor(time, levels = c(4, 1, 2, 3, 5)))
然后，您可以在没有交互的情况下使用 timefac 并获得没有 NA 的回归 table。
summary(felm(
y ~ x + timefac | id + time, data = df
))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x 0.98548 0.05028 19.599 < 2e-16 ***
time_fac1 -0.01335 0.27553 -0.048 0.961
time_fac2 -0.10332 0.27661 -0.374 0.709
time_fac3 0.24169 0.27575 0.876 0.381
time_fac5 -1.16305 0.27557 -4.221 3.03e-05 ***
这个想法来自：https://blogs.worldbank.org/impactevaluations/econometrics-sandbox-event-study-designs-co
包 fixest 执行固定效应估计（如 lfe）并包含处理交互的实用程序。函数 i（或 interact）就是您要找的。
这里是一个例子，其中治疗与第 5 年相互作用，第 5 年被淘汰：
library(fixest)
data(base_did)
est_did = feols(y ~ x1 + i(treat, period, 5) | id + period, base_did)
est_did
#> OLS estimation, Dep. Var.: y
#> Observations: 1,080
#> Fixed-effects: id: 108, period: 10
#> Standard-errors: Clustered (id)
#> Estimate Std. Error t value Pr(>|t|)
#> x1 0.973490 0.045678 21.312000 < 2.2e-16 ***
#> treat:period::1 -1.403000 1.110300 -1.263700 0.206646
#> treat:period::2 -1.247500 1.093100 -1.141200 0.254068
#> treat:period::3 -0.273206 1.106900 -0.246813 0.805106
#> treat:period::4 -1.795700 1.088000 -1.650500 0.099166 .
#> treat:period::6 0.784452 1.028400 0.762798 0.445773
#> treat:period::7 3.598900 1.101600 3.267100 0.001125 **
#> treat:period::8 3.811800 1.247500 3.055500 0.002309 **
#> treat:period::9 4.731400 1.097100 4.312600 1.8e-05 ***
#> treat:period::10 6.606200 1.120500 5.895800 5.17e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Log-likelihood: -2,984.58 Adj. R2: 0.48783
好的一点是，您可以从估计中绘制出交互系数，以快速直观地表示结果（如果您觉得图表过于简洁，不用担心，您可以自定义其中的几乎所有内容）：
coefplot(est_did)
如果您不想使用 fixest 进行估算，您仍然可以使用函数 i 来创建交互。它的语法是 i(var, f, ref, drop, keep)：它将变量 var 与 f 中每个值的虚拟变量交互。您可以 select 使用参数 ref、drop 和 keep 保留 f 的哪些值。 drop 嗯...从 f 和 ref 中删除值与 drop 相同，但引用显示在 coefplot 中（而值在drop 没有出现在图表中。
下面是 i 的一个例子：
head(with(base_did, i(treat, period, keep = 3:7)))
#> treat:period::3 treat:period::4 treat:period::5 treat:period::6 treat:period::7
#> 1 0 0 0 0 0
#> 2 0 0 0 0 0
#> 3 1 0 0 0 0
#> 4 0 1 0 0 0
#> 5 0 0 1 0 0
#> 6 0 0 0 1 0
head(with(base_did, i(treat, period, drop = 3:7)))
#> treat:period::1 treat:period::2 treat:period::8 treat:period::9 treat:period::10
#> 1 1 0 0 0 0
#> 2 0 1 0 0 0
#> 3 0 0 0 0 0
#> 4 0 0 0 0 0
#> 5 0 0 0 0 0
#> 6 0 0 0 0 0
您可以在 fixest here.