重复测量数据在临床研究和实验科学里出现的频率极高——同一个受试者在不同时间点被反复观测，比如术后第1天、第3天、第7天的某项指标，或者用药前后多个随访节点的量表得分。这类数据的特点是观测值之间不独立，同一个体内部的多次测量存在相关性，如果直接拿去做普通的独立样本t检验或者单因素方差分析，p值大概率是错的，结论也就站不住脚。我见过不少做临床数据分析的朋友，数据明明是重复测量的结构，却硬生生把每个时间点拆开做多次t检验，结果要么假阳性泛滥，要么把本来显著的处理效应给淹没了。这篇内容就围绕“基于R进行两组患者重复测量数据的分析”这个主题，把重复测量方差分析的完整思路、R语言落地细节、以及我踩过的坑，一次性讲透。适合有一定R基础、正在处理纵向随访数据或者交叉实验数据的从业者参考，也适合刚接触重复测量设计、想搞清楚“为什么不能拆开做t检验”的研究人员从头看起。

1. 重复测量方差分析到底在解决什么问题

1.1 从数据结构说起：为什么普通方差分析不够用

先把数据结构摆清楚。假设我们做了一个两组对照的临床观察，A组用新方案，B组用常规方案，每组各若干名患者，在基线、第4周、第8周、第12周四个时间点测量同一个连续型指标。数据排成“长格式”之后，每一行是“某患者在某时间点的一次测量”，包含患者ID、分组、时间点、测量值四个核心字段。这种结构下，同一个患者贡献了4行数据，这4个观测值天然相关——一个基线值偏高的患者，后续几个时间点大概率也偏高。普通单因素方差分析的前提是观测独立，这个前提在这里直接被打破了。

如果无视相关性，把4个时间点当成4个独立组来做方差分析，误差项里混入了“患者个体差异”这部分变异，导致误差均方被高估，检验效能下降。更麻烦的是，时间点之间的比较会失真，因为你没法区分“时间效应”和“个体差异”。重复测量方差分析的核心思路，就是把总变异拆成三块：组间变异（处理效应）、组内变异（时间效应）、以及交互效应（处理×时间），同时把“个体差异”从误差项里单独剥离出来，让误差估计更干净。这就是它比普通方差分析更适合纵向数据的根本原因。

用生活化的类比：假设你想比较两个班级学生的考试成绩变化，如果每个学生只考一次，那是独立样本问题；但如果是同一批学生期中、期末各考一次，学生自身的水平差异就是一个巨大的干扰项。重复测量方差分析相当于先把“这个学生本来就行/本来就不行”这部分变异扣掉，再去看“时间带来的变化”和“班级之间的差异”，结论自然更靠谱。

1.2 两组设计的特殊性：组间、组内、交互三重效应

两组重复测量设计里，我们关心的效应有三个层面。第一是组间效应，也就是A组和B组整体上有没有差异，这个效应在统计上等价于“把每个患者的多次测量取平均，再做两组比较”，但它用的是合并后的误差项，比直接做t检验更稳健。第二是组内效应，也就是时间主效应，所有患者合在一起看，指标随时间有没有显著变化。第三是交互效应，这是最有临床价值的一项——两组随时间变化的趋势是否不同。举个例子，如果新方案组在第8周之后指标快速下降，而常规组下降平缓，交互效应就会显著，说明处理方案的效果依赖于时间，这比单纯看“两组有没有差异”信息量大得多。

很多初学者只盯着组间p值看，忽略了交互效应，这是很可惜的。临床上真正有意义的往往是“处理组起效更快”或者“处理组在后期维持得更好”这类趋势差异，而这些恰恰体现在交互项上。所以在R里跑重复测量方差分析时，交互项的显著性一定要单独拎出来解读，不能只看主效应。

1.3 球形假设：这个前提不满足会怎样

重复测量方差分析有一个绕不开的假设——球形假设（sphericity）。简单说，它要求任意两个时间点之间差值的方差相等。四个时间点的话，就有6对差值，这6个差值的方差理论上应该差不多。如果球形假设被违反，组内效应和交互效应的F检验自由度需要校正，否则p值偏小，容易假阳性。

