如何在lmer/lm模型中以所有因子水平为对照组并生成对比矩阵？

在使用线性模型（lm）或线性混合效应模型（lmer）时，若模型包含无序分类变量，函数默认会将该变量的第一个因子水平作为对照组进行对比。我现在有一个多水平的分类变量，希望让这个变量的每个水平都依次作为对照基准组，并且能自动化完成这个过程，最终生成包含所有组间对比p值的矩阵。

以下是使用diamonds数据集的示例代码：

library(lmer)
library(lmerTest)    

# 创建无序分类变量
diamonds$color = factor(sample(c('red','white','blue','green','black'), nrow(diamonds), replace=T))
# 构建包含固定效应分类变量的lmer模型
mod = lmer(data=diamonds, carat~color+(1|clarity))
summary(mod, corr=F)

当前模型摘要显示，默认以black作为对照组，我需要将其他所有颜色水平也分别设为对照组：

采用REML拟合的线性混合模型，t检验使用Satterthwaite方法 [lmerModLmerTest]
公式: carat ~ color + (1 | clarity)
数据: diamonds

收敛时的REML准则值: 64684

标准化残差:
Min 1Q Median 3Q Max
-2.228 -0.740 -0.224 0.540 8.471

随机效应:
Groups Name Variance Std.Dev.
clarity (Intercept) 0.0763 0.276
Residual 0.1939 0.440
观测数: 53940，分组数: clarity, 8

固定效应:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.786709 0.097774 7.005805 8.05 8.7e-05 ***
colorblue -0.000479 0.005989 53927.996020 -0.08 0.94
colorgreen 0.007455 0.005998 53927.990722 1.24 0.21
colorred 0.000746 0.005986 53927.988909 0.12 0.90
colorwhite 0.000449 0.005971 53927.993708 0.08 0.94 #

显著性代码: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

方法1：使用emmeans包直接生成所有组间对比的p值矩阵 #

emmeans包可以直接对模型的边际均值进行两两对比，自动生成所有组间的p值，无需手动切换对照组，效率更高：

library(emmeans)

# 提取color变量的边际均值，并执行两两对比
emm_result = emmeans(mod, pairwise ~ color)
# 提取对比结果中的p值，整理为矩阵形式
p_matrix = matrix(
  summary(emm_result$contrasts)$p.value,
  nrow = length(levels(diamonds$color)),
  dimnames = list(levels(diamonds$color), levels(diamonds$color))
)
# 保留上三角（避免重复对比）
p_matrix[lower.tri(p_matrix)] = NA
print(p_matrix)

方法2：自定义循环切换对照组提取p值 #

如果需要手动控制每个水平作为对照组的过程，可以通过循环修改因子参考水平、重新拟合模型来提取p值：

# 获取color变量的所有水平
color_levels = levels(diamonds$color)
# 初始化p值矩阵
p_matrix = matrix(
  NA,
  nrow = length(color_levels),
  ncol = length(color_levels),
  dimnames = list(color_levels, color_levels)
)

# 循环将每个水平设为对照组
for (ref_level in color_levels) {
  # 修改因子的参考水平
  diamonds$color_ref = relevel(diamonds$color, ref = ref_level)
  # 重新拟合模型
  mod_ref = lmer(data = diamonds, carat~color_ref+(1|clarity))
  # 提取固定效应的p值（排除截距项）
  p_vals = summary(mod_ref, corr=F)$coefficients[-1, "Pr(>|t|)"]
  # 将p值填入矩阵对应位置
  p_matrix[ref_level, names(p_vals)] = p_vals
}

print(p_matrix)

注意：方法2需要多次拟合模型，在数据集较大时（如示例中的5万+观测）速度会明显慢于方法1，优先推荐使用emmeans方案。

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