r - 带有 felm (lfe) 回归的三明治标准误差：动态定义公式

Question

我使用 lfe 包的 felm 函数来运行具有许多固定效果的回归，并且我需要使用 Sandwich 包来估计标准误差。回归中使用的公式是动态构建的，我将其保存在变量“表单”中。当我调用三明治函数时，它不理解将公式保存在变量中的模型。

这是一个简单的例子：

# gen process
set.seed(42)
nn = 10
n1 = 3

x <- rnorm(nn)
f1 <- sample(n1, length(x), replace=TRUE)
y <- 2.13*x + cos(f1) + rnorm(length(x), sd=0.5)

# sandwich working
est <- lfe::felm(y ~ x | f1)
summary(est)
sandwich::vcovPL(est)

# sandwich not working
form <- as.formula("y ~ x | f1")
est <- lfe::felm(form)
summary(est)
sandwich::vcovPL(est)

尽管回归的结果是相同的，但在第二种情况下，我不能使用三明治函数，最后一行给出了一个错误，内容为：错误：'symbol' 类型的对象不是子集

关于如何解决这个问题的任何线索？

非常感谢，

标签： rregression

Answer 1

更新

关于调用的第二个版本不起作用的确切原因：它与model.matrix哪个不起作用有关（对不起，我没有解决方案）。请注意 和 之间的兼容性lfe，sandwich正如我在下面的第一个答案中强调的那样。

---

有两件事。首先，（我认为但有待确认）felm对象似乎与方差不直接兼容sandwich，从而导致错误的结果。其次，计算标准误差涉及许多细节，特别是关于要考虑的自由度的决定——这是软件之间差异的主要原因。

以下是如何（几乎）获得felmVCOV的示例sandwich：

base = iris
names(base) = c("y", "x1", "x2", "x3", "species")
 
library(lfe) ; library(sandwich)
 
est_felm = felm(y ~ x1 + x2 + x3 | species | 0 | species, base)

# vcovCL does not work appropriately when applied to felm objects
vcov(est_felm) / vcovCL(est_felm, cluster = base$species)
#>           x1        x2        x3
#> x1 0.1126600 0.1011106 0.7028052
#> x2 0.1011106 0.1074147 0.1702584
#> x3 0.7028052 0.1702584 0.2571940
 
# Equivalent lm estimation
est_lm = lm(y ~ x1 + x2 + x3 + species, base)

# Now: almost equivalence. They are proportional.
vcov(est_felm) / vcovCL(est_lm, cluster = base$species)[2:4, 2:4]
#>           x1        x2        x3
#> x1 0.9863014 0.9863014 0.9863014
#> x2 0.9863014 0.9863014 0.9863014
#> x3 0.9863014 0.9863014 0.9863014

如您所见，您需要首先估计一个lm模型，然后计算集群 VCOV sandwich。您最终会得到成比例的 VCOV，它们仅在小样本调整中有所不同。

如果您想为单个模型提供多个 SE，您可以查看一个替代包，以使用多个 FE 执行 OLS 和 GLM 估计：fixst。

使用 fixst

该软件包fixest与sandwich. 下面是一个等价的例子：

library(fixest)

est_feols = feols(y ~ x1 + x2 + x3 | species, base)

# feols automatically clusters the SEs
# => proportional VCOVs
vcov(est_feols) / vcovCL(est_feols, cluster = base$species)
#>          x1       x2       x3
#> x1 1.020548 1.020548 1.020548
#> x2 1.020548 1.020548 1.020548
#> x3 1.020548 1.020548 1.020548

# Equivalences -- I:
vcov(est_feols, dof = dof(fixef.K = "none")) / vcovCL(est_feols, cluster = base$species,
                                                      type = "HC1")
#>    x1 x2 x3
#> x1  1  1  1
#> x2  1  1  1
#> x3  1  1  1 

# Equivalences -- II:
vcov(est_feols, dof = dof(adj = FALSE)) / vcovCL(est_feols, cluster = base$species)
#>    x1 x2 x3
#> x1  1  1  1
#> x2  1  1  1
#> x3  1  1  1

有关如何计算 SEfixest以及如何从替代软件获得等效 SE 的详细信息，请查看此小插图：关于标准错误。