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

Question

I am running a regression on two independent variables and their interaction, along with district fixed effects and clustering the standard errors at the observation level using lm_robust command from estimatr (version 0.22.0) with R version 3.6.0.

I want to visualize the predicted regression results using the plot command from ggpredict (version 0.14.3) but am getting an error that appears to be due to the inclusion of the fixed effects.

The specific error I get is: Error in X[, !beta_na, drop = FALSE] %*% coefs[!beta_na, ] : non-conformable arguments

If I use ggpredict after running a regression that only clusters the standard errors but does not include fixed effects, the code runs just fine. I get the same error when using the wrapper commands from sjPlot instead of ggpredict.

Below is a 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

Wondering either how to get ggeffects/sjPlot to work lm_robust models containing both clustered standard errors and fixed effects - or alternative packages? I already switched from using felm in the lfe library for the fixed effects and clustering because ggeffects doesn't work on felm objects.

score 0 · Accepted Answer · answered Apr 28 '20 at 02:35

The error seems to be with estimatr's predict method. See a similar discussion involving plm here.

If we try to predict using estimatr:::predict.lm_robust, the same error is generated:

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

If it is not too onerous, you could do a workaround by fitting the fixed effects model in lm(), cluster the standard errors in the ggeffects() call. I know that estimating fixed effects with lm() is not ideal, but it works!

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)

^{Created on 2020-04-27 by the reprex package (v0.3.0)}

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

1 Answers1