R中如何获取转换后glm.nb各因子双变量均值及95%置信区间

模型适配性前置说明 #

负二项回归仅适用过离散计数型因变量场景：因变量为计数数据、且方差显著大于均值时，glm.nb选型合理；如果目标变量为连续型，请替换为普通线性模型lm()。示例代码中rnorm()生成的连续型x不符合负二项回归使用前提，请根据实际数据类型调整模型。

分组调整均值与95%置信区间计算 #

使用emmeans包计算控制协变量后的各组边际均值，输出结构完全匹配需要的df.mean格式：

# 加载依赖包
library(MASS)
library(ggplot2)
library(dplyr)
library(emmeans)

# 示例数据（实际使用替换为自有数据即可）
set.seed(123)
df <- data.frame(
  x = rpois(12, lambda = 4), # 示例用计数数据适配负二项模型
  y = rnorm(12, 6, 4), 
  tmt = rep(c("A","B","C"), each = 4)
)

# 拟合负二项回归模型
mod_x <- glm.nb(x ~ y + tmt, data = df)
# 计算x维度各组调整均值与95%CI（控制y为总体均值，结果返回原始尺度）
x_res <- as.data.frame(emmeans(mod_x, ~ tmt, type = "response")) %>%
  rename(
    x.mean.data = response,
    xaxis.CI.low = asymp.LCL,
    xaxis.CI.high = asymp.UCL
  ) %>%
  select(tmt, x.mean.data, xaxis.CI.low, xaxis.CI.high)

# 拟合y维度对应模型（y为连续型故使用线性模型）
mod_y <- lm(y ~ tmt, data = df)
# 计算y维度各组调整均值与95%CI
y_res <- as.data.frame(emmeans(mod_y, ~ tmt)) %>%
  rename(
    y.mean.data = emmean,
    yaxis.CI.low = lower.CL,
    yaxis.CI.up = upper.CL
  ) %>%
  select(tmt, y.mean.data, yaxis.CI.low, yaxis.CI.up)

# 合并为目标结果表df.mean，字段顺序与要求完全一致
df.mean <- left_join(x_res, y_res, by = "tmt") %>%
  select(x.mean.data, y.mean.data, yaxis.CI.low, yaxis.CI.up, xaxis.CI.low, xaxis.CI.high, tmt)

运行后输出的df.mean包含每个处理组的x/y调整均值、x/y维度的95%置信区间上下限、处理组标签，和给出的示例表结构完全匹配。

可视化实现 #

以下代码可直接生成和参考样式一致的带双变量置信区间的图形，支持添加分组联合置信椭圆：

ggplot(df.mean, aes(x = x.mean.data, y = y.mean.data, color = tmt)) +
  # 绘制y方向置信区间误差线
  geom_errorbar(aes(ymin = yaxis.CI.low, ymax = yaxis.CI.up), width = 0.1, linewidth = 1) +
  # 绘制x方向置信区间误差线
  geom_errorbarh(aes(xmin = xaxis.CI.low, xmax = xaxis.CI.high), height = 0.1, linewidth = 1) +
  # 绘制分组均值点
  geom_point(size = 3) +
  # 可选：添加和参考图一致的95%联合置信椭圆
  stat_ellipse(data = df, aes(x = x, y = y, color = tmt), linetype = 2, linewidth = 0.8) +
  labs(x = "X调整后均值", y = "Y调整后均值", color = "处理组") +
  theme_bw(base_size = 12)

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