Ошибка ggpredict с моделью lm_robust, включая фиксированные эффекты и кластерные стандартные ошибки

Question

Я выполняю регрессию по двум независимым переменным и их взаимодействию, наряду с фиксированными эффектами округа и кластеризации стандартных ошибок на уровне наблюдения, используя команду lm_robust из estimatr (версия 0.22.0) с версией R 3.6. 0.

Я хочу визуализировать предсказанные результаты регрессии, используя команду plot из ggpredict (версия 0.14.3), но получаю ошибку, которая, по-видимому, связана с включением фиксированных эффектов.

Указанная ошибка c, которую я получаю: 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 модели, содержащие как кластерные стандартные ошибки, так и фиксированные эффекты - или альтернативные пакеты? Я уже перешел с использования felm в библиотеке lfe для фиксированных эффектов и кластеризации, потому что ggeffects не работает на felm объектах.

Answer 1

Ошибка, похоже, связана с методом прогнозирования estimatr. Смотрите подобное обсуждение с участием 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)

^{Создано в 2020-04-27 представьте пакет (v0.3.0)}