首页 > 解决方案 > lm_robust 模型的 ggpredict 错误,包括固定效应和聚集标准错误

问题描述

我正在使用lm_robust来自estimatr(版本 0.22.0)和 R 版本 3.6.0 的命令对两个自变量及其交互作用以及区域固定效应和观测级别的标准误差进行回归分析。

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

我得到的具体错误是: 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包含集群标准错误和固定效果的模型 - 或替代包?我已经从felmlfe库中使用固定效果和聚类切换,因为ggeffects它不适用于felm对象。

标签: rsjplot

解决方案


错误似乎与estimatr's predict 方法有关。plm 请参阅此处涉及的类似讨论。

如果我们尝试使用 进行预测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 包(v0.3.0)于 2020-04-27 创建


推荐阅读