如何用R的lm函数估计多类别下组间直接均值差异？

需要用R中的lm函数估计多类别下处理组与对照组的均值差异（即A_treat - A_control、B_treat - B_control、C_treat - C_control）。目前用y~Categ:Group -1能得到各组单独均值，但用y~Categ*Group -1时，仅A类的差异正确，其余为相对差异，需要调整公式或找到更简便的方法。

示例数据与代码 #

## 创建数据，A组处理与对照组差异为0.5
Categ <- rep(LETTERS[1:3], times=400)
Group <- rep(c("Control", "Treat"), each = 600)

set.seed(123)
y <- as.numeric(as.factor(Categ))+
  rnorm(length(Group)) +
  ifelse(Categ=="A" & Group =="Treat",rnorm(mean=0.5, length(Group)),0)

df <- data.frame(Categ, Group, y)

## 估计各组单独均值
reg_full <- lm(y~Categ:Group -1, data=df)
coef(reg_full)
#> CategA:GroupControl CategB:GroupControl CategC:GroupControl   CategA:GroupTreat 
#>            1.002340            1.976993            3.085993            1.363267 
#>   CategB:GroupTreat   CategC:GroupTreat 
#>            2.039785            3.016188

## 尝试估计组间差异，但结果不符合预期
reg_diff <- lm(y~Categ*Group -1, data=df)
coef(reg_diff)
#>            CategA            CategB            CategC        GroupTreat 
#>         1.0023399         1.9769931         3.0859931         0.3609268 
#> CategB:GroupTreat CategC:GroupTreat 
#>        -0.2981345        -0.4307316

## 期望得到的直接差异：处理组减对照组的均值差
coef(reg_full)[4:6]-coef(reg_full)[1:3]
#> CategA:GroupTreat CategB:GroupTreat CategC:GroupTreat 
#>         0.3609268         0.0627923        -0.0698048
 
## 当前得到的结果是相对差异
coef(reg_diff)[4:6]
#>        GroupTreat CategB:GroupTreat CategC:GroupTreat 
#>         0.3609268        -0.2981345        -0.4307316

方法1：用嵌套公式直接估计组内差异 #

使用Categ/Group -1的公式，本质是拟合Categ + Categ:Group -1，其中Categ:Group的系数就是每个类别下处理组与对照组的直接差异：

model <- lm(y ~ Categ/Group - 1, data = df)
coef(model)

输出中CategA:GroupTreat、CategB:GroupTreat、CategC:GroupTreat的系数就是对应类别Treat - Control的均值差，和你期望的结果完全一致。

方法2：用分组对比工具直接提取差异 #

借助emmeans包可以更直观地提取每个类别下的组间差异，同时给出统计显著性：

# 先拟合标准交互模型（不需要去掉截距）
df$Group <- relevel(factor(df$Group), ref = "Control")
model <- lm(y ~ Categ * Group, data = df)

# 提取每个类别下处理组与对照组的对比结果
library(emmeans)
emmeans(model, pairwise ~ Group | Categ)$contrasts

方法3：手动构造对比矩阵计算差异 #

如果不想用额外包，可以基于全均值模型构造对比矩阵，直接计算目标差异：

reg_full <- lm(y ~ Categ:Group - 1, data = df)

# 构造对比矩阵：每行对应一个类别Treat - Control的对比
contrasts <- matrix(
  c(-1, 0, 0, 1, 0, 0,   # A_Treat - A_Control
    0, -1, 0, 0, 1, 0,   # B_Treat - B_Control
    0, 0, -1, 0, 0, 1),  # C_Treat - C_Control
  nrow = 3, byrow = TRUE,
  dimnames = list(c("A_Treat-Control", "B_Treat-Control", "C_Treat-Control"),
                  names(coef(reg_full)))
)

# 计算差异及显著性
library(multcomp)
summary(glht(reg_full, linfct = contrasts))

方法4：分组计算均值差（最简方式） #

用dplyr直接按类别分组计算均值差，适合快速查看结果：

library(dplyr)
df %>%
  group_by(Categ, Group) %>%
  summarise(mean_y = mean(y), .groups = "drop") %>%
  tidyr::pivot_wider(names_from = Group, values_from = mean_y) %>%
  mutate(treat_control_diff = Treat - Control)

如果需要统计检验，可嵌套t.test：

df %>%
  group_by(Categ) %>%
  summarise(
    diff = t.test(y ~ Group)$estimate[2] - t.test(y ~ Group)$estimate[1],
    p_value = t.test(y ~ Group)$p.value,
    .groups = "drop"
  )

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