Графики на основе моделей с кластерными СЭ

Question

Я пытаюсь создать график на основе регрессионных моделей с кластерными стандартными ошибками. В зависимости от размера набора данных я обычно объединяю стандартные ошибки, используя felm из пакета "lfe" или cluster.vcov и coeftest из пакетов "multiwayvcov" и "lmtest".

Вот пример того, как могут выглядеть регрессии с этими двумя методами:

# 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 с кластерными устойчивыми стандартными ошибками и, похоже, напоминает результаты felm.

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

The 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)))