如何获得带聚类标准误差的lm_robust()后的边际效应？

Question

我正在按年运行一个具有聚集标准错误的回归。这对于Stata来说很容易，但是我必须用R来完成，所以我使用来自lm_robust()包的estimatr函数来运行它。问题是，我现在必须得到一些变量的边际效应，但我做不到，我猜是因为集群标准错误。我遵循了lm_robust()手册上的内容，并且我看到他们只对其他函数使用了seen命令，没有集群标准错误.有没有人知道我怎样才能得到和画出边际效果？

set.seed(42)
library(fabricatr)
library(randomizr)
dat <- fabricate(
  N = 100,                        # sample size
  x = runif(N, 0, 1),             # pre-treatment covariate
  y0 = rnorm(N, mean = x),        # control potential outcome
  y1 = y0 + 0.35,                 # treatment potential outcome
  z = complete_ra(N),             # complete random assignment to treatment
  y = ifelse(z, y1, y0),          # observed outcome

  # We will also consider clustered data
  clust = sample(rep(letters[1:20], each = 5)),
  z_clust = cluster_ra(clust),
  y_clust = ifelse(z_clust, y1, y0)
)

然后，当我使用lm_robust()函数运行回归时：

library(estimatr)
lmout_cl <- lm_robust(
  y_clust ~ z_clust + x,
  data = dat,
  clusters = clust
)

最后，我试着获得利润.

library(margins)
mar_cl <- margins(lmout_cl)

但这会导致一个错误：

Error in attributes(.Data) <- c(attributes(.Data), attrib) :'names' attribute 
[1] must be the same length as the vector [0]

Answer 1

问题是，estimatr::lm_robust()产生了一个"lm_robust"对象，而margins()目前似乎不支持这个对象。我们可以使用miceadds::lm.cluster()，并获得与Stata相同的标准错误。

library(miceadds)

lmout_cl <- lm.cluster(y_clust ~ z_clust + x, data=dat, cluster=dat$clust)

这导致一个包含两个元素的列表，其中普通的lm-object存储在第一个元素中，而方差协方差矩阵的标准错误聚集在第二个元素中(参见str(lmout_cl))：

> names(lmout_cl)
[1] "lm_res" "vcov"  

margins()现在可以指定为margins(model=model, vcov=vcov)，所以我们说：

mar_cl <- with(lmout_cl, margins(lm_res, vcov=vcov))

产量

> mar_cl
Average marginal effects
stats::lm(formula = formula, data = data)

 z_clust     x
  0.6558 1.444

和

> summary(mar_cl)
  factor    AME     SE      z      p  lower  upper
       x 1.4445 0.3547 4.0728 0.0000 0.7494 2.1396
 z_clust 0.6558 0.1950 3.3633 0.0008 0.2736 1.0379

带有聚集的标准错误。

与Stata的比较

R

foreign::write.dta(dat, "dat.dta")  # export as Stata data to wd

Stata

. use dat, clear
(Written by R.              )

. quietly regress y_clust z_clust x, vce(cluster clust)

. mfx

Marginal effects after regress
      y  = Fitted values (predict)
         =  .67420391
------------------------------------------------------------------------------
variable |      dy/dx    Std. Err.     z    P>|z|  [    95% C.I.   ]      X
---------+--------------------------------------------------------------------
 z_clust*|   .6557558      .19498    3.36   0.001   .273609   1.0379        .5
       x |   1.444481      .35466    4.07   0.000   .749352  2.13961   .524479
------------------------------------------------------------------------------
(*) dy/dx is for discrete change of dummy variable from 0 to 1

. 

我们可以清楚地看到，在这样做的过程中，R产生的结果与Stata关于标准误差和边际效应的结果相同。

Answer 2

对此错误表示歉意，它阻止margins()使用estimatr版本0.10及更早版本中具有非数字集群的lm_robust()对象。这是由estimatr::lm_robust()和margins::margins()处理模型中哪些变量的内部方式创建的。

这个bug已经解决了，所以您在estimatr中有两个解决方案。

首先让我生成数据。

library(fabricatr)
library(randomizr)
dat <- fabricate(
  N = 100,
  x = runif(N),
  clust = sample(rep(letters[1:20], each = 5)),
  y_clust = rnorm(N),
  z_clust = cluster_ra(clust),
)

estimatr (v0.11.0)的最新版本

https://declaredesign.org/r/estimatr上的dev版本已经修复了这个bug，并且它将在下个月左右在CRAN上运行。

install.packages("estimatr", dependencies = TRUE,
                 repos = c("http://r.declaredesign.org", "https://cloud.r-project.org"))
library(estimatr)
lmout_cl <- lm_robust(
  y_clust ~ z_clust + x,
  data = dat,
  clusters = clust
)
library(margins)
mar_cl <- margins(lmout_cl)

estimatr (v0.10.0)的CRAN版本使用数字集群

在CRAN上使用现有版本的estimatr的一个解决办法是使用数字集群而不是字符集群。

dat <- fabricate(
  N = 100,
  x = runif(N),
  clust = sample(rep(1:20, each = 5)),
  y_clust = rnorm(N),
  z_clust = cluster_ra(clust),
)
install.packages("estimatr")
library(estimatr)
lmout_cl <- lm_robust(
  y_clust ~ z_clust + x,
  data = dat,
  clusters = clust
)
mar_cl <- margins(lmout_cl)