球形假设在时间点较多、或者时间间隔不均匀时特别容易出问题。我的经验是：三个时间点以内，球形假设通常还过得去；四个及以上时间点，最好默认做校正。R里car包的Anova()函数配合idata参数可以自动输出GG校正（Greenhouse-Geisser）和HF校正（Huynh-Feldt）的结果，直接看校正后的p值就行。这里有个细节：球形假设只影响组内效应和交互效应，组间效应不受影响，因为组间比较用的是合并的个体均值，不涉及时间点之间的差值结构。这一点很多人搞混，以为球形假设不满足整个分析就废了，其实不是。

2. R语言落地的整体方案与工具选型

2.1 宽格式还是长格式：数据整形的关键决策

在R里做重复测量方差分析，数据格式的选择直接影响后续代码怎么写。宽格式是每个患者一行，不同时间点的测量值各占一列，比如id, group, t1, t2, t3, t4。长格式是每个患者每个时间点一行，比如id, group, time, value。两种格式各有用途：宽格式适合用aov()配合Error()项做经典的重复测量方差分析，长格式适合用混合效应模型或者car::Anova()配合idata参数。

我个人的建议是：如果时间点固定且不多（比如四个随访节点），用宽格式配合aov()最直接，代码短、输出清晰；如果时间点较多、或者有缺失值、或者想顺便做趋势分析，长格式更灵活。实际项目里我通常两份都准备——用tidyr::pivot_longer()和pivot_wider()互相转换，几行代码的事。需要注意的是，宽格式转长格式时列名要规范，比如week0, week4, week8, week12，转完之后时间因子要显式设定水平顺序，否则R默认按字母排序，week12会排在week4前面，图和分析全乱套。

2.2 核心函数对比：aov、Anova与混合模型

R里能做重复测量方差分析的函数不止一个，选哪个取决于你的设计复杂度和对缺失值的容忍度。stats::aov()配合Error()项是最经典的做法，语法是aov(value ~ group * time + Error(id/time))，其中Error(id/time)告诉R“个体是随机效应，时间是个体内因子”。这个函数的好处是输出直接给出组间、组内、交互三张方差分析表，球形假设的校正需要额外用summary()配合Mauchly检验来看。

car::Anova()是另一个常用选择，它需要把宽格式数据转成mlm对象（多元线性模型），然后用idata指定个体内因子结构。它的优势是自动输出球形校正结果，而且对不平衡设计更稳健。缺点是数据准备稍麻烦，需要构造idata数据框。

如果数据有缺失值（临床随访里太常见了），上面两个函数都会直接丢掉含缺失的整行，样本量缩水严重。这时候混合效应模型（lme4::lmer()或nlme::lme()）是更优解，它用限制性极大似然估计，能利用部分信息，不需要完整案例。不过混合模型的方差结构设定更复杂，随机斜率、随机截距的选择需要结合设计来定，不是一句两句能说完的。这篇内容聚焦经典重复测量方差分析，混合模型留到后续单独展开。

2.3 必备包清单与版本注意事项

跑完整套分析，我常用的包有这么几个：car负责Anova()和球形校正，emmeans负责事后两两比较和简单效应分析，tidyr和dplyr负责数据整形，ggplot2负责可视化，rstatix提供了一些便捷的管道式统计函数。版本上要注意，car包3.0之后的Anova()对idata参数的校验更严格，idata的行数必须和个体内因子的水平数完全一致，顺序也要对上，否则报错信息很隐晦。emmeans在1.7版本之后对重复测量模型的mode参数有调整，做简单效应时记得显式指定mode = "df.error"或者根据模型类型选对自由度近似方式，不然置信区间会偏窄。

提示：装包时建议锁定版本，尤其是car和emmeans，这两个包更新频繁，不同版本之间输出格式和默认行为有差异，团队协作时最好在项目里记录sessionInfo()。

3. 完整实操流程：从数据准备到结果解读

3.1 模拟一份两组四时点的临床数据

为了把流程跑通，我先构造一份模拟数据。两组各20名患者，四个时间点，指标值设定为：A组基线均值50，随时间逐步下降到第12周的38；B组基线均值50，下降到第12周的45。个体随机效应标准差设为5，残差标准差设为3。这样设定下，交互效应应该是显著的，因为两组下降斜率不同。

