r - 使用 R 的事件研究设计中具有交互项的公式

Question

我正在估计 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 s 因为几个系数被单位和时间效应吸收了。理想情况下，我可以指定没有这些系数的模型......

# 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 前缀。)

Answer 1

fixest 包执行固定效应估计(如 lfe )并包含处理交互的实用程序。函数 i (或 interact )正是您要找的。
下面是一个例子，其中治疗与年份相互作用并且第 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 中每个值的虚拟变量进行交互。您可以使用参数 f 、 ref 和 drop 选择要保留的 keep 值。 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 的更多信息。