Regression confidence bands from plm models

Question

I am currently running a fixed effects model with the use of the plm function from the package of the same name. I have clustering in the data, so I am using the vcovCR function from the clubSandwich package to get the standard errors of the estimated regression coefficients.

What I want to do is the following: Create a graph that shows the predicted value of Y given X1 (when the other variables are held constant at their mean values) and also add 95% confidence bands to the plot.

What I have tried:

1. Using the ggpredict() function from the ggeffects package. This did not work due to the plm package apparently being problematic for ggeffects.
2. Manually using the "group centered" approach of fixed effects estimation, and then running the usual lm() function on the group centered data. Running ggpredict() on that data does give me nice confidence bands, but the values of X and Y are then group centered (surprise, surprise!) and I want the predictions and values of X to be on the original scale, not on the mean centered scale.

Any suggestions on how to produce confidence bands for fixed effect estimation, when using cluster robust standard errors would be much appreciated.

Answer 1

I was trying something similar, however, I could not find a straight forward solution. I am estimating an event study with lags and leads on a fe model. The following worked for me:

ES1<-plm(total_thefts ~ lag6 + lag5 + lag4 + lag3 + lag2 + E + lead1 + lead2 + lead3 + lead4 + lead5 + lead6 + factor(month) + factor(year), data = panel, model="within")

Since, plm does not have a straight forward way to estimate robust standard errors, I do that manually and replace the original ones:

G <- length(unique(panel$id))
c <- G/(G - 1)
robust_se <- coeftest(ES1, c * vcovHC(ES1, type = 'HC1', cluster = 'group'))
ES1<-summary(ES1)
ES1$coefficients[,2:4] <- robust_se[,2:4]

After that, I separate my coefficients of interest and estimate the CI bands (95%), following this post:

coefs <- tidy(ES1[["coefficients"]])
coefs$conf.low <- coefs$Estimate+c(-1)*coefs$Std..Error*qt(0.975,42)
coefs$conf.high <- coefs$Estimate+c(1)*coefs$Std..Error*qt(0.975,42)

Finally, after that I keep my parameters of interest and plot the event study:

interest <- c("lag6","lag5","lag4","lag3","lag2","E","lead1","lead2","lead3","lead4","lead5","lead6")
coefs <- subset(coefs,coefs$.rownames%in%interest)
coefs$time <- c(seq(-6,-2,1),seq(0,6,1))

ggplot(coefs, aes(time, Estimate))+
           geom_line() +
           geom_point()+
           geom_pointrange(aes(ymin = conf.low, ymax = conf.high))+
           geom_hline(yintercept = 0, linetype = 2) +
           ggtitle(paste0("N: ",nobs(ES1),", Sample: ",pdim(ES1)$nT$n,", Buff: 1000m"))) +
           theme(plot.title = element_text(hjust = 0.5)))