set.seed(2024) n_per_group <- 20 times <- c(0, 4, 8, 12) groups <- rep(c("A", "B"), each = n_per_group) id <- 1:(2 * n_per_group) # 个体随机截距 rand_intercept <- rnorm(2 * n_per_group, 0, 5) # 构造长格式 dat_long <- expand.grid(id = id, time = times) dat_long$group <- groups[match(dat_long$id, id)] dat_long$value <- 50 + rand_intercept[dat_long$id] + ifelse(dat_long$group == "A", -1.0 * dat_long$time, -0.4 * dat_long$time) + rnorm(nrow(dat_long), 0, 3) dat_long$time <- factor(dat_long$time, levels = times) dat_long$group <- factor(dat_long$group, levels = c("A", "B"))

这段代码里有两个细节值得说。第一，time转因子时显式指定了levels = times，保证顺序是0、4、8、12，而不是默认的字符串排序。第二，个体随机截距用rand_intercept[dat_long$id]按id索引，确保同一个体的四次测量共享同一个随机截距，这是模拟重复测量数据的关键，漏了这一步数据就退化成独立观测了。

转成宽格式用于aov()分析：

library(tidyr) dat_wide <- pivot_wider(dat_long, names_from = time, values_from = value, names_prefix = "t")

转完之后列名是t0, t4, t8, t12，每个患者一行，配合group列，正好是aov()需要的结构。

3.2 用aov加Error项跑经典重复测量方差分析

宽格式数据准备好之后，核心代码就一行：

library(car) mod_aov <- aov(cbind(t0, t4, t8, t12) ~ group, data = dat_wide) summary(mod_aov)

注意这里用的是cbind()把四个时间点绑成一个矩阵作为因变量，group是组间因子。这个写法下，summary()输出的是组间效应的检验结果。要拿到组内效应和交互效应，需要把模型转成mlm对象再用Anova()：

mod_mlm <- lm(cbind(t0, t4, t8, t12) ~ group, data = dat_wide) idata <- data.frame(time = factor(c("t0", "t4", "t8", "t12"), levels = c("t0", "t4", "t8", "t12"))) Anova(mod_mlm, idata = idata, idesign = ~time, type = "III")

idesign = ~time指定个体内因子是time，idata的行数和顺序必须和cbind()里列的排列完全一致，这是最容易出错的地方。输出里会包含time的球形检验结果（Mauchly's test），以及未校正和GG/HF校正后的F值和p值。组间效应看group那一行，交互效应看group:time那一行。

实测这份模拟数据，交互项的GG校正p值应该在0.05以下，说明两组的时间趋势确实不同。组间主效应可能不显著，因为两组在第12周的均值差距被基线到中期的重叠稀释了——这正好说明为什么不能只看组间p值。

3.3 球形检验与校正结果的取舍

Anova()输出里会有一行Mauchly's test for sphericity，给出W统计量和p值。如果p < 0.05，说明球形假设被违反，需要看校正后的结果。GG校正（Greenhouse-Geisser）的epsilon值通常偏保守，HF校正（Huynh-Feldt）偏宽松。我的取舍原则是：epsilon在0.75以上，用HF校正；0.75以下，用GG校正。这个阈值不是绝对的，但实践中比较稳妥。

还有一个细节：球形检验本身对样本量敏感，小样本时检验效能低，可能漏检；大样本时又过于敏感，轻微偏离就报显著。所以不能完全依赖p值判断，结合epsilon值一起看更靠谱。如果epsilon接近1（比如0.95以上），即使Mauchly检验显著，未校正结果也基本可用。

注意：球形假设只约束组内效应和交互效应，组间效应的F检验不需要校正。解读结果时别把组间p值也拿去做校正，那是画蛇添足。

3.4 事后比较：简单效应分析与两两对比

交互效应显著之后，下一步是拆解——在每个时间点上比较两组，或者在每组内部比较各时间点。这就是简单效应分析。用emmeans包可以很干净地完成：

library(emmeans) emm <- emmeans(mod_mlm, ~ group | time, mode = "df.error") pairs(emm)

~ group | time的意思是“在每个time水平上比较group”，输出的是四个时间点各自的两组对比。反过来，如果想在每组内部比较时间点，写成~ time | group。mode = "df.error"指定用误差项自由度做近似，对于经典的重复测量设计这是合适的；如果用的是混合模型，要改成mode = "satterthwaite"或"kenward-roger"。

多重比较校正别忘了。四个时间点各做一次两组比较，一共4次检验，用Bonferroni校正的话alpha除以4。emmeans里加adjust = "bonferroni"或者"holm"都行，我一般用Holm，因为它对检验之间的相关性不敏感，功效比Bonferroni略高。

如果交互不显著但时间主效应显著，那就退一步，只看所有患者合并后的时间点两两比较，用emmeans(mod_mlm, ~ time)配合pairs()即可。这时候不需要按组分拆，因为交互不显著意味着两组的趋势可以认为是一致的。

4. 常见问题与排查技巧实录

4.1 报错“number of items to replace is not a multiple of replacement length”

这个报错在构造idata或者模拟数据时特别常见，根源是向量长度不匹配。比如idata写了4行但cbind()里只绑了3列，或者模拟数据时rand_intercept的长度和id对不上。排查方法很简单：在报错的那一行之前，用length()和nrow()把涉及的每个对象的维度打印出来，逐个核对。我习惯在关键步骤后加一句stopifnot(length(x) == expected)做断言，早发现早解决，比事后debug省时间。

另一个隐蔽的变体是因子水平顺序不一致。idata里time的水平是t0, t4, t8, t12，但cbind()里列的顺序是t0, t8, t4, t12，R不会报错，但结果全错——因为它按位置匹配，不按名字匹配。这种错误最难查，因为代码能跑通，输出看起来也正常。我的做法是：cbind()里的列顺序和idata的行顺序，都用同一个变量控制，比如先定义time_levels <- c("t0","t4","t8","t12")，两边都引用它，从源头上杜绝错位。

4.2 缺失值导致样本量骤减怎么办

临床随访数据几乎不可能没有缺失，患者失访、漏检、样本溶血，各种原因都会造成某个时间点缺值。aov()和Anova()默认做完整案例分析，只要一个患者四个时间点里缺了任何一个，整行就被丢掉。20人的组缺了5个，样本量直接掉四分之一，检验效能大打折扣。

应对策略分两层。第一层是分析前的缺失机制判断：如果缺失完全随机（MCAR），完整案例分析虽然浪费信息但结论无偏；如果缺失与观测值有关（MAR），完整案例分析就有偏了，这时候应该上混合效应模型或者多重插补。第二层是实操层面的补救：用mice包做多重插补，插补后对每个插补数据集跑一遍重复测量方差分析，再用Rubin规则合并结果。这套流程代码量不小，但mice和mitml配合起来还算顺手。如果缺失比例很低（比如5%以内），完整案例分析的影响可以接受，但要在方法部分如实报告缺失情况，别藏着。

4.3 交互效应显著但事后比较全不显著

这种情况我遇到过好几次，第一次碰到时一头雾水——交互p值0.03，但四个时间点分别做两组比较，校正后一个都不显著。原因通常是交互效应体现在“趋势”上而非“某个具体时间点的截距差异”上。比如两组都是线性下降，但斜率不同，在中间某个时间点可能刚好交叉，单点比较自然不显著，但整体趋势差异是真实的。

解决办法有两个方向。一是做趋势分析，把时间因子拆成线性、二次、三次正交多项式对比，看交互效应主要落在哪个趋势成分上。emmeans里用contrast(emm, "poly")或者手动指定正交对比系数都行。二是做区域检验或者模型-based inference，比如拟合一个时间×组的线性回归，直接检验斜率差异。这两种方法都比逐点比较更贴合“趋势差异”的本质。我的经验是：如果研究假设本来就是“起效速度不同”，那从一开始就该规划趋势对比，而不是等到事后比较不显著了再补救。

4.4 问题速查表

4.5 几个容易被忽略的实操心得

第一个心得：先画图再跑模型。把每个患者的轨迹画出来（ggplot2里用geom_line(aes(group = id))），再叠加组均值轨迹，数据的形态一目了然。有没有离群个体、缺失集中在哪个时间点、两组轨迹是平行还是交叉，图上一看便知。很多统计上的困惑，画个图就解开了。我现在的习惯是，任何重复测量分析，图一定跑在模型前面。

第二个心得：样本量估算要在设计阶段做。重复测量设计的样本量估算和普通两组比较不一样，要考虑时间点之间的相关性（用预实验的相关系数估计）和球形假设的epsilon。pwr包和WebPower包都能做，但参数设定需要预实验数据支撑。事后补算意义不大，但方法部分最好报告一下估算依据。

第三个心得：报告结果时把效应量和置信区间带上。光报p值信息量不够，尤其是交互效应，报个偏eta方或者两组均值差的置信区间，读者才能判断临床意义。emmeans的confint()直接给置信区间，effectsize包可以算偏eta方，几行代码的事，别省。

第四个心得：时间因子的编码方式影响解读。如果时间点间隔不均匀（比如0、4、8、12周是均匀的，但0、2、6、12就不均匀），把时间当分类因子做方差分析会损失“间隔信息”。这时候可以考虑把时间当连续变量做趋势分析，或者用正交多项式对比。分类因子适合“每个时间点都关心”的场景，连续编码适合“关心变化速率”的场景，选哪个取决于研究问题。

5. 结果呈现与可视化建议

5.1 均值轨迹图：让交互效应一眼可见

重复测量分析的结果图，核心是展示“组均值随时间的变化轨迹”，同时把个体变异或置信区间叠上去。我常用的画法是：x轴时间点，y轴指标均值，两条线分别代表两组，每个点加误差棒（标准误或95%置信区间），底层用半透明细线画出每个个体的轨迹。这样既有群体趋势，又有个体异质性的信息。

library(ggplot2) dat_summary <- dat_long %>% group_by(group, time) %>% summarise(mean = mean(value), se = sd(value)/sqrt(n()), .groups = "drop") ggplot() + geom_line(data = dat_long, aes(x = time, y = value, group = id), color = "grey80", alpha = 0.4) + geom_line(data = dat_summary, aes(x = time, y = mean, group = group, color = group), size = 1.2) + geom_point(data = dat_summary, aes(x = time, y = mean, color = group), size = 3) + geom_errorbar(data = dat_summary, aes(x = time, ymin = mean - se, ymax = mean + se, color = group), width = 0.2) + theme_bw()

这张图的信息密度比单纯的两条均值线高得多。审稿人或者读者一眼就能看出：两组基线是否可比、轨迹是否交叉、个体变异有多大、哪个时间点开始分化。如果交互显著，图上两条线的分离趋势会很明显；如果交互不显著但时间主效应显著，两条线基本平行下降。

5.2 球形校正结果在表格里怎么摆

结果表格的呈现有个讲究：未校正和校正后的p值都要报，但主结论以校正后为准。我通常的表格结构是：效应项（组间、时间、交互）、F值、未校正自由度、未校正p值、epsilon值、校正方法、校正后p值。这样读者既能看原始结果，也能判断校正的影响程度。如果epsilon接近1，可以在表注里说明“球形假设满足，报告未校正结果”，避免表格过于臃肿。

效应量方面，偏eta方（partial eta squared）是重复测量方差分析里最常报的，计算方式是效应平方和除以（效应平方和+误差平方和）。组间、组内、交互各报一个。解释阈值参考：0.01小、0.06中、0.14大，但这是通用参考，具体领域可能有自己的惯例。

5.3 简单效应结果的文字化表达

事后比较的结果别只甩一张表，要用文字把关键对比说清楚。比如“在第8周，A组均值比B组低6.2个单位（95% CI: -9.8 至 -2.6，校正后p = 0.002），而基线时两组差异无统计学意义（p = 0.71）”。这种写法把时间点、组别、效应量、置信区间、校正后p值全带上了，信息完整且可读。如果对比很多，挑研究假设相关的几个重点写，其余的放补充材料。

提示：报告简单效应时，务必注明多重比较的校正方法。不同校正方法下p值差异可能很大，读者需要知道你的结论建立在哪种校正之上。

6. 写在最后的一点个人体会

这套流程我在好几个项目里反复用过，从最初的照葫芦画瓢到后来能快速判断“这个数据该用哪种方法”，中间踩的坑基本都写在上面了。如果只让我留一句话给刚接触重复测量分析的朋友，那就是：先把数据结构想清楚，再动手写代码。宽格式长格式怎么转、因子水平怎么设、缺失值怎么处理，这些“脏活”决定了后面模型跑出来可不可信。统计方法本身反而不复杂，aov()加Error()项或者Anova()加idata，核心代码就那么几行，难的是对设计逻辑的理解和对结果的诚实解读。

另外提一句，重复测量方差分析不是万能的。时间点很多、缺失较多、或者个体轨迹形态差异很大的时候，混合效应模型或者广义加性模型往往更合适。经典方法适合设计规整、时间点少、数据相对完整的情形，超出这个范围就该考虑升级工具了。工具没有高低之分，匹配研究问题才是关键。