Search code examples
rsjplot

Adjust Confidence Intervals in SJ Plot to reflect clustered standard errors

I am trying to use the fuction plot_modelto to graph an interaction between two variables. I need to include clusters for the variables, I have done so using coeftest and it works well and can even get it into an stargazer table. However the coeftest object does not seem to work too well with plot_model. The result is suuch that the confidence intervals are not adjusted for the cluster as I have only been able to plot the following:

model4 <- glm(TotalOrPartialVetoDummy ~  mwd*CoalescenceAmorim+CoalitionPercentageAmorim+Disciplinados +HowLongIntoPresidency + as.factor(honeymoonDummy100days) + as.factor(VotNominal),
             data = gis_dat,
             family = binomial,
             weights = ifelse(TIPOLEI == 1, 1, 0))
cluster_model4 <- coeftest(model4, vcovCL(model4, cluster = gis_dat$NomeDaCoalizao))

plot_model(model4, type = "pred", terms = c("mwd", "CoalescenceAmorim [0,1]"), title = "Figure 6: Predicted Probability of Partial Vetoes on MPVs", axis.title = c("Coalition Heterogeneity", "Predicted Probability"), legend.title = "COALITION COALESCENCE")+
  scale_color_discrete(labels = c("0 - Lowest Coalescence", "1 - Highest Coalescence"))

Does anyone know how I might be able to add the clusters.

I have tried adding vcov.fun to the plot_modelbut it does not seem to produce any changes whatsoever on the confidence intervals and the information within the plot.

Solution

  • This is not a plot_model() solution, but perhaps some readers will be interested in this approach using the marginaleffects package.

    Here is an example with predicted probabilities along with confidence intervals built using a clustered standard errors by firm:

    library(marginaleffects)
dat <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/plm/EmplUK.csv")
dat <- transform(dat, BigW = emp > median(emp), BigK = wage > median(wage))
mod <- glm(BigW ~ output * BigK, data = dat, family = binomial)

plot_predictions(mod, 
    condition = c("output", "BigK"),
    vcov = ~ firm)