如何在sjPlot的plot_model中调整置信区间以反映聚类标准误

我用plot_model绘制两个变量的交互效应图，需要给置信区间加上聚类调整。已经用coeftest实现了聚类标准误的计算，还能导入stargazer表格，但coeftest对象和plot_model不兼容，导致画出的置信区间没做聚类调整。目前只能运行以下代码出图：

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"))

之前试过在plot_model里加vcov.fun参数，但置信区间完全没变化，怎么才能让图里的置信区间用上聚类调整？

方法1：给plot_model补全vcov参数 #

之前单独加vcov.fun没用，是因为没传递聚类变量的参数。plot_model支持直接调用vcovCL，需要同时指定vcov.fun和vcov.args：

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",
           # 添加聚类调整参数
           vcov.fun = "vcovCL",
           vcov.args = list(cluster = gis_dat$NomeDaCoalizao)) +
  scale_color_discrete(labels = c("0 - Lowest Coalescence", "1 - Highest Coalescence"))

这样plot_model会用vcovCL计算聚类调整后的标准误，进而生成正确的置信区间。

方法2：用marginaleffects包替代（更灵活） #

如果plot_model的参数设置还是有问题，可以换用marginaleffects包，它对聚类标准误的支持更直观，还能结合ggplot自由绘图：

library(marginaleffects)

# 计算带聚类调整的预测值
preds <- predictions(model4, 
                     # 指定要绘制的变量范围
                     variables = list(mwd = "grid", CoalescenceAmorim = c(0, 1)),
                     # 直接指定聚类变量
                     vcov = ~NomeDaCoalizao, 
                     data = gis_dat)

# 用ggplot绘图
ggplot(preds, aes(x = mwd, y = estimate, color = factor(CoalescenceAmorim))) +
  geom_line() +
  # 添加置信区间带
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = factor(CoalescenceAmorim)), 
              alpha = 0.2, color = NA) +
  labs(title = "Figure 6: Predicted Probability of Partial Vetoes on MPVs",
       x = "Coalition Heterogeneity",
       y = "Predicted Probability",
       color = "COALITION COALESCENCE",
       fill = "COALITION COALESCENCE") +
  scale_color_discrete(labels = c("0 - Lowest Coalescence", "1 - Highest Coalescence")) +
  scale_fill_discrete(labels = c("0 - Lowest Coalescence", "1 - Highest Coalescence"))

这个方法生成的预测值自带聚类调整的置信区间，绘图自由度更高。

内容的提问来源于stack exchange，提问作者Jair Moreira