问题描述
我使用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)
尽管回归的结果相同，但在第二种情况下，我无法使用三明治函数，最后一行给出了以下错误信息：错误：“符号”类型的对象不可子集化
关于如何解决此问题的任何线索吗？
非常感谢
解决方法
更新
确切说明为什么第二个版本的呼叫不起作用：它与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，则可以查看一个替代软件包，以使用多个FE执行OLS和GLM估算：fixest。
使用固定的
软件包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的详细信息，请查看以下插图：On standard-errors。