重复测量数据在临床研究和实验科学里出现的频率极高——同一个受试者在不同时间点被反复观测比如术后第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周的38B组基线均值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指定个体内因子是timeidata的行数和顺序必须和cbind()里列的排列完全一致这是最容易出错的地方。输出里会包含time的球形检验结果Mauchlys test以及未校正和GG/HF校正后的F值和p值。组间效应看group那一行交互效应看group:time那一行。实测这份模拟数据交互项的GG校正p值应该在0.05以下说明两组的时间趋势确实不同。组间主效应可能不显著因为两组在第12周的均值差距被基线到中期的重叠稀释了——这正好说明为什么不能只看组间p值。3.3 球形检验与校正结果的取舍Anova()输出里会有一行Mauchlys 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, t12R不会报错但结果全错——因为它按位置匹配不按名字匹配。这种错误最难查因为代码能跑通输出看起来也正常。我的做法是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 问题速查表问题现象可能原因排查与解决Anova()报错维度不匹配idata行数与cbind()列数不一致用length()核对统一用变量控制顺序结果与预期完全相反因子水平顺序错位检查factor(levels)设定idata与列顺序对齐组内效应p值异常小球形假设违反未校正看Mauchly检验和epsilon选GG或HF校正事后比较全不显著交互体现在趋势而非单点改用正交多项式对比或斜率检验样本量比预期少很多缺失值完整案例删除评估缺失比例考虑多重插补或混合模型emmeans置信区间过窄自由度近似方式不对显式指定mode参数匹配模型类型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核心代码就那么几行难的是对设计逻辑的理解和对结果的诚实解读。另外提一句重复测量方差分析不是万能的。时间点很多、缺失较多、或者个体轨迹形态差异很大的时候混合效应模型或者广义加性模型往往更合适。经典方法适合设计规整、时间点少、数据相对完整的情形超出这个范围就该考虑升级工具了。工具没有高低之分匹配研究问题才是关键。