在R中用stat_summary与lm模型绘制均值和CI的差异原因

问题描述 #

我需要为正态分布数据集绘制均值及95%置信区间，分别采用两种方式：一是使用ggplot2的stat_summary函数，二是拟合lm模型并提取系数绘图。但两种方式得到的图形存在差异，请问该差异的产生原因是什么？

方法一：使用ggplot2内置函数绘图 #

代码实现：

set.seed(123)
df <- data.frame(x.axis = rep(c("A","B","C","D"), each = 50),
y.axis = rnorm(1000, 1.2, 2))

ggplot(aes(x = x.axis, y = y.axis), data = df) + 
  geom_point(aes(x = x.axis, y = y.axis), color = "gray")+
  stat_summary(fun.data=mean_cl_boot, geom="errorbar", width=0.3, colour="black") + 
  stat_summary(fun.y = mean, color = "black", geom ="point", size = 5,show.legend = FALSE) + 
  theme_classic() +
  theme(aspect.ratio = 7/5)  + geom_hline(yintercept = 1) 

使用stat_summary绘制的图形：（图形显示各组均值及基于bootstrap的置信区间）

方法二：使用lm模型拟合绘图 #

代码实现：

mod <- lm(y.axis~x.axis, df)
summary(mod)
confint(mod)

summary_df <- summary(mod)$coefficients %>% 
  as_tibble() %>% 
  dplyr::mutate(x.axis = c("A","B","C","D")) %>%
  dplyr::rename(y.axis = Estimate) %>%
  mutate(ymin = y.axis - 1.96 * `Std. Error`,
         ymax = y.axis + 1.96 * `Std. Error`)

ggplot(aes(x = x.axis, y = y.axis,), data = df) + 
  geom_point(aes(x = x.axis, y = y.axis), color = "grey", data = df, inherit.aes=FALSE) + 
  geom_point(aes(x = x.axis, y = y.axis), color = 'black', data = summary_df,size = 5, inherit.aes=FALSE) +
  geom_errorbar(aes(x = x.axis, ymin = ymin, ymax = ymax,y = NULL), data = summary_df, inherit.aes=FALSE, width = 0.3) + 
  theme_classic() +
  theme(aspect.ratio = 7/5)  + geom_hline(yintercept = 1) 

使用lm模型绘制的图形：（图形显示各组均值及基于参数化方法的置信区间）

差异产生的核心原因 #

两种方法的差异本质是置信区间的计算逻辑完全不同：

• stat_summary的mean_cl_boot：该函数采用bootstrap自助法生成置信区间——对每组数据进行多次重复抽样，计算大量样本均值后取95%分位数区间，属于非参数估计，不依赖正态分布、方差齐性等假设，区间宽度由每组自身的数据分布决定。
• lm模型的置信区间：拟合的lm(y.axis~x.axis)本质是单因素方差分析模型，置信区间基于参数化方法计算：标准误由整体模型的合并残差方差推导而来，你手动使用的1.96是正态分布的95%临界值（正确做法应该用对应自由度的t分布临界值），这类区间依赖数据服从正态分布、各组方差齐性的假设。

此外还有两个细节放大了差异：

1. 你在lm部分误用了正态分布临界值1.96，实际小样本下应使用t分布临界值（不过本次每组50个样本，t值接近1.97，这不是差异的核心原因）。
2. lm模型的标准误基于合并方差（所有组残差方差的平均值）计算，而mean_cl_boot是每组单独计算方差，当各组方差存在差异时，两者的标准误会明显不同。

让两种方法结果一致的修正方案 #

如果要使两种方法的置信区间结果一致，可以选择以下两种方式：

1. 修改stat_summary的置信区间计算逻辑，改用参数化的mean_cl_normal：

ggplot(aes(x = x.axis, y = y.axis), data = df) + 
  geom_point(color = "gray")+
  stat_summary(fun.data=mean_cl_normal, geom="errorbar", width=0.3, colour="black") + 
  stat_summary(fun.y = mean, color = "black", geom ="point", size = 5,show.legend = FALSE) + 
  theme_classic() +
  theme(aspect.ratio = 7/5)  + geom_hline(yintercept = 1) 

2. 修改lm部分的置信区间提取方式，直接用模型自带的t分布置信区间：

library(broom)
summary_df <- tidy(mod, conf.int = TRUE) %>% 
  dplyr::mutate(x.axis = c("A","B","C","D")) %>%
  dplyr::rename(y.axis = estimate, ymin = conf.low, ymax = conf.high)

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