基于具有聚类 SE 的模型的图

Question

我正在尝试基于具有聚集标准误差的回归模型创建图。根据数据集的大小，我通常使用来自"lfe" 包的felm 或来自"multiwayvcov" 和"lmtest" 包的cluster.vcov 和coeftest 对标准错误进行聚类。

以下是使用这两种方法可能会出现的回归示例：

# Load Packages
library('lfe')
library('multiwayvcov')
library('lmtest')
data(mtcars)

# Using felm
m1 <- felm(mpg ~ wt + disp + wt*disp | 0 | 0 | carb, data=mtcars)

# multiwayvcov
m2 <- glm(mpg ~ wt + disp + wt*disp, data=mtcars)
er2 <- cluster.vcov(m2, mtcars, 'carb')
m2ct <- coeftest(m2, er2)

问题在于这两个选项产生的对象类型与 R 中的大多数绘图功能不兼容。

例如，如果我试图创建这样的情节：

interplot(m=m2, var1='wt', var2='disp', ci=.9, ralpha=.3, rfill='dodgerblue1') +
  geom_hline(yintercept=0, linetype="dashed") +
  guides(colour=guide_legend(override.aes=list(size=3)))

如何从这些具有聚集标准误差的模型中创建具有正确置信区间的图？

Answer 1

miceadds 包中有一个 glm.cluster 函数，它适合具有集群稳健标准错误的 GLM，并且似乎类似于 felmresults。

library(miceadds)
m.glm.cl <- glm.cluster(mpg ~ wt + disp + wt*disp, cluster="carb", data=mtcars)
summary(m.glm.cl)
#                Estimate  Std. Error    t value      Pr(>|t|)
# (Intercept) 44.08199770 1.960767641  22.482010 6.225590e-112
# wt          -6.49567966 0.635946018 -10.214200  1.712866e-24
# disp        -0.05635816 0.009436760  -5.972194  2.340843e-09
# wt:disp      0.01170542 0.001965829   5.954445  2.609566e-09

library('lfe')
m.felm <- felm(mpg ~ wt + disp + wt*disp | 0 | 0 | carb, data=mtcars)
summary(m.felm)$coef
#                Estimate Cluster s.e.    t value     Pr(>|t|)
# (Intercept) 44.08199770  1.960767641  22.482010 1.846712e-19
# wt          -6.49567966  0.635946018 -10.214200 6.015931e-11
# disp        -0.05635816  0.009436760  -5.972194 1.972266e-06
# wt:disp      0.01170542  0.001965829   5.954445 2.068825e-06

p 值不同，但是，我不知道它们是否/如何在 interplot 中使用（您可能可以通过查看 interplot 的方法来弄清楚这一点使用：interplot::interplot.default)。

无论如何，glm.cluster 返回一个列表，其中glm_res 是我们正在寻找的，

names(m.glm.cl)
# [1] "glm_res" "vcov"

以及我们可以为interplot 提供什么。

library("interplot")
interplot(m=m.glm.cl$glm_res, var1='wt', var2='disp', ci=.9, ralpha=.3, rfill='dodgerblue1') +
  geom_hline(yintercept=0, linetype="dashed") +
  guides(colour=guide_legend(override.aes=list(size=3)))