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

【问题标题】：Formula with interaction terms in event-study designs using R使用 R 的事件研究设计中具有交互项的公式
【发布时间】：2020-11-02 23:57:35
【问题描述】：

我正在估计 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
))

这看起来不错，但会产生很多NAs，因为有几个系数被单位和时间效应吸收了。理想情况下，我可以指定没有这些系数的模型...

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

【问题讨论】：

标签： r formula dummy-variable

【解决方案1】：

包fixest 执行固定效应估计（如lfe）并包含处理交互的实用程序。函数i（或interact）就是你要找的。p>

这是一个例子，其中治疗与年份交互，第 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 中的每个值的虚拟变量交互。您可以使用参数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

您可以在fixesthere找到更多信息。

【讨论】：

这看起来很棒！我喜欢明确指定基期和相互作用（在固定效应和感兴趣的系数中）。

【解决方案2】：

您可以重新定义 timefac，以便将未处理的观测值编码为省略的时间类别。

df %>% 
  mutate(time = ifelse(treat == 0, 4, time),
         timefac = factor(time, levels = c(4, 1, 2, 3, 5)))

然后，您可以在没有交互的情况下使用 timefac，并获得一个没有 NA 的回归表。

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

【讨论】：