lm_robust 模型的 ggpredict 错误，包括固定效应和聚集标准错误答案

【问题标题】：ggpredict error with lm_robust model including fixed effects and clustered standard errorslm_robust 模型的 ggpredict 错误，包括固定效应和聚集标准错误
【发布时间】：2020-08-11 15:03:10
【问题描述】：

我正在使用来自estimatr（版本 0.22.0）和 R 版本 3.6 的lm_robust 命令对两个自变量及其交互作用进行回归，以及区域固定效应和聚类观察级别的标准误差。 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 对象。

【问题讨论】：

标签： r sjplot

【解决方案1】：

错误似乎与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 年 4 月 27 日创建}

【讨论】：