R plm 与 Fixst 包 - 结果不同？

我试图理解为什么 R 包”plm" and "fixest“当我使用异方差稳健标准误差（“HC1”）和状态固定效应估计面板模型时，给我不同的标准误差。

有人给我提示吗？

这是代码：

library(AER)       # For the Fatality Dataset
library(plm)       # PLM 
library(fixest)    # Fixest
library(tidyverse) # Data Management 


data("Fatalities")

# Create new variable : fatality rate
Fatalities <- Fatalities %>% 
  mutate(fatality_rate = (fatal/pop)*10000)

# Estimate Fixed Effects model using the plm package
plm_reg <- plm(fatality_rate ~ beertax,
               data = Fatalities,
               index = c("state", "year"),
               effect = "individual")

# Print Table with adjusted standard errors
coeftest(plm_reg, vcov. = vcovHC, type = "HC1")
# Output 
>t test of coefficients:

        Estimate Std. Error t value Pr(>|t|)  
beertax -0.65587    0.28880  -2.271  0.02388 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Estimate the very same model using the fixest package
# fixest is much faster and user friendly (in my opinion) 
fixest_reg <- feols(fatality_rate ~ beertax | state , 
                    data = Fatalities,
                    vcov = "HC1",
                    panel.id = ~ state + year)

# print table 
etable(fixest_reg)

#output
>                         fixest_reg
Dependent Var.:      fatality_rate
                                  
beertax         -0.6559** (0.2033)
Fixed-Effects:  ------------------
state                          Yes
_______________ __________________
S.E. type       Heteroskedas.-rob.
Observations                   336
R2                         0.90501
Within R2                  0.04075

在此示例中，使用时标准误差较大plm与固定结果相比（如果state+year使用固定效应）。有谁知道发生这种情况的原因？

实际上 VCOV 是不同的。

In plm vcovHC默认为 Arellano (1987)，它也考虑了序列相关性。查看文档here http://vps58738.lws-hosting.com/Rpkg/plm/reference/vcovHC.plm.html.

如果添加参数method = "white1"，您最终会得到相同类型的 VCOV。

最后，您还需要更改固定效应的计算方式fixest获得相同的标准误差（参见小样本校正的详细信息here https://lrberge.github.io/fixest/articles/standard_errors.html).

结果如下：

# Requesting "White" VCOV
coeftest(plm_reg, vcov. = vcovHC, type = "HC1", method = "white1")
#> 
#> t test of coefficients:
#> 
#>         Estimate Std. Error t value  Pr(>|t|)    
#> beertax -0.65587    0.18815 -3.4858 0.0005673 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
# Changing the small sample correction in fixest (discarding the fixed-effects)
etable(fixest_reg, vcov = list("hc1", hc1 ~ ssc(fixef.K = "none")), fitstat = NA)
#>                         fixest_reg          fixest_reg
#> Dependent Var.:      fatality_rate       fatality_rate
#>                                                       
#> beertax         -0.6559** (0.2033) -0.6559*** (0.1882)
#> Fixed-Effects:  ------------------ -------------------
#> state                          Yes                 Yes
#> _______________ __________________ ___________________
#> S.E. type       Heteroskedas.-rob. Heteroskedast.-rob.

# Final comparison
rbind(se(vcovHC(plm_reg, type = "HC1", method = "white1")),
      se(fixest_reg, hc1 ~ ssc(fixef.K = "none")))
#>        beertax
#> [1,] 0.1881536
#> [2,] 0.1881536