如何在 lm_robust() 后获得具有聚集标准误差的边际效应?

How to get the marginal effects after lm_robust() with clustered standard errors?

我正在 运行 进行回归分析,其中包含按年份分类的标准误差。使用 Stata 很容易做到这一点,但我必须使用 R 来做到这一点,所以我 运行 它使用 estimatr 包中的 lm_robust() 函数。问题是我现在必须得到一些变量的边际效应,但我做不到,我猜这是因为集群标准错误。我按照 lm_robust() 手册上的内容进行操作,我看到他们只将 margins 包中的 margins 命令用于其他功能而没有聚集标准错误......有没有人知道我如何获得和绘制边际效应?

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]

问题是 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 对象存储在第一个元素中,具有聚类标准误差的方差-协方差矩阵存储在第二个元素中(参见 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 产生的结果相同。

对于此错误表示歉意,该错误阻止 margins()estimatr 版本 0.10 和更早版本中使用具有非数字簇的 lm_robust() 对象。这是通过 estimatr::lm_robust()margins::margins() 处理模型中的变量的内部方式创建的。

此错误已得到解决,因此您在 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)

使用 CRAN 版本 estimatr (v0.10.0)

的数字簇

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)