一文搞定：现代化高级R包afex与重复测量方差分析（代码+可视化+结果解读+多重检验）
一文搞定：现代化高级R包afex与重复测量方差分析（代码+可视化+结果解读+多重检验）
医学和生信笔记
发布于 2026-04-16 15:29:37
发布于 2026-04-16 15:29:37
本文目录
- afex包简介
- 单因素重复测量
- 数据说明
- 加载R包
- 数据整理
- 进行方差分析
- 结果解读
- 两因素多水平（混合设计）
- 数据说明
- 数据整理
- 进行方差分析
- 结果解读
- 多重检验
- 为什么需要多重检验
- 什么是边际均值
- 三种层次的多重比较
- 可视化
- 使用afex_plot()
- 使用ggplot2自定义绘图
之前给大家详细介绍过如何使用R自带的经典函数aov做重复测量方差分析以及如何使用car::Anova做重复测量方差分析。这两种方法各有优缺点：aov函数简单易用，但在处理违反球对称假设时较为繁琐；car::Anova在自由度校正上更加灵活，但语法相对复杂。今天给大家介绍的afex包则综合了两者的优点，同时还能与emmeans、multcomp等包无缝衔接，完成事后多重比较，是目前进行重复测量方差分析的更优选择。
关于重复测量方差分析前的球对称检验，大家可以参考之前的文章R语言球对称检验。简单来说，球对称假设（Sphericity assumption）是重复测量方差分析的前提条件之一，它要求各个重复测量条件之间的差值具有相同的方差。一旦该假设被违反，F检验的结果会偏于乐观（即更容易出现假阳性），此时需要对自由度进行校正。
数据还是使用孙振球《医学统计学》第4版的数据，该数据大家可以加群免费下载，数据在QQ群文件（加群后别问数据在哪里，谢谢大家，数据在QQ群文件！）。
unsetunsetafex包简介unsetunset
afex是一个现代化的方差分析R包，名称来源于”Analysis of Factorial EXperiments”（析因实验分析），它整合并改进了car、ez等R包的核心功能。其主要功能包括：
- 三种ANOVA分析接口：提供
aov_car()、aov_ez()和aov_4()三个函数，用于估计标准方差分析模型，可灵活处理任意数量或任意组合的组内变量（within-subjects，即重复测量因素）和组间变量（between-subjects，即分组因素）。这三个函数拟合的是完全相同的统计模型，仅在模型写法（公式语法）上有所不同，输出结果完全一致，读者可根据自身习惯任选其一。 - 混合效应模型接口：
mixed()函数提供了一个基于lme4包（lmer/glmer）的混合效应模型分析接口，并能自动为所有固定效应项（包括主效应和交互作用）计算p值，弥补了lme4本身不直接输出p值的不足。 - 结果可视化：
afex_plot()函数可直接基于方差分析结果进行可视化，前景展示估计的边际均值（estimated marginal means）及其置信区间，背景呈现原始数据的分布，能够直观反映各组/各时间点的整体趋势与个体差异。 - 与emmeans包的无缝衔接：所有
afex模型对象（ANOVA模型和混合模型均可）都可以直接传递给emmeans包，用于进行事后检验（post-hoc tests）、计划对比（planned contrasts）或其他均值比较，极大简化了后续分析的流程。
unsetunset单因素重复测量unsetunset
数据说明
本节使用孙振球《医学统计学》第4版表12-3的数据。该数据只有一组受试者（即没有分组变量），8名受试者分别在4个时间点（t0、t45、t90、t135，单位：分钟）测量了血糖值。需要特别说明的是，这是一个不符合球对称假设的数据，因此在分析时需要对自由度进行校正，这也是我们在后文中选择afex包（而非直接用aov函数）的重要原因之一。
先读取数据：
data12_3b <- foreign::read.spss(
"../datasets/表12-3重复测量ANOVA.sav",
to.data.frame = T,
reencode = "utf-8"
)
## re-encoding from utf-8
str(data12_3b)
## 'data.frame': 8 obs. of 4 variables:
## $ t0 : num 5.32 5.32 5.94 5.49 5.71 6.27 5.88 5.32
## $ t45 : num 5.32 5.26 5.88 5.43 5.49 6.27 5.77 5.15
## $ t90 : num 4.98 4.93 5.43 5.32 5.43 5.66 5.43 5.04
## $ t135: num 4.65 4.7 5.04 5.04 4.93 5.26 4.93 4.48
## - attr(*, "variable.labels")= Named chr(0)
## ..- attr(*, "names")= chr(0)
## - attr(*, "codepage")= int 936
head(data12_3b)
## t0 t45 t90 t135
## 1 5.32 5.32 4.98 4.65
## 2 5.32 5.26 4.93 4.70
## 3 5.94 5.88 5.43 5.04
## 4 5.49 5.43 5.32 5.04
## 5 5.71 5.49 5.43 4.93
## 6 6.27 6.27 5.66 5.26
加载R包
library(afex)
## Warning: package 'afex' was built under R version 4.5.2
## Loading required package: lme4
## Loading required package: Matrix
## ************
## Welcome to afex. For support visit: http://afex.singmann.science/
## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - Get and set global package options with: afex_options()
## - Set sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************
##
## Attaching package: 'afex'
## The following object is masked from 'package:lme4':
##
## lmer
数据整理
afex包要求数据采用长格式（long format），即每行代表一个观测值（一个受试者在一个时间点的测量结果），而不是常见的宽格式（每行代表一个受试者，每列代表一个时间点）。
此外，afex的ANOVA函数还要求明确指定一个被试标识列（即id列），用于区分不同受试者。这是因为：如果某位受试者在同一实验条件下有多次测量值（即每格多于一个观测），afex默认会将这些重复值聚合（取平均）后再进行分析，而id列正是实现这一聚合的依据。
我们用dplyr和tidyr包将数据转换为长格式，并添加id列：
library(dplyr)
library(tidyr)
data12_3b_long <- data12_3b %>%
# 新增id列：用行号作为受试者编号，并转为因子型
mutate(id = row_number(), id = factor(id)) %>%
# 将4个时间点从宽格式转为长格式：列名存入"times"，数值存入"bg"
pivot_longer(cols = 1:4, names_to = "times", values_to = "bg") %>%
# 将时间点设为有序因子，确保t0→t45→t90→t135的顺序正确
mutate(times = factor(times, levels = c("t0", "t45", "t90", "t135")))
head(data12_3b_long)
## # A tibble: 6 × 3
## id times bg
## <fct> <fct> <dbl>
## 1 1 t0 5.32
## 2 1 t45 5.32
## 3 1 t90 4.98
## 4 1 t135 4.65
## 5 2 t0 5.32
## 6 2 t45 5.26
转换后，数据变成了 8 × 4 = 32 行，每行记录一名受试者在一个时间点的血糖值，这正是afex所要求的格式。
进行方差分析
afex提供了三种等价的模型写法，以下两种写法给出的结果完全相同，任选一种即可：
# 写法一：公式写法（类似aov()中Error()的语法）
# aov_car(bg ~ times + Error(id/times), data = data12_3b_long)
# 写法二：分别指定组间变量和组内变量（语义更清晰）
f <- aov_ez(
"id",
"bg",
data12_3b_long,
between = NULL, # 该数据没有分组变量（组间变量）
within = "times"# 时间点是重复测量的组内变量
)
f
## Anova Table (Type 3 tests)
##
## Response: bg
## Effect df MSE F ges p.value
## 1 times 1.58, 11.09 0.02 79.14 *** .514 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
结果解读
无论数据是否符合球对称假设，afex都会默认采用Greenhouse-Geisser（GG）校正法自动调整自由度和p值（输出中标注为Sphericity correction method: GG），从而避免因球对称假设违反而产生的I型错误膨胀。
输出结果中各指标的含义如下：
times的p值 < 0.05：说明不同时间点之间的血糖值存在显著差异，即存在显著的时间主效应。换句话说，受试者的血糖水平随时间推移发生了统计学上有意义的变化。- **
df(1.58, 11.09)**：这是经过GG校正后的自由度，括号内依次为分子自由度（效应自由度）和分母自由度（误差自由度）。未校正时，4个时间点的分子自由度本应为3，GG校正将其缩小至1.58，自由度的缩小使p值更为保守（即更难达到显著），从而更好地控制了假阳性风险。 MSE（Mean Square Error，均方误差）：误差均方，代表无法被模型解释的随机变异量，是计算F值的分母组成部分。MSE越小，说明数据的组内变异越小，统计检验的精度越高。- **
F = 79.14**：方差分析的F统计量（此处为GG校正后的结果）。F值越大，说明组间变异相对于组内变异越显著。F值 = 效应均方 / 误差均方，直观理解就是”信号强度”与”噪声强度”之比。 - **
ges = 0.514：Generalized Eta-Squared（广义eta平方），是afex默认报告的效应量指标。效应量用于衡量统计效应的实际大小，不受样本量的影响。ges = 0.514表示约51.4%**的总变异可以由”时间”这一因素来解释，属于大效应（通常 ηG2 > 0.14 即为大效应）。
如果你希望得到与car::Anova()输出格式相同的详细结果（包括球对称检验Mauchly’s test以及GG、HF两种校正结果），可以使用summary()：
summary(f)
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 914.53 1 2.53170 7 2528.623 3.221e-10 ***
## times 2.96 3 0.26182 21 79.141 1.304e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## times 0.06273 0.0082069
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## times 0.52827 5.368e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## times 0.6573076 2.890609e-08
也可以使用anova(f)，将结果输出为data.frame格式，方便后续导出或进一步处理：
# 这个结果是数据框格式，方便导出到Excel等
anova(f)
## Anova Table (Type 3 tests)
##
## Response: bg
## num Df den Df MSE F ges Pr(>F)
## times 1.5848 11.094 0.023601 79.141 0.51447 5.368e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
unsetunset两因素多水平（混合设计）unsetunset
数据说明
本节使用孙振球《医学统计学》第4版例12-3的数据。研究设计如下：将手术要求基本相同的15名患者随机分为3组，分别采用A、B、C三种麻醉诱导方法（这是组间变量：不同患者接受不同的处理）；同时在T0（诱导前）、T1、T2、T3、T4共5个时相测量患者的收缩压（这是组内变量：同一患者在多个时间点被反复测量）。这种同时含有组间和组内因素的设计，称为混合设计重复测量方差分析（mixed-design repeated measures ANOVA）。
data12_3 <- foreign::read.spss(
"../datasets/例12-03.sav",
to.data.frame = T,
reencode = "gbk"
)
## re-encoding from gbk
str(data12_3)
## 'data.frame': 15 obs. of 7 variables:
## $ No : num 1 2 3 4 5 6 7 8 9 10 ...
## $ group: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 2 2 2 2 2 ...
## $ t0 : num 120 118 119 121 127 121 122 128 117 118 ...
## $ t1 : num 108 109 112 112 121 120 121 129 115 114 ...
## $ t2 : num 112 115 119 119 127 118 119 126 111 116 ...
## $ t3 : num 120 126 124 126 133 131 129 135 123 123 ...
## $ t4 : num 117 123 118 120 126 137 133 142 131 133 ...
## - attr(*, "variable.labels")= Named chr [1:7] "序号" "组别" "" "" ...
## ..- attr(*, "names")= chr [1:7] "No" "group" "t0" "t1" ...
head(data12_3)
## No group t0 t1 t2 t3 t4
## 1 1 A 120 108 112 120 117
## 2 2 A 118 109 115 126 123
## 3 3 A 119 112 119 124 118
## 4 4 A 121 112 119 126 120
## 5 5 A 127 121 127 133 126
## 6 6 B 121 120 118 131 137
数据共7列：第1列是患者编号（No），第2列是诱导方法（group，即A/B/C三组），第3至7列是5个时间点（T0～T4）的收缩压。
数据整理
与单因素重复测量相同，先将数据从宽格式转为长格式：
library(dplyr)
library(tidyr)
data12_3_long <- data12_3 %>%
# 将患者编号转为因子型，作为被试id
mutate(No = factor(No)) %>%
# 保留No和group列，将其余时间点列转为长格式
pivot_longer(
cols = -c("No", "group"),
names_to = "times",
values_to = "hp"
) %>%
mutate(times = factor(times))
head(data12_3_long)
## # A tibble: 6 × 4
## No group times hp
## <fct> <fct> <fct> <dbl>
## 1 1 A t0 120
## 2 1 A t1 108
## 3 1 A t2 112
## 4 1 A t3 120
## 5 1 A t4 117
## 6 2 A t0 118
进行方差分析
# 经典的公式写法（等价）：
# f <- aov_car(hp ~ times * group + Error(No/times), data = data12_3_long)
# 推荐写法：明确指定组间变量（between）和组内变量（within）
f <- aov_ez(
"No",
"hp",
data12_3_long,
between = "group", # 组间变量：麻醉诱导方法（A/B/C）
within = "times"# 组内变量：测量时相（T0~T4）
)
## Contrasts set to contr.sum for the following variables: group
f
## Anova Table (Type 3 tests)
##
## Response: hp
## Effect df MSE F ges p.value
## 1 group 2, 12 78.87 5.78 * .430 .017
## 2 times 2.71, 32.58 8.08 106.56 *** .659 <.001
## 3 group:times 5.43, 32.58 8.08 19.10 *** .409 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG
结果解读
输出结果包含三个效应的检验：
group（麻醉诱导方法的主效应）：F = 5.78，p < 0.05，说明三种麻醉诱导方法对患者收缩压的总体影响存在显著差异。也就是说，不考虑时间因素，A、B、C三组患者的平均收缩压有显著不同。times（时间的主效应）：说明不考虑分组因素，患者的收缩压在5个时相之间存在显著变化（具体结果可由输出查看）。group:times（交互作用）：F = 19.10，p < 0.01，这是最重要的结果。交互作用显著意味着：三种麻醉诱导方法对收缩压的影响随时间推移呈现出不同的变化趋势——例如，A组的血压可能先降后升，而B组则持续下降。简单来说，“哪种方法效果更好”这个问题的答案会随时间点的不同而改变。当交互作用显著时，通常需要进一步做简单效应分析（即分组或分时间点的事后比较），才能得出有临床意义的结论。
unsetunset多重检验unsetunset
为什么需要多重检验
方差分析的F检验只能告诉我们”各组/各时间点之间存在差异”，但无法告诉我们”哪两组之间有差异”。当有3个及以上的水平需要两两比较时，如果直接反复做t检验，会导致I型错误（假阳性）急剧膨胀：例如做10次比较时，犯至少一次错误的概率高达 1 − (1 − 0.05)10 ≈ 40%。因此需要采用多重检验校正（multiple comparison correction）的方法来控制整体错误率。
关于重复测量数据多重检验的更多背景，可参考重复测量方差分析的多重检验。afex包的所有分析结果均可直接传递给emmeans包进行事后检验，下面逐步演示。
什么是边际均值
在进行多重比较之前，先理解一个核心概念——边际均值（marginal means），在emmeans包中也称为估计边际均值（estimated marginal means，EMMs）。
边际均值是指在控制（平衡）其他因素的影响后，某一特定因素各水平上的预测平均值。它与我们直接计算的”样本均值”有所不同：
- 样本均值：直接对观测数据求平均，受各组样本量不均衡、其他协变量等因素的影响。
- 边际均值：基于统计模型推导，通过对其他因素（如时间点、协变量）进行边际化（平均化）处理，得到在”理想平衡条件下”各组的平均估计值。
举个例子：在我们的麻醉数据中，A组的边际均值是指：对A组患者在所有时间点（T0～T4）上的预测收缩压取平均，代表了A组”整体上”的收缩压水平，而不受各时间点样本量差异的干扰。边际均值更能公平地反映各组的真实效应，因此是事后多重比较的基础。
三种层次的多重比较
加载emmeans包，逐步演示三种层次的多重比较：
library(emmeans)
## Warning: package 'emmeans' was built under R version 4.5.2
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
第一层：比较不同组（麻醉方法）之间的差异
首先计算每个麻醉诱导方法组的边际均值（此时已对所有时间点取均值）：
# emmeans(模型对象, 指定要比较的因素)
# ~ group 表示我们关注的是"group"这个因素上各水平的均值
m1 <- emmeans(f, ~group)
m1
## group emmean SE df lower.CL upper.CL
## A 120 1.78 12 116 124
## B 124 1.78 12 121 128
## C 128 1.78 12 124 132
##
## Results are averaged over the levels of: times
## Confidence level used: 0.95
输出结果中，emmean列即为各组的边际均值，SE为标准误，lower.CL和upper.CL为95%置信区间。
然后进行A、B、C三组之间的两两比较：
# pairs()函数默认使用Tukey法对p值进行多重校正
pairs(m1)
## contrast estimate SE df t.ratio p.value
## A - B -4.80 2.51 12 -1.911 0.1780
## A - C -8.52 2.51 12 -3.392 0.0137
## B - C -3.72 2.51 12 -1.481 0.3338
##
## Results are averaged over the levels of: times
## P value adjustment: tukey method for comparing a family of 3 estimates
# 上述代码等价于：
# contrast(m1, method = "pairwise")
pairs(m1)默认采用Tukey法进行多重校正，输出每对比较的估计差值（estimate）、标准误（SE）、自由度（df）、t统计量及校正后的p值（p.value）。Tukey法是最常用的事后检验方法之一，在各组样本量相等时表现最优，能在保证检验效能的同时有效控制族错误率（familywise error rate）。
使用multcomp包获得更强大的p值校正
multcomp包（Multiple Comparisons using Linear Models，基于线性模型的多重比较）是R语言中进行多重比较和同步推断的专业工具包，由Frank Bretz、Torsten Hothorn、Peter Westfall等统计学家开发。
相比emmeans内置的校正方法，multcomp的核心优势在于：
- 更丰富的校正方法：除了Tukey、Bonferroni等经典方法外，还提供了
"free"、"Westfall"等更先进的校正策略，能够更好地利用检验之间的相关性结构，在控制错误率的同时保留更高的检验效能。 - 同步置信区间：
multcomp能计算与多重校正后的p值完全配套的同步置信区间（simultaneous confidence intervals），而不是单独的逐对置信区间，在理论上更加严格。 - 广泛适用性：可用于线性模型、广义线性模型、混合模型等多种统计模型框架的多重比较。
下面的代码将emmeans的比较结果转入multcomp，并使用"free"方法进行p值校正：
library(multcomp)
summary(as.glht(pairs(m1)), test = adjusted("free"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## A - B == 0 -4.800 2.512 -1.911 0.1392
## A - C == 0 -8.520 2.512 -3.392 0.0137 *
## B - C == 0 -3.720 2.512 -1.481 0.1644
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- free method)
这行代码分三步执行，我们逐一拆解：
第一步：pairs(m1)
对m1（A、B、C三组的边际均值）进行所有两两配对比较，生成3组对比（A-B、A-C、B-C）的检验结果。此时结果仍在emmeans框架内，pairs()默认使用Tukey法进行校正——但请注意，这个默认校正会在下一步被覆盖，最终p值完全由summary()中指定的adjusted("free")方法决定。
第二步：as.glht(pairs(m1))
as.glht()函数将emmeans生成的配对比较结果，转换为multcomp包识别的广义线性假设检验（general linear hypotheses test，glht）对象。这一步是连接emmeans与multcomp的桥梁，使我们可以在multcomp的框架下对同样的比较应用更灵活的校正方法。
第三步：summary(..., test = adjusted("free"))
对转换后的glht对象进行汇总，并指定p值校正方法为"free"。
"free"方法（也称Free Step-Down Resampling方法）是一种基于多元检验统计量相关性结构的逐步校正方法。其核心思想是：在进行多个比较时，如果各检验之间存在正相关（例如同一组数据中的多个比较往往是相关的），则Bonferroni校正会过于保守（p值被压得过低，导致本可显著的比较变得不显著）。"free"方法通过利用这种相关性，给出比Bonferroni更宽松、比不校正更保守的调整p值，在控制族错误率的前提下尽可能保留检验效能。需要说明的是，"free"与Tukey法针对的都是强族错误率（strong FWER）的控制，两者并无简单的优劣之分；"free"的优势主要体现在比较数量较多、检验之间相关性较强的场景下，此时它通常比Bonferroni有更高的检验效能，结果也可能比Tukey更宽松——但具体效果取决于数据的相关结构。
第二层：分组内的时间点两两比较
计算在每个麻醉方法组内，各时间点各自的边际均值（即在每个组别中，分别估计T0～T4五个时间点的均值，组别之间互不混合）：
# "times" 是我们关注的因素，by = "group" 表示分组进行
m2 <- emmeans(f, "times", by = "group")
# 等价写法：emmeans(f, ~ times | group)
m2
## group = A:
## times emmean SE df lower.CL upper.CL
## t0 121 1.72 12 117 125
## t1 112 2.21 12 108 117
## t2 118 2.08 12 114 123
## t3 126 2.20 12 121 131
## t4 121 1.76 12 117 125
##
## group = B:
## times emmean SE df lower.CL upper.CL
## t0 121 1.72 12 117 125
## t1 120 2.21 12 115 125
## t2 118 2.08 12 113 123
## t3 128 2.20 12 123 133
## t4 135 1.76 12 131 139
##
## group = C:
## times emmean SE df lower.CL upper.CL
## t0 126 1.72 12 122 130
## t1 123 2.21 12 118 128
## t2 119 2.08 12 114 123
## t3 143 2.20 12 138 147
## t4 131 1.76 12 127 134
##
## Confidence level used: 0.95
这里by = "group"的含义是：分别在A组、B组、C组内，计算每个时间点（T0～T4）的边际均值。输出会呈现3组×5时间点 = 15个估计值。
然后分别在每个组内对5个时间点进行两两比较：
# 每个组内分别进行时间点的两两比较，各自做Tukey校正
pairs(m2)
## group = A:
## contrast estimate SE df t.ratio p.value
## t0 - t1 8.6 1.490 12 5.772 0.0007
## t0 - t2 2.6 1.320 12 1.964 0.3382
## t0 - t3 -4.8 2.060 12 -2.333 0.2001
## t0 - t4 0.2 1.680 12 0.119 0.9999
## t1 - t2 -6.0 0.913 12 -6.573 0.0002
## t1 - t3 -13.4 1.060 12 -12.624 <0.0001
## t1 - t4 -8.4 1.530 12 -5.507 0.0010
## t2 - t3 -7.4 1.460 12 -5.066 0.0021
## t2 - t4 -2.4 1.340 12 -1.789 0.4223
## t3 - t4 5.0 1.630 12 3.062 0.0618
##
## group = B:
## contrast estimate SE df t.ratio p.value
## t0 - t1 1.4 1.490 12 0.940 0.8760
## t0 - t2 3.2 1.320 12 2.417 0.1762
## t0 - t3 -7.0 2.060 12 -3.402 0.0347
## t0 - t4 -14.0 1.680 12 -8.317 <0.0001
## t1 - t2 1.8 0.913 12 1.972 0.3345
## t1 - t3 -8.4 1.060 12 -7.914 <0.0001
## t1 - t4 -15.4 1.530 12 -10.096 <0.0001
## t2 - t3 -10.2 1.460 12 -6.983 0.0001
## t2 - t4 -17.2 1.340 12 -12.820 <0.0001
## t3 - t4 -7.0 1.630 12 -4.287 0.0076
##
## group = C:
## contrast estimate SE df t.ratio p.value
## t0 - t1 3.2 1.490 12 2.148 0.2625
## t0 - t2 7.6 1.320 12 5.740 0.0007
## t0 - t3 -16.4 2.060 12 -7.971 <0.0001
## t0 - t4 -4.4 1.680 12 -2.614 0.1293
## t1 - t2 4.4 0.913 12 4.820 0.0031
## t1 - t3 -19.6 1.060 12 -18.465 <0.0001
## t1 - t4 -7.6 1.530 12 -4.982 0.0024
## t2 - t3 -24.0 1.460 12 -16.432 <0.0001
## t2 - t4 -12.0 1.340 12 -8.944 <0.0001
## t3 - t4 12.0 1.630 12 7.348 <0.0001
##
## P value adjustment: tukey method for comparing a family of 5 estimates
输出会分三个部分（group = A、group = B、group = C），每部分包含5个时间点的两两比较（共C52 = 10对）。这种分析方式可以回答：“在A组内，T0和T1的收缩压有没有显著差异？”等问题，对临床解读交互作用非常有帮助。
第三层：所有组别×时间点组合的两两比较
如果需要进行更全面的比较（例如”A组T0”与”B组T1”是否有差异），可以计算所有组合下的边际均值：
# 同时考虑times和group两个因素的所有组合
m3 <- emmeans(f, c("times", "group"))
# 等价写法：emmeans(f, ~ times * group)
m3
## times group emmean SE df lower.CL upper.CL
## t0 A 121 1.72 12 117 125
## t1 A 112 2.21 12 108 117
## t2 A 118 2.08 12 114 123
## t3 A 126 2.20 12 121 131
## t4 A 121 1.76 12 117 125
## t0 B 121 1.72 12 117 125
## t1 B 120 2.21 12 115 125
## t2 B 118 2.08 12 113 123
## t3 B 128 2.20 12 123 133
## t4 B 135 1.76 12 131 139
## t0 C 126 1.72 12 122 130
## t1 C 123 2.21 12 118 128
## t2 C 119 2.08 12 114 123
## t3 C 143 2.20 12 138 147
## t4 C 131 1.76 12 127 134
##
## Confidence level used: 0.95
此时共有 5 × 3 = 15 个水平组合，两两比较会产生 C152 = 105 对结果，通常只在有明确研究假设时才进行全面比较，否则结果过多，临床解读困难：
# 结果共105行，数量较多，仅在有明确需要时运行
pairs(m3)
unsetunset可视化unsetunset
使用afex_plot()
afex自带的afex_plot()函数可以直接对接方差分析结果进行可视化，前景展示的是各组各时间点的边际均值（而非原始均值），背景展示个体数据的分布。这种设计可以帮助读者同时看到总体趋势和个体差异：
# x指定横轴变量（时间点），trace指定用不同线条区分的变量（诱导方法）
# error = "none" 表示不显示误差线（可改为"mean_se"等）
afex_plot(f, x = "times", trace = "group", error = "none")
使用ggplot2自定义绘图
afex_plot()本质上是对ggplot2的封装，如果希望对图形进行更灵活的自定义（如修改颜色、添加标注、调整主题等），完全可以用ggplot2实现同样乃至更精美的图形：
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.2
library(dplyr)
# 先计算各组各时间点的样本均值（用于绘制折线）
mm <- data12_3_long %>%
summarise(hp = mean(hp), .by = c(times, group))
ggplot() +
# 绘制各组均值折线（不同组用不同线型区分）
geom_line(
data = mm,
aes(x = times, y = hp, group = group, linetype = group)
) +
# 绘制各组均值点（不同组用不同形状区分）
geom_point(data = mm, aes(x = times, y = hp, shape = group)) +
# 以半透明散点展示原始数据（背景层）
geom_point(
data = data12_3_long,
aes(x = times, y = hp, shape = group),
alpha = 0.2
) +
labs(
x = "测量时相",
y = "收缩压 (mmHg)",
linetype = "诱导方法",
shape = "诱导方法"
) +
theme_bw()
对比两图可以发现，afex_plot()展示的是边际均值（模型估计值），而ggplot2绘制的是样本均值（直接计算值）。在样本量均衡的情况下两者差异不大，但若设计不均衡（各组样本量不等），推荐使用边际均值进行展示，结果更具参考价值。
参考资料:
- https://cran.r-project.org/web/packages/afex/vignettes/afex_anova_example.html
- Singmann, H., & Kellen, D. (2020). An Introduction to Mixed Models for Experimental Psychology. In D. Spieler & E. Schumacher (Eds.),New Methods in Cognitive Psychology. Psychology Press.
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2026-04-08，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号