用felm (lfe)回归的三明治标准误差:动态定义公式

Question

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

下面是一个简单的例子：

# 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)

尽管回归的结果是相同的，但在第二种情况下，我不能使用三明治函数，最后一行给出了一个错误，即: Error:类型为‘符号’的对象不可设置。

有什么我能解决的线索吗？

非常感谢,

Answer 1

更新

关于调用的第二个版本不起作用的确切原因:它与model.matrix有关，它不起作用(对不起，我没有解决方案)。但是，要注意lfe和sandwich之间的兼容性，正如下面的第一个答案所强调的那样。

有两件事。首先，(我认为需要确认) felm对象似乎与sandwich变量不直接兼容，从而导致错误的结果。第二，在计算标准时涉及到的许多细节--错误，特别是关于考虑自由度的决定--这是软件之间差异的主要原因。

下面是一个示例，说明如何(几乎)使用felm获得sandwich VCOV

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模型，然后使用sandwich计算集群VCOV。最后，VCOV是成比例的，它们只是在小样本调整中有所不同。

如果您想要一个模型有几个SE，您可以查看另一个包来使用多个FEs：最固定来执行OLS和GLM估计。

使用固定

包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

有关在fixest中计算SE的详细信息，以及如何从替代软件中获取等效SE，请查看这个小片段：论标准误差。