R语言那些事 统计系列（十四） rstatix ANOVA分析 单因素 双因素 三因素

缘由

统计是生信的基石之一

碎碎念

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

Ref

ANOVA in R: The Ultimate Guide - Datanovia

单因素方差分析

library(tidyverse)
library(ggpubr)
library(rstatix)
 
data("PlantGrowth")
set.seed(1234)
PlantGrowth %>% sample_n_by(group, size = 1)
 
levels(PlantGrowth$group)
 
PlantGrowth <- PlantGrowth %>%
	#重排因子水平，确保 ctrl 为基线组，方便解释 ANOVA
	reorder_levels(group, order = c("ctrl", "trt1", "trt2"))
 
 
# boxplot 检查分布、极端值、是否对称等
PlantGrowth %>%
	group_by(group) %>%
	# 描述性统计（均值和标准差）
	get_summary_stats(weight, type = "mean_sd")
ggboxplot(PlantGrowth, x = "group", y = "weight")
 
#  检测组内异常值
PlantGrowth %>% 
	group_by(group) %>%
	identify_outliers(weight)
 
model  <- lm(weight ~ group, data = PlantGrowth)
# 残差分布 QQ 图检查正态性
ggqqplot(residuals(model))
# 全局 Shapiro 正态性检验（检验残差是否正态）
shapiro_test(residuals(model))
#每组单独做 Shapiro
PlantGrowth %>%
	group_by(group) %>%
	shapiro_test(weight)
ggqqplot(PlantGrowth, "weight", facet.by = "group")
 
#方差齐性检验（对 ANOVA 很关键）
plot(model, 1)
PlantGrowth %>% levene_test(weight ~ group)
 
# 单因素 ANOVA（方差齐时采用）
res.aov <- PlantGrowth %>% anova_test(weight ~ group)
res.aov
#ges 是广义效应大小（由于因子引起的变异量） 26 （26%） 的效应大小意味着 变化 weight 的 26% 可以解释治疗条件。
 
#Tukey 事后比较（方差齐）
pwc <- PlantGrowth %>% tukey_hsd(weight ~ group)
pwc
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)
	)
 
 
#如果方差不齐 → 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)
	)
 
# >>> pairwise_t_test(pool.sd=FALSE) 也是 Games-Howell 的等价替代
pwc3 <- PlantGrowth %>% 
	pairwise_t_test(
		weight ~ group, pool.sd = FALSE,
		p.adjust.method = "bonferroni"
	)
pwc3
 

双向 ANOVA

set.seed(123)
data("jobsatisfaction", package = "datarium")
jobsatisfaction %>% sample_n_by(gender, education_level, size = 1)
 
jobsatisfaction %>%
	group_by(gender, education_level) %>%
	get_summary_stats(score, type = "mean_sd")
 
# boxplot 按 gender 分面、education_level 着色
bxp <- ggboxplot(
	jobsatisfaction, x = "gender", y = "score",
	color = "education_level", palette = "jco"
)
bxp
 
# >>> 检查离群点
jobsatisfaction %>%
	group_by(gender, education_level) %>%
	identify_outliers(score)
 
model  <- lm(score ~ gender*education_level,
			 data = jobsatisfaction)
# >残差正态性
ggqqplot(residuals(model))
shapiro_test(residuals(model))
 
jobsatisfaction %>%
	group_by(gender, education_level) %>%
	shapiro_test(score)
 
ggqqplot(jobsatisfaction, "score", ggtheme = theme_bw()) +
	facet_grid(gender ~ education_level)
 
# 方差齐性
jobsatisfaction %>% levene_test(score ~ gender*education_level)
 
#分析
res.aov <- jobsatisfaction %>% anova_test(score ~ gender * education_level)
res.aov
 
# 如果交互显著，要分性别做简单主效应：
# 在男性内部比较 education_level；在女性内部比较 education_level
model <- lm(score ~ gender * education_level, data = jobsatisfaction)
jobsatisfaction %>%
	group_by(gender) %>%
	anova_test(score ~ education_level, error = model)
 
# pairwise emmeans，在每个性别组内做教育水平的事后比较 emmeans 用边际均值法，适合不平衡数据
pwc <- jobsatisfaction %>% 
	group_by(gender) %>%
	emmeans_test(score ~ education_level, p.adjust.method = "bonferroni") 
pwc
 
 
# 如果交互不显著 → 直接跑主效应的简单事后比较
jobsatisfaction %>%
	pairwise_t_test(
		score ~ education_level, 
		p.adjust.method = "bonferroni"
	)
 
# 也可使用 emmeans（不 group_by → 表示全局边际效应）
model <- lm(score ~ gender * education_level, data = jobsatisfaction)
jobsatisfaction %>% 
	emmeans_test(
		score ~ education_level, p.adjust.method = "bonferroni",
		model = model
	)
 
# >>> 两种 emmeans 的区别：
# >>> group_by(gender) = “简单主效应”
# >>> 不分组 = “整体主效应” （忽略 gender）
pwc <- pwc %>% add_xy_position(x = "gender")
bxp +
	stat_pvalue_manual(pwc) +
	labs(
		subtitle = get_test_label(res.aov, detailed = TRUE),
		caption = get_pwc_label(pwc)
	)

三因素方差分析

set.seed(123)
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
 
#异常值检测
headache %>%
	group_by(gender, risk, treatment) %>%
	identify_outliers(pain_score)
 
 
model  <- lm(pain_score ~ gender*risk*treatment, data = headache)
# 检查残差正态性
ggqqplot(residuals(model))
shapiro_test(residuals(model))
 
 
headache %>%
	group_by(gender, risk, treatment) %>%
	shapiro_test(pain_score)
 
ggqqplot(headache, "pain_score", ggtheme = theme_bw()) +
	facet_grid(gender + risk ~ treatment, labeller = "label_both")
#正态性
headache %>% levene_test(pain_score ~ gender*risk*treatment)
 
# 三因素 ANOVA：含 3 主效应 + 3 两两交互 + 1 三阶交互
res.aov <- headache %>% anova_test(pain_score ~ gender*risk*treatment)
res.aov
 
 
# ===========================================================
# 重点 Part：三阶交互显著 → 必须分析 “简单简单主效应”
# ===========================================================
 
model  <- lm(pain_score ~ gender*risk*treatment, data = headache)
 
headache %>%
	group_by(gender) %>%
	# 第一步：在性别内部做一个二因素 ANOVA（risk * treatment）
	# 检查交互：risk*treatment 是否依 gender 不同
	anova_test(pain_score ~ risk*treatment, error = model)
 
 
# 第二步：在 gender × risk 每个子组中比较 treatment
# 这是三因素 ANOVA 最关键、最容易搞混的步骤
treatment.effect <- headache %>%
	group_by(gender, risk) %>%
	#在男性-高风险、男性-低风险、女性-高风险、女性-低风险各自内部比较 treatment
	anova_test(pain_score ~ treatment, error = model)
treatment.effect %>% as_tibble() %>% dplyr::filter(gender == "male")
 
 
library(emmeans)
pwc <- headache %>%
	group_by(gender, risk) %>%
	#在每个 gender×risk 内做 treatment 的两两比较 这叫“简单简单比较（simple-simple comparisons）”
	emmeans_test(pain_score ~ treatment, p.adjust.method = "bonferroni") %>%
	select(-df, -statistic, -p) # Remove details
 
# 例如：取男性 + 高风险
pwc %>% filter(gender == "male", risk == "high")
 
# 把对应的边际均值提取出来（用于画图）
get_emmeans(pwc) %>% filter(gender == "male", risk == "high")
 
pwc <- pwc %>% add_xy_position(x = "treatment")
pwc.filtered <- pwc %>% filter(gender == "male", risk == "high")
 
bxp +
	stat_pvalue_manual(
		pwc.filtered, color = "risk", linetype = "risk", hide.ns = TRUE,
		tip.length = 0, step.increase = 0.1, step.group.by = "gender"
	) +
	labs(
		subtitle = get_test_label(res.aov, detailed = TRUE),
		caption = get_pwc_label(pwc)
	)

当三阶交互（ABC）显著时，你必须做：

第一层：对每个 A（gender）做 B C 的 ANOVA

group_by(gender) %>% anova_test(pain_score ~ risk * treatment)

---

第二层：在每个 A×B（gender×risk）子组内检查 treatment

group_by(gender, risk) %>% anova_test(pain_score ~ treatment)

---

第三层：做事后比较（pairwise）

group_by(gender, risk) %>% emmeans_test(pain_score ~ treatment)