缘由

统计是生信的基石之一

碎碎念

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

Ref

Repeated Measures ANOVA in R: The Ultimate Guide - Datanovia

解读

一、准备阶段

1. 加载所需包

library(tidyverse)  # 数据整理
library(ggpubr)     # 图形展示（对 ggplot2 的高级封装）
library(rstatix)    # 统计分析核心包

二、单因素重复测量 ANOVA

1. 数据来源

我们使用 datarium 包中的 selfesteem 数据集，研究 12 个被试在 3 个时间点（t1, t2, t3）的自尊分数变化。

data("selfesteem", package = "datarium")
head(selfesteem, 3)

输出（宽格式）：

# A tibble: 3 × 4
     id    t1    t2    t3
  <int> <dbl> <dbl> <dbl>
1     1  4.01  5.18  7.11
2     2  2.56  6.91  6.31
3     3  3.24  4.44  9.78

2. 转换成长格式

重复测量 ANOVA 需要“长格式”（每一行是一次观测）：

selfesteem <- selfesteem %>%
  gather(key = "time", value = "score", t1, t2, t3) %>%
  convert_as_factor(id, time)
 
head(selfesteem, 3)

3. 描述统计

selfesteem %>%
  group_by(time) %>%
  get_summary_stats(score, type = "mean_sd")

4. 可视化

bxp <- ggboxplot(selfesteem, x="time", y="score", add="point")
bxp

可以直观看出不同时间点的分数变化。

5. 检查假设

(1) 异常值

selfesteem %>% group_by(time) %>% identify_outliers(score)

(2) 正态性

selfesteem %>% group_by(time) %>% shapiro_test(score)
ggqqplot(selfesteem, "score", facet.by="time")

6. 方差分析

res.aov <- anova_test(data = selfesteem, dv = score, wid = id, within = time)
get_anova_table(res.aov)

输出示例：

Effect DFn DFd      F        p p<.05   ges
time     2  18  55.5 2.01e-08     *  0.83

说明时间点对自尊分数有显著影响。

7. 事后检验（配对 t 检验）

pwc <- selfesteem %>%
  pairwise_t_test(score ~ time, paired = TRUE, p.adjust.method = "bonferroni")
pwc <- pwc %>% add_xy_position(x="time")
 
bxp + 
  stat_pvalue_manual(pwc) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
  )

三、双因素重复测量 ANOVA

现在我们加入第二个重复测量因素：处理（treatment）。

1. 数据准备

data("selfesteem2", package = "datarium")
 
selfesteem2 <- selfesteem2 %>%
  gather(key = "time", value = "score", t1, t2, t3) %>%
  convert_as_factor(id, time)
 
head(selfesteem2, 3)

2. 描述统计

selfesteem2 %>%
  group_by(treatment, time) %>%
  get_summary_stats(score, type = "mean_sd")

3. 可视化

bxp2 <- ggboxplot(selfesteem2, x="time", y="score",
                  color="treatment", palette="jco")
bxp2

4. 检查假设

selfesteem2 %>% group_by(treatment, time) %>% identify_outliers(score)
selfesteem2 %>% group_by(treatment, time) %>% shapiro_test(score)
ggqqplot(selfesteem2, "score", ggtheme = theme_bw()) +
  facet_grid(time ~ treatment, labeller = "label_both")

5. 双因素重复测量 ANOVA

res.aov2 <- anova_test(
  data = selfesteem2, dv = score, wid = id,
  within = c(treatment, time)
)
get_anova_table(res.aov2)

输出示例：

ANOVA Table (type III tests)

          Effect  DFn   DFd      F        p p<.05   ges
1      treatment 1.00 11.00 15.541 2.00e-03     * 0.059
2           time 1.31 14.37 27.369 5.03e-05     * 0.049
3 treatment:time 2.00 22.00 30.424 4.63e-07     * 0.050

三个效应都显著，说明不同处理组随时间变化趋势不同。

6. 事后检验：针对显著交互作用

因为 treatment:time 显著，我们需要看 不同 treatment 水平下时间的变化。

pwc2 <- selfesteem2 %>%
  group_by(treatment) %>%
  pairwise_t_test(score ~ time, paired = TRUE, p.adjust.method = "bonferroni")
 
pwc2 <- pwc2 %>% add_xy_position(x = "time")
 
bxp2 + 
  stat_pvalue_manual(pwc2, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov2, detailed = TRUE),
    caption = get_pwc_label(pwc2)
  )

四、三因素重复测量 ANOVA

使用 weightloss 数据集。
包含三种因素：

饮食（diet：yes/no）
运动（exercises：yes/no）
时间（time：t1, t2, t3）

1. 数据准备

data("weightloss", package = "datarium")
 
weightloss <- weightloss %>%
  gather(key = "time", value = "score", t1, t2, t3) %>%
  convert_as_factor(id, time)
 
head(weightloss, 3)

2. 描述统计

weightloss %>%
  group_by(diet, exercises, time) %>%
  get_summary_stats(score, type = "mean_sd")

3. 可视化

bxp3 <- ggboxplot(
  weightloss, x = "exercises", y = "score",
  color = "time", palette = "jco",
  facet.by = "diet", short.panel.labs = FALSE
)
bxp3

4. 检查假设

weightloss %>%
  group_by(diet, exercises, time) %>%
  identify_outliers(score)
 
weightloss %>%
  group_by(diet, exercises, time) %>%
  shapiro_test(score)
 
ggqqplot(weightloss, "score", ggtheme = theme_bw()) +
  facet_grid(diet + exercises ~ time, labeller = "label_both")

5. 三因素重复测量 ANOVA

res.aov3 <- anova_test(
  data = weightloss, dv = score, wid = id,
  within = c(diet, exercises, time)
)
get_anova_table(res.aov3)

示例输出：

              Effect  DFn   DFd       F        p p<.05   ges
1                diet 1.00 11.00   6.021 3.20e-02     * 0.028
2           exercises 1.00 11.00  58.928 9.65e-06     * 0.284
3                time 2.00 22.00 110.942 3.22e-12     * 0.541
4      diet:exercises 1.00 11.00  75.356 2.98e-06     * 0.157
5           diet:time 1.38 15.17   0.603 5.01e-01       0.013
6      exercises:time 2.00 22.00  20.826 8.41e-06     * 0.274
7 diet:exercises:time 2.00 22.00  14.246 1.07e-04     * 0.147
>

6. 显著与不显著交互的处理逻辑（重点）

情况	说明	该怎么做
三阶交互显著	三因素相互影响	分别固定一个因素 (`group_by`) 看另两个因素
某二阶交互显著	仅这两个因素有交互	仅对该交互做分层分析
某二阶交互不显著	无需分层	只报告主效应
主效应显著	做事后比较（pairwise_t_test）	逐水平比较

例如：

# 如果 exercises:time 显著 → 分层
weightloss %>%
  group_by(exercises) %>%
  anova_test(dv = score, wid = id, within = time)

⚠️ 若交互不显著，不建议做分层（统计意义不大，只可作探索）。

7. 分层示例（simple effects）

# 固定 diet，分析 "运动 × 时间"
two.way <- weightloss %>%
  group_by(diet) %>%
  anova_test(dv = score, wid = id, within = c(exercises, time))
get_anova_table(two.way)

# 固定 diet + exercises，分析时间趋势
time.effect <- weightloss %>%
  group_by(diet, exercises) %>%
  anova_test(dv = score, wid = id, within = time)
get_anova_table(time.effect)

8. 事后检验

交互显著 → 主效应不能直接解释 → 必须分层分析；
分层显著 → 可以在该层内部做事后检验。
仅针对显著交互部分做事后比较：

pwc3 <- weightloss %>%
  group_by(diet, exercises) %>%
  pairwise_t_test(score ~ time, paired = TRUE, p.adjust.method = "bonferroni")
pwc3 <- pwc3 %>% add_xy_position(x = "exercises")
 
bxp3 +
  stat_pvalue_manual(pwc3, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov3, detailed = TRUE),
    caption = get_pwc_label(pwc3)
  )

五、为什么不需要显式建立模型？

传统写法（老派 aov）：

aov(score ~ treatment * time + Error(id/time), data = df)

现代写法（rstatix）：

anova_test(data = df, dv = score, wid = id, within = c(treatment, time))

区别在于：

aov 需要手动写出误差项；
anova_test 自动构造误差结构；
参数 wid = id 指明重复测量单位；
并自动进行球形度检验与校正。

所以重复测量 ANOVA 不需要显式模型，因为 wid 已定义“个体内重复”。

六、总结与建议

层次	函数写法	分析目的
单因素	`anova_test(dv=score, wid=id, within=time)`	检查时间效应
双因素	`anova_test(dv=score, wid=id, within=c(treatment,time))`	检查处理 × 时间交互
三因素	`anova_test(dv=score, wid=id, within=c(diet,exercises,time))`	检查三因素交互
分层分析	`group_by(...) %>% anova_test(...)`	检查简单效应
事后比较	`pairwise_t_test(..., paired=TRUE)`	检查具体差异

七分层分析（simple effects analysis）与事后检验（post-hoc test）。

一、基本区别概览

对比项	分层分析（Simple Effects Analysis）	事后检验（Post-hoc Test）
目的	探索交互效应的内部结构（即一个因素在另一个因素不同水平下的效应）	探索主效应显著后，具体哪些组之间存在差异
前提条件	通常在存在显著交互效应时进行	通常在主效应显著但无交互或交互不显著时进行
分析逻辑	“拆解交互” → 在一个因素的每个水平下，考察另一个因素的效应	“细化主效应” → 对同一因素不同水平进行两两比较
统计模型	group_by一个或多个因素后，再对剩下的因素做anova_test（如 group_by(A, B) → 检查 C 的效应）	pairwise_t_test()，在显著主效应的因素内部进行两两比较
反映的信息	哪个因素之间的交互作用导致总体交互显著；即条件性效应	哪些水平之间的均值差异导致主效应显著；即平均效应
图形表现	交互折线图（lines crossing/interacting）	分组箱线图或均值±误差条图

二、用三因素方差分析举例说明

假设我们研究的是 饮食类型（Diet） × 运动方式（Exercise） × 时间（Time） 对 减重效果（Score） 的影响：

1️ 如果三者有显著交互（Diet × Exercise × Time）：

说明“时间的减重趋势”取决于饮食和运动组合。
此时我们需要 分层分析（Simple Effects Analysis）：

例如，group_by(Diet, Exercise) 后在每个组合下做 anova_test(dv = score, within = time)，看看不同饮食+运动条件下，时间效应是否显著；
或者 group_by(Diet, Time) 检查在不同时间点下，不同饮食的效果是否一致。

这样能告诉我们：

“哪一组的时间变化最明显？”
“是否只有有氧+高蛋白饮食组在第3周显著下降？”

分层分析揭示的是 交互的来源和方向性。

2️ 如果交互效应不显著，但某个主效应显著（例如 time 显著）：

表示总体上“时间对减重的影响”存在显著差异，且这种影响与饮食或运动无关。
此时适合做 事后检验（Post-hoc Test）：

使用 pairwise_t_test(score ~ time, paired = TRUE, p.adjust.method = "bonferroni")
结果告诉我们哪个时间点之间的均值不同（比如 t1 与 t3 显著差异，t2 与 t3 不显著）

事后检验揭示的是 主效应中显著组间差异的具体来源。

lvzhaopeng的数字花园

探索

R语言那些事 统计系列（十） rstatix 重复测量方差分析

缘由

碎碎念

Ref

解读

一、准备阶段

1. 加载所需包

二、单因素重复测量 ANOVA

1. 数据来源

2. 转换成长格式

3. 描述统计

4. 可视化

5. 检查假设

(1) 异常值

(2) 正态性

6. 方差分析

7. 事后检验（配对 t 检验）

三、双因素重复测量 ANOVA

1. 数据准备

2. 描述统计

3. 可视化

4. 检查假设

5. 双因素重复测量 ANOVA

6. 事后检验：针对显著交互作用

四、三因素重复测量 ANOVA

1. 数据准备

2. 描述统计

3. 可视化

4. 检查假设

5. 三因素重复测量 ANOVA

6. 显著与不显著交互的处理逻辑（重点）

7. 分层示例（simple effects）

8. 事后检验

五、为什么不需要显式建立模型？

六、总结与建议

七 分层分析（simple effects analysis） 与 事后检验（post-hoc test）。

一、基本区别概览

二、用三因素方差分析举例说明

1️ 如果三者有显著交互（Diet × Exercise × Time）：

2️ 如果交互效应不显著，但某个主效应显著（例如 time 显著）：

关系图谱

目录

反向链接

R语言那些事统计系列（十） rstatix 重复测量方差分析

七分层分析（simple effects analysis）与事后检验（post-hoc test）。