lm_robust 模型的 ggpredict 误差包括固定效应和聚类标准误差

ggpredict error with lm_robust model including fixed effects and clustered standard errors

我是 运行 对两个自变量及其交互作用的回归,以及区域固定效应和使用来自 estimatrlm_robust 命令在观察级别对标准误差进行聚类(版本 0.22.0) 和 R 版本 3.6.0。

我想使用 ggpredict(版本 0.14.3)中的 plot 命令可视化预测的回归结果,但我收到一个错误,似乎是由于包含固定效应.

我得到的具体错误是: Error in X[, !beta_na, drop = FALSE] %*% coefs[!beta_na, ] : non-conformable arguments

如果我在 运行 后使用 ggpredict 回归,只聚集标准错误但不包括固定效应,代码运行得很好。使用来自 sjPlot 而不是 ggpredict.

的包装器命令时出现相同的错误

下面是一个MWE:

library(ggeffects)
library(estimatr)
library(sjPlot)
N <- 1000
df <- data.frame(id = rep(1:N),
                 district = as.factor(rep(1:20, times = 50)),
                 x = rpois(N, lambda = 4),
                 y = rnorm(N),
                 z = factor(rbinom(N, 1, prob = 0.5)))

mod1 <- lm_robust(y ~ x*z,
                  clusters = id,
                  fixed_effects = ~district,
                  data = df)
summary(mod1)
predDF <- ggpredict(mod1, terms = c("x", "z")) # use ggpredict from ggeffects
plot_model(mod1, type = "pred", terms = c("x", "z")) # using plot_model from sjPlot

想知道如何让 ggeffects/sjPlot 工作 lm_robust 包含聚类标准误差和固定效应的模型 - 或者替代包?我已经从 lfe 库中的 felm 切换到固定效果和聚类,因为 ggeffects 不适用于 felm 对象。

错误似乎与 estimatr 的预测方法有关。请参阅涉及 plm here.

的类似讨论

如果我们尝试使用 estimatr:::predict.lm_robust 进行预测,则会生成相同的错误:

library(ggeffects)
library(estimatr)
library(sjPlot)
N <- 1000
df <- data.frame(id = rep(1:N),
                 district = as.factor(rep(1:20, times = 50)),
                 x = rpois(N, lambda = 4),
                 y = rnorm(N),
                 z = factor(rbinom(N, 1, prob = 0.5)))

mod1 <- lm_robust(y ~ x*z,
                  clusters = id,
                  fixed_effects = ~district,
                  data = df)

new_df <- data.frame(x = rep(0:10, 2), z = factor(c(rep(0, 11), rep(1, 11))))

predict(mod1, newdata = new_df)
#> Error in X[, !beta_na, drop = FALSE] %*% coefs[!beta_na, ]: non-conformable arguments

如果不太麻烦,您可以通过在 lm() 中拟合固定效应模型来解决问题,在 ggeffects() 调用中对标准误差进行聚类。我知道用 lm() 估计固定效应并不理想,但它确实有效!

mod1 <- lm(y ~ x*z + factor(district), data = df)

predDF <- ggpredict(mod1, c("x", "z"), 
                    vcov.fun = "vcovCL", 
                    vcov.type = "HC1",
                    vcov.args = list(cluster = df$id))

plot(predDF)

reprex package (v0.3.0)

于 2020-04-27 创建