子样本之间的系数检验比较

Question

我正在尝试对两个子样本的系数检验进行比较。为了实现这一目标，我执行以下操作：

full_model <- lm(y ~ v1*subsample_dummy + fixed_effects, data=df)
reduced_model <- lm(y ~ v1 + subsample_dummy + fixed_effects, data=df)

test <- anova(full_model, reduced_model)

上面给了我结果。

但是，在必须通过 year 变量对模型进行

聚类

的情况下，我不确定如何执行相同的操作。

我可以使用以下代码对 lm 模型进行聚类：

library(sandwich)
# cluster by year
clustered_se <- vcovCL(full_model, ~ year) 
clustered_se1 <- vcovCL(reduced_model, ~ year)

# generate summaries with clustered standard errors 
a <- coeftest(full_model, vcov. = clustered_se) 
b <- coeftest(reduced_model, vcov. = clustered_se1) 

但是，问题仍然存在，因为我仍然无法做到：

anova(a, b)

模型需要标准误聚类时，如何实现跨子样本的系数检验比较？

Answer 1

我们可以使用

sandwich::vcovCL

来获得与

lfe::felm

基本相同的标准错误。让我们来估计一些模型。

> est1 <- lfe::felm(y ~ x1 + x2 | id + firm | 0 | firm, data=d)  ## id + firm FE, clustered by firm
> est2 <- lfe::felm(y ~ x1 | id + firm | 0 | firm, data=d)  ## same, restricted model
> est3 <- lm(y ~ x1 + x2 + as.factor(id) + as.factor(firm), data=d)  ## as est1 w/o clustering
> est4 <- lm(y ~ x1 + as.factor(id) + as.factor(firm), data=d)  ## as restricted est3

比较 est3 与 est1,

的标准误差

> lmtest::coeftest(est3, vcov.=\(x) 
+                  sandwich::vcovCL(x, cluster=d$firm, type='HC0'))[1:3, ]
             Estimate Std. Error  t value      Pr(>|t|)
(Intercept) 3.7416642 0.10554929 35.44945 5.454680e-177
x1          1.0432612 0.02965723 35.17730 3.631800e-175
x2          0.4904104 0.03186679 15.38939  5.822822e-48
> coef(summary(est1))
    Estimate Cluster s.e.  t value     Pr(>|t|)
x1 1.0432612   0.02968696 35.14207 1.799173e-13
x2 0.4904104   0.03189874 15.37398 2.931632e-09

产量基本相同。

因此，根据交叉验证[1]上的帖子，我们可以使用lmtest::waldtest来比较模型est1和

est2

（还有一个工作方式不同的

lfe::waldtest

）。

> lmtest::waldtest(est3, est4, vcov=\(x) 
+                  sandwich::vcovCL(x, cluster=d$firm, type='HC0'))
Wald test

Model 1: y ~ x1 + x2 + as.factor(id) + as.factor(firm)
Model 2: y ~ x1 + as.factor(id) + as.factor(firm)
  Res.Df Df      F    Pr(>F)    
1    966                        
2    967 -1 236.83 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

希望这能让您更进一步。

顺便说一句：您绝对可以就此问题提出功能请求，作为作者GitHub的问题。

---

数据：

> set.seed(42)
> n <- 1e3
> d <- data.frame(x1=rnorm(n), x2=rnorm(n), id=factor(sample(20, n, replace=TRUE)),
+                 firm=factor(sample(13, n, replace=TRUE)), u=rnorm(n))
> id.eff <- rnorm(nlevels(d$id))
> firm.eff <- rnorm(nlevels(d$firm))
> d$y <- d$x1 + 0.5 * d$x2 + id.eff[d$id] + firm.eff[d$firm] + d$u