如何在R中按分组执行ANCOVA亚组分析并生成Anova结果？

错误原因 #

你的代码中do(model = lm(six_PHQ ~ Study.Arm + base_PHQ, data = PHQ1))直接调用了完整数据集PHQ1，而非group_by后的分组子集，导致两个婚姻状况组拟合的是同一个全局模型，所以系数和p值完全一致。

正确实现步骤 #

1. 加载所需包 #

先安装并加载必要的R包：

install.packages(c("dplyr", "car", "broom"))
library(dplyr)
library(car)
library(broom)

2. 分组拟合ANCOVA模型并提取结果 #

用.指代分组后的子集，结合Anova()生成Type III的ANCOVA结果：

Ancova_Marital <- PHQ1 %>%
  group_by(Marital_Status) %>%
  do(
    # 对每个亚组拟合ANCOVA模型
    model = lm(six_PHQ ~ Study.Arm + base_PHQ, data = .),
    # 提取Type III方差分析结果并整理为整洁数据框
    ancova_out = tidy(Anova(.$model, type = 3))
  ) %>%
  unnest(ancova_out) %>%
  # 筛选出干预组vs对照组的效应结果（可保留base_PHQ查看协变量效应）
  filter(term == "Study.Arm") %>%
  select(Marital_Status, term, estimate, p.value)

3. 查看最终结果 #

运行后打印输出：

print(Ancova_Marital)

此时会得到不同婚姻状况亚组中，调整基线PHQ评分后，干预组与对照组6个月抑郁评分的效应差异及对应的p值。

补充说明 #

• 若需要查看协变量base_PHQ的效应，删除filter(term == "Study.Arm")即可。
• Type III方差分析适用于不平衡数据集，若你的数据是平衡设计，Type I/II结果也可，但Type III的结果更稳健。
• broom包的作用是将统计模型输出转换为易读的数据框格式，方便后续整理。

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