如何修正R中emmeans包计算估计边际均值(EMMs)的标准误异常问题

你观测到的标准误完全一致是整体线性模型的方差齐性默认假设导致的：
你当前拟合的lm默认假设所有观测组的残差方差完全相同，计算标准误时调用的是全数据集统一估算的残差均方（MSE），而你每个Age*Location*Treatment分组的样本量完全相等（都是3），因此两组差值的标准误计算公式为SE = √(MSE*(1/n1 + 1/n2))，所有分组算出来的结果自然完全相同。

如果你希望每个分组的标准误仅基于该分组的观测方差计算，和单独做t检验的逻辑一致，可以参考以下两种修正方案：

方案1：按分组独立拟合模型（完全匹配独立t检验结果） #

用lme4的lmList对每个Age*Location子集单独拟合Value ~ Treatment模型，再用emmeans计算差值：

library(lme4)
# 按Age和Location分组，每组单独拟合模型
m_list <- lmList(Value ~ Treatment | Age + Location, data = df1)
# 计算组间差值
emm_list <- emmeans(m_list, ~ Treatment | Age + Location)
emm_pairs_list <- pairs(emm_list)
# 此时的SE和每组单独做t检验的结果完全一致
emm_pairs_list
# 绘图
plot(emm_pairs_list, by = "Location", CIs = TRUE) +
  labs(x = "均值差值 (± SE)",
       y = "年龄") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_y_discrete(labels = c("1","2","3")) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        text = element_text(size = 16),
        axis.text.x = element_text(size = 14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"),
        strip.background =element_rect(fill="white"))

方案2：在整体模型中纳入异方差结构 #

如果你需要保留整体模型的框架，同时允许不同组方差不同，可以用nlme的gls拟合带异方差的模型：

library(nlme)
# 拟合允许不同Age*Treatment*Location组方差不同的模型
m_gls <- gls(Value ~ Age * Treatment * Location, 
             weights = varIdent(form = ~1 | Age*Treatment*Location),
             data = df1)
# 计算emmeans和差值
emm_gls <- emmeans(m_gls, ~ Treatment | Age, by = "Location")
emm_pairs_gls <- pairs(emm_gls, by = c("Age","Location"))
emm_pairs_gls

两种方案输出的差值标准误都会随各组的实际方差变化，不再是统一值。

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