缘由

统计是生信的基石之一

碎碎念

要爱一个具体的人，不要爱一个抽象的人；要爱生活，不要爱生活的意义。 ——《卡拉马佐夫兄弟》

Ref

ANOVA in R: The Ultimate Guide - Datanovia

解释

一、单因素方差分析（One-way ANOVA）

我们先使用内置数据集 PlantGrowth 来进行单因素方差分析，逐步掌握统计流程。

library(tidyverse)
library(ggpubr)
library(rstatix)
 
data("PlantGrowth")
set.seed(1234)
PlantGrowth %>% sample_n_by(group, size = 1)
levels(PlantGrowth$group)
 
PlantGrowth %>%
  reorder_levels(group, order = c("ctrl", "trt1", "trt2"))

（1）描述性统计

PlantGrowth %>%
  group_by(group) %>%
  get_summary_stats(weight, type = "mean_sd")

（2）初步可视化

ggboxplot(PlantGrowth, x = "group", y = "weight")

通过箱线图我们可以快速观察三组植物的生长重量分布情况。方差分析的核心思想就是——比较三组均值是否存在显著差异。

（3）异常值检测

PlantGrowth %>% 
  group_by(group) %>%
  identify_outliers(weight)

如果存在明显异常值，建议使用 WRS2 包进行稳健方差分析（robust ANOVA），以防结果偏差。

（4）正态性检验

整体残差检验：

model  <- lm(weight ~ group, data = PlantGrowth)
ggqqplot(residuals(model))
shapiro_test(residuals(model))

分组正态性检验：

PlantGrowth %>%
  group_by(group) %>%
  shapiro_test(weight)
ggqqplot(PlantGrowth, "weight", facet.by = "group")

（5）方差齐性检验

plot(model, 1)
PlantGrowth %>% levene_test(weight ~ group)

（6）方差齐时的单因素方差分析

res.aov <- PlantGrowth %>% anova_test(weight ~ group)
res.aov

输出中的 ges 表示广义 Eta 平方（效应量），反映了因变量变异中可由分组解释的比例。例如 ges = 0.26 表示 26% 的植物重量变异可由处理方式解释。

（7）Tukey 事后比较与可视化

pwc <- PlantGrowth %>% tukey_hsd(weight ~ group)
pwc <- pwc %>% add_xy_position(x = "group")
 
ggboxplot(PlantGrowth, x = "group", y = "weight") +
  stat_pvalue_manual(pwc, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )

（8）方差不齐时：Welch ANOVA + Games-Howell 检验

res.aov2 <- PlantGrowth %>% welch_anova_test(weight ~ group)
pwc2 <- PlantGrowth %>% games_howell_test(weight ~ group)
pwc2 <- pwc2 %>% add_xy_position(x = "group", step.increase = 1)
 
ggboxplot(PlantGrowth, x = "group", y = "weight") +
  stat_pvalue_manual(pwc2, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov2, detailed = TRUE),
    caption = get_pwc_label(pwc2)
  )

二、双因素方差分析（Two-way ANOVA）

数据集：jobsatisfaction
因变量：score
自变量：gender（性别）与 education_level（教育水平）

set.seed(123)
data("jobsatisfaction", package = "datarium")
jobsatisfaction %>% sample_n_by(gender, education_level, size = 1)

（1）描述性与可视化

jobsatisfaction %>%
  group_by(gender, education_level) %>%
  get_summary_stats(score, type = "mean_sd")
 
bxp <- ggboxplot(
  jobsatisfaction, x = "gender", y = "score",
  color = "education_level", palette = "jco"
)
bxp

（2）模型与正态性检验

model  <- lm(score ~ gender*education_level, data = jobsatisfaction)
ggqqplot(residuals(model))
shapiro_test(residuals(model))
jobsatisfaction %>% levene_test(score ~ gender*education_level)

（3）双因素 ANOVA

res.aov <- jobsatisfaction %>% anova_test(score ~ gender * education_level)
res.aov

如果交互项 gender:education_level 显著（p < 0.05），说明教育水平的影响依赖于性别。

（4）为什么需要指定 `error = model`？

jobsatisfaction %>%
  group_by(gender) %>%
  anova_test(score ~ education_level, error = model)

在多因素设计中，model 决定了误差来源。
不加 error = model，R 会在组内重新拟合简单模型，忽略其他主效应与交互作用，导致方差划分错误。

三、为什么使用 `emmeans` 而不是 `pairwise_t_test`

pairwise_t_test()：
直接进行两两比较，不考虑整体模型结构，适合单因素或独立样本设计。

emmeans_test()：
基于线性模型的估计均值（Estimated Marginal Means），可处理：

不平衡样本；
多因素交互；
协变量；
精确控制误差来源。

library(emmeans)
pwc <- jobsatisfaction %>% 
  group_by(gender) %>%
  emmeans_test(score ~ education_level, p.adjust.method = "bonferroni")
pwc

四、主效应与交互效应的关系

效应类型	含义	举例解释
主效应	一个因素对因变量的平均影响	性别总体上是否影响工作满意度
交互效应	一个因素的影响依赖于另一个因素	性别对满意度的影响是否因教育水平不同而不同

解释规则：

1️交互效应显著 → 进行简单效应分析
2️交互效应不显著 → 直接解释主效应
3️主效应显著 → 做事后比较
4️ 主效应不显著 → 不做进一步分析

五、简单效应分析详解

5.1 什么是简单效应分析？

简单效应分析是指在多因素 ANOVA 中，当交互作用显著时，在另一因素的特定水平上分别考察一个因素的效应。

5.2 理解示例：`jobsatisfaction` 数据

jobsatisfaction %>%
  group_by(gender, education_level) %>%
  get_summary_stats(score, type = "mean_sd")

假设发现：

男性不同教育水平间满意度差异大；
女性不同教育水平间差异小。

→ 教育水平的效应依赖于性别 → 存在交互作用。

5.3 实际操作步骤

Step 1：检查交互显著性

res.aov <- jobsatisfaction %>% anova_test(score ~ gender * education_level)
res.aov

Step 2：在性别的不同水平上分析教育水平

model <- lm(score ~ gender * education_level, data = jobsatisfaction)
simple_effects <- jobsatisfaction %>%
  group_by(gender) %>%
  anova_test(score ~ education_level, error = model)
simple_effects

Step 3：结果解释示例

gender  Effect           p.value
female  education_level  0.850
male    education_level  0.003

解释：

女性组中教育水平无显著影响；
男性组中教育水平显著影响满意度。

5.4 进一步的事后比较（以男性为例）

male_comparisons <- jobsatisfaction %>%
  filter(gender == "male") %>%
  emmeans_test(score ~ education_level, p.adjust.method = "bonferroni")
male_comparisons

5.5 可视化

ggline(jobsatisfaction, x = "education_level", y = "score", 
       color = "gender", palette = "jco",
       add = c("mean_se", "dotplot"),
       facet.by = "gender") +
  labs(title = "简单效应分析可视化",
       subtitle = "教育水平对满意度的影响在男女中不同")

5.6 主效应 vs 简单效应

比喻：

研究“学习方法(A、B)”和“教学形式(线上、线下)”对考试成绩的影响。

分析	结果	含义
主效应	学习方法不显著	平均来看 A 与 B 没差异
交互效应	显著	教学形式影响学习方法的效果
简单效应	线下A>B（显著）线上A=B（不显著）	学习方法效果取决于教学形式

5.7 交互作用不显著时的正确做法

如果交互不显著：

不做简单效应分析；
直接解释主效应。

res.aov <- jobsatisfaction %>% anova_test(score ~ gender * education_level)
if(res.aov$p[3] > 0.05) {
  res.main <- jobsatisfaction %>% anova_test(score ~ gender + education_level)
  res.main
}

再做主效应的事后比较：

education_main <- jobsatisfaction %>% 
  pairwise_t_test(score ~ education_level, p.adjust.method = "bonferroni")
 
gender_main <- jobsatisfaction %>% 
  pairwise_t_test(score ~ gender, p.adjust.method = "bonferroni")

5.8 完整决策逻辑

判断条件	下一步分析	说明
交互显著	简单效应分析	分层比较，探讨“何时有差异”
交互不显著	主效应分析	直接报告总体影响
主效应显著	事后比较	找出具体差异
主效应不显著	无需后续分析	结束分析

六、三因素方差分析（Three-way ANOVA）

数据集：headache
自变量：gender、risk、treatment
因变量：pain_score

data("headache", package = "datarium")
headache %>% sample_n_by(gender, risk, treatment, size = 1)

描述与可视化

headache %>%
  group_by(gender, risk, treatment) %>%
  get_summary_stats(pain_score, type = "mean_sd")
 
bxp <- ggboxplot(headache, x = "treatment", y = "pain_score", 
                 color = "risk", palette = "jco", facet.by = "gender")
bxp

三因素 ANOVA

res.aov <- headache %>% anova_test(pain_score ~ gender*risk*treatment)
res.aov

若性别 × 风险 × 治疗显著，则说明治疗效果受性别和风险的双重影响。

分层分析与事后比较

model <- lm(pain_score ~ gender*risk*treatment, data = headache)
 
headache %>%
  group_by(gender) %>%
  anova_test(pain_score ~ risk*treatment, error = model)
 
pwc <- headache %>%
  group_by(gender, risk) %>%
  emmeans_test(pain_score ~ treatment, p.adjust.method = "bonferroni")
# Show comparison results for male at high risk
pwc %>% filter(gender == "male", risk == "high")
get_emmeans(pwc) %>% filter(gender == "male", risk == "high")

七、总结与核心理解

model 参数决定误差划分
必须指定，否则多因素设计的方差解释会错误。

交互显著 → 简单效应分析；不显著 → 主效应解释
切记不要在交互不显著时强行做简单效应。

多因素实验分析逻辑：

主效应显著？——> 是 → 做事后比较
交互显著？——> 是 → 做简单效应
           否 → 解释主效应

lvzhaopeng的数字花园

探索

R语言那些事 统计系列（九） rstatix 从单因素到三因素方差分析 含简单效应分析与模型解释

缘由

碎碎念

Ref

解释