如何在R中为多项逻辑回归模型的效应图添加置信区间？

我需要为多项逻辑回归（multinomial logistic regression）模型的效应图添加置信区间（confidence intervals，或误差线）。以下是可复现的数据、模型及绘图代码：

install.packages("effects")
library("effects")
ToyData2 <- data.table(Answer = c("Remained incorrect", "Changed to correct", "Changed to correct", "Remained incorrect",
                             "Remained correct", "Changed to incorrect", "Remained correct", "Changed to correct",
                             "Changed to correct", "Changed to incorrect", "Remained incorrect", "Remained incorrect",
                             "Remained correct", "Changed to incorrect", "Remained correct", "Changed to correct"),
                   Sex = c("Male", "Male", "Female", "Male", "Female", "Female", "Male", "Male", 
                           "Female", "Male", "Female", "Female", NA, "Male", "Female", "Male"),
                   Education = c("Less than Bachelors", "Less than Bachelors", "Bachelors",NA,
                                 "Post Bachelors","Post Bachelors","Bachelors", "Less than Bachelors",
                                 "Bachelors",NA,"Post Bachelors","Post Bachelors",
                                 "Post Bachelors", "Less than Bachelors", "Bachelors","Post Bachelors"),
                   Minutes = c(1, 5, 1, 29, NA, 2, 1, 3, 4, 5, 1, 2, 2, 1, 3, 4),
                   Problem = c("A","A","A","A", "B","B","B","B", "C","C","C","C", "D","D","D","D"))
ToyData2

ToyData2$Answer <- factor(ToyData2$Answer, ordered = FALSE)
summary(ToyData2$Answer)
ToyData2$Answer <- relevel(ToyData2$Answer, ref = 'Changed to incorrect')

ToyData2$Education <- factor(ToyData2$Education, levels = c("Less than Bachelors","Bachelors","Post Bachelors"))

ToyOutcomeLogit <- multinom(Answer ~ Sex + Education + Minutes + Problem, data = ToyData2)
summary(ToyOutcomeLogit)
table(ToyData2$Minutes)
plot(Effect('Minutes', ToyOutcomeLogit), multiline=T)
table(ToyData2$Education)
plot(Effect('Sex', ToyOutcomeLogit), multiline=T)

上述代码生成两个效应图：

1. 时间（Minutes）的效应图
2. 教育程度（Education）的效应图

我当前使用的R包如下：

library('pacman')
p_load(
  'data.table',
  'DescTools',
  'effects',
  'ggpubr',
  'ggsignif',
  'glue',
  'Hmisc',
  'irr',
  'lm.beta',
  'nnet',
  'openxlsx',
  'psych',
  'scales',
  'sjPlot',
  'stats',
  'tidyr',
  'tidyverse'
)

具体问题：
A. 如何调整现有代码以添加置信区间？
B. 还有哪些函数可用于绘制带置信区间的效应图？
C. 应使用哪些其他包的函数？

（注：当前R版本为4.2.2）

A. 调整现有代码添加置信区间 #

你使用的effects包本身支持在效应图中添加置信区间，只需在Effect()中指定置信水平，再在plot()里开启ci=TRUE即可：

基础调整示例 #

# 绘制Minutes的效应图并添加95%置信区间
plot(Effect('Minutes', ToyOutcomeLogit, confidence.level = 0.95), multiline=T, ci=TRUE)

# 绘制Sex的效应图并添加95%置信区间
plot(Effect('Sex', ToyOutcomeLogit, confidence.level = 0.95), multiline=T, ci=TRUE)

自定义绘图（ggplot2） #

如果需要更灵活的样式控制，可以先提取效应数据，再用ggplot2手动绘制：

# 提取Education的效应数据（包含置信区间上下限）
eff_edu <- Effect('Education', ToyOutcomeLogit, confidence.level = 0.95)
eff_edu_df <- as.data.frame(eff_edu)

# ggplot2绘制带置信区间的分组效应图
ggplot(eff_edu_df, aes(x=Education, y=fitted, color=Answer)) +
  geom_point(position=position_dodge(width=0.5)) +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=0.2, position=position_dodge(width=0.5)) +
  labs(title="教育程度对回答结果的效应（95%置信区间）", x="教育程度", y="拟合概率") +
  theme_minimal()

B. 其他可绘制带置信区间效应图的函数 #

除了effects包的原生绘图函数，还有这些实用选项：

• sjPlot包的plot_model()：专为统计模型可视化设计，支持多项逻辑回归，一键生成带置信区间的效应图
• ggeffects包的ggpredict()：生成模型预测值及置信区间，直接输出ggplot风格图形
• emmeans结合ggplot2：先计算边际均值与置信区间，再手动绘图，灵活性拉满

C. 推荐使用的其他包及函数 #

1. sjPlot包（快速可视化） #

plot_model()代码简洁，支持直接指定分组和置信水平：

# 绘制Minutes的连续变量效应图
plot_model(ToyOutcomeLogit, type = "eff", terms = "Minutes", ci.lvl = 0.95)

# 绘制Sex的分类变量效应图，按Answer分组
plot_model(ToyOutcomeLogit, type = "eff", terms = c("Sex", "Answer"), ci.lvl = 0.95)

2. ggeffects包（ggplot风格） #

ggpredict()返回ggplot对象，方便后续自定义修改：

library(ggeffects)

# 获取Minutes的预测值及置信区间
ggp_minutes <- ggpredict(ToyOutcomeLogit, terms = "Minutes", ci.lvl = 0.95)
plot(ggp_minutes) + theme_minimal() + labs(x="用时（分钟）", y="拟合概率")

# 获取Education的分组预测值
ggp_edu <- ggpredict(ToyOutcomeLogit, terms = c("Education", "Answer"), ci.lvl = 0.95)
plot(ggp_edu) + theme_minimal()

3. emmeans + ggplot2（高度自定义） #

先通过emmeans计算边际均值，再用ggplot2自由调整样式：

library(emmeans)

# 计算Education分组下的边际概率及置信区间
emm_edu <- emmeans(ToyOutcomeLogit, ~ Answer | Education, type = "response")
emm_edu_df <- as.data.frame(emm_edu)

# 绘制带置信区间的点线图
ggplot(emm_edu_df, aes(x=Education, y=prob, color=Answer, group=Answer)) +
  geom_line(position=position_dodge(0.5)) +
  geom_point(position=position_dodge(0.5)) +
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=0.2, position=position_dodge(0.5)) +
  theme_minimal()

内容的提问来源于stack exchange，提问作者Nick Byrd