重复测量方差分析Repeated Measures ANOVA是我在带学生做项目时被问到最多的一类统计问题。之前有位同学拿着两轮的随访数据来找我说想比较同一批患者治疗前后三个时间点的血压变化问我用配对t检验行不行——我说第一轮、第二轮做配对自然没问题但三个时间点摆在一起再两两配对一类错误就失控了。这才是重复测量方差分析真正的用武之地处理同一批对象在多个时间点或多种条件下的连续测量结果。这篇文章我会尽量不用教科书语气把概念、原理、适用条件、Stata实现、结果解读串在一起讲清楚末尾还整理了一份零基础入门Stata的路径导读。不管你是临床医学、公共卫生、心理学还是经管类的研究生看完之后拿着自己的数据直接套用应该不会有障碍。1. 重复测量方差分析到底解决什么问题1.1 从实战案例理解“重复测量”四个字先说个最典型的场景你做了一项运动干预实验找24名高血压前期受试者分别在干预前、干预4周、干预8周、干预12周测量收缩压。现在问题来了——四个时间点、同一批人、同一个指标希望知道血压随时间有没有显著变化。如果换成完全随机设计把24人随机分成4组每组只在某一个时间点测一次血压那就是普通的单因素方差分析。但重复测量设计不是这样它是同一组人从第0周到第12周被追着测了4次。这组数据里除了“不同时间点造成差异”还藏着“同一个人的特质在重复出现”——张三就是偏高李四就是偏低这种来自个体本身的差异会在每次测量中都冒出来。重复测量方差分析的核心任务就是把这部分“老熟人带来的稳定差异”从误差中分离出去。说得直白点普通方差分析就是把总变异拆成组间和组内而重复测量方差分析把组内变异再拆成了“时间效应”和“个体随机差异”然后用更小的误差项去检验时间效应检验效能自然比不分离要高出不少。1.2 什么时候才该用重复测量方差分析判断标准实际上很朴素你只需要回答三个问题。第一因变量是否连续变量血压、得分、反应时、浓度这些都可以如果是二分类或者计数数据那就得走广义估计方程或者混合效应逻辑回归不在本方法射程内。第二是否来自同一批对象的多次测量这是重复测量方差分析与普通方差分析最根本的分水岭。干预前、中、后的纵向随访数据同一患者左右两只眼睛的术后视力数据同一批被试在三种实验条件下的正确率数据都属于这个家族。第三测量次数是否有限且相对较少一般建议测量时间点控制在2到5个左右比较稳妥。如果随访了10个时间点甚至更多重复测量方差分析的球形假设会非常脆弱直接上混合效应模型反而更耐造。还有一个容易被忽视的细节重复测量方差分析要求测量时间点或条件对所有对象一致。也就是说所有人都是第0、4、8、12周测张三不能拖到第13周李四不能缺第8周。如果随访时间漂移严重或者有人中途脱落数据变成不平衡状态这时更靠谱的选择是混合效应模型而不是硬凹重复测量方差分析。2. 原理拆解方差从哪里来F值是怎么算的2.1 变异被拆成了几块想真正看明白Stata输出的那张方差分析表你脑子里得有一张变异的分解地图。对于最简单的一个组内因子时间、一个被试因子编号的设计总变异可以拆成三部分被试间变异不同人的平均水平差异、时间变异不同时间点的整体变化、误差变异就是去掉人和时间效应后剩下的随机波动。这里有个容易绕晕的概念被试间变异虽然在重复测量方差分析里算一个方差来源但我们通常不关心它对应的F检验是否显著因为“人本身有差异”是既定事实。我们真正盯着的是时间效应这个F值——它的分子是时间变异分母是误差变异。相比普通方差分析把被试间变异混进误差这里的误差项被单独拎了出来数值明显更小这也就是为什么重复测量设计更容易检测到微小的时间效应。如果实验设计里还有分组变量比如实验组与对照组变异分解就会再多一个被试间的分组效应以及“时间×分组”的交互效应。交互效应稍微绕一点它回答的是“两组人的血压随时间变化的轨迹到底一不一样”这个问题在干预研究里往往比时间主效应还要关键。2.2 球形假设为什么这一条是灵魂重复测量方差分析能顺利运行有个前提条件叫球形性sphericity也叫复合对称性的推广版本。你可以把它理解成“不同时间点之间的差值的方差应该是差不多的”。换句话说第0周到第4周的血压差值、第4周到第8周的差值、第8周到第12周的差值这些差值的离散程度不能差太远。为什么这个条件这么重要因为我们前面说过误差项是把所有时间点的数据合在一起凑出来的。如果时间点之间的差值变异相差悬殊那这个合并的误差项就无法真实代表每个时间点成对比较的误差算出来的F值就会偏大一类错误率就失控了。球形假设的判断方法Stata会在anova命令的重复测量输出里自动给出Mauchlys W检验以及Greenhouse-Geisser、Huynh-Feldt两组校正系数。W检验的p值如果小于0.05说明球形假设不成立这时候你必须看校正后的p值不能继续看未校正的结果。这是初学者翻车率最高的一个环节后面讲结果解读时单开一节细说。2.3 检验效能与样本量直觉很多人会问重复测量方差分析最少要多少样本量。虽然没有绝对铁律但有个经验公式可以感受一下因为每个对象提供了多次测量同等检验效能下重复测量设计的样本量往往只需要完全随机设计的六到七成。但还是得提个醒样本量的计算并没有想象中那么简单。你还需要知道组内相关系数或者预期的效应量这时可以用G*Power之类的软件去估算。如果只给一个粗糙参考临床研究中一个两个水平的时间因子、想检测中等效应每组30例左右比较常见。时间点越多单组样本量可以稍微降一降但样本量过低时球形检验的把握度也会很低结果的可信度会打折扣。3. 实操从零跑通一个完整案例3.1 数据长表还是宽表reshape必须会绝大多数新手栽在数据分析的第一关不是不会跑命令而是数据摆得不对。重复测量数据在Stata里有两种常见形态宽表比较好理解一个人占一行第0周、第4周、第8周、第12周的血压分别占四列。长表则是一个人占四行用time变量区分是第几个时间点再加一列血压值。anova repeated这类方差分析命令通常直接接受宽表即可。但后续画趋势图、做缺失值处理、运行混合效应模型时几乎全是围绕长表工作。所以我会建议你务必学会reshape命令。* 假设当前是宽表: id sbp0 sbp4 sbp8 sbp12 reshape long sbp, i(id) j(time) * 转换后变成: id time sbp这里i()是每个人唯一的编号j()是时间点的变量名。转换完成后会多出一个time变量取值0、4、8、12。想从长表再转回宽表用reshape wide即可。这条命令几乎可以算Stata数据管理中最重要的命令之一宁可多花十分钟练熟也不要到后面分析时手忙脚乱。3.2 anova命令配合repeated()选项在Stata里跑重复测量方差分析经典命令其实很简洁。假设因变量是sbp被试编号是id时间因子time有四个水平。anova sbp id time, repeated(time)这一行命令背后做的事情是将id视为被试间因子time作为被试内因子采用重复测量的误差估计方式来检验时间主效应。运行结束后Stata会输出一份方差分析表并在表下方额外给出球形检验以及校正后的显著性结果。如果实验里还有组别group想检验组别主效应、时间主效应和交互效应命令要略微调整anova sbp id group time group*time, repeated(time)这里注意写法被试间部分要包含id和group但id通常是嵌套在group里的Stata会自动处理。实际项目中我的习惯是先跑一遍全模型看交互项是否显著。交互效应不显著时再考虑精简模型或只解读主效应。交互效应显著时我会直接用下一节讲的contrast或margins去做简单效应分析而不是停留在主效应上。3.3 更稳妥的进阶路子mixed命令说句掏心窝的话我在自己日常研究里现在已经很少用anova repeated了。数据平衡时两者结果接近但从零基础入门角度来看anova repeated的计算逻辑更直观适合学习原理。读了研究生以后或者随访时间点较多时我更推荐用mixed命令拟合线性混合模型因为它的灵活性和稳健性确实更好。mixed sbp time || id:, repeated(time) covariance(unstructured)这行命令将时间time作为固定效应使用非结构化协方差矩阵来刻画重复测量之间的相关性。相比方差分析它不要求球形假设也能自然处理个别缺失数据。不过做混合模型需要理解固定效应和随机效应的区别刚开始学的时候不必强求先把自己课题数据用anova repeated跑通再过渡到mixed会顺畅不少。如果你要在论文里详细描述方法学部分建议把分析策略交代清楚“采用重复测量方差分析检验不同时间点之间的差异球形假设不满足时采用Greenhouse-Geisser校正事后两两比较采用Bonferroni校正。”这一句话虽然简短审稿人看了就知道你基本功扎实。3.4 事后两两比较怎么做仅仅知道时间主效应显著还不够你得知道到底是哪个时间点和哪个时间点不一样。以干预前为基准的情况可以直接用contrast命令contrast {time 1 -1 0 0} {time 1 0 -1 0} {time 1 0 0 -1}, nowald这里每组花括号代表一次比较1 -1 0 0表示第一个时间点与第二个时间点相减。如果觉得这种写法太费脑子更简单的方式是切换到混合效应模型做边际均值比较mixed sbp time || id: margins time, pwcompare(bonferroni) effectspwcompare会给出所有时间点两两比较的结果加上bonferroni选项后多重比较校正直接完成。注意混合模型的随机斜率设定可以更加灵活mixed sbp time || id: time, covariance(unstructured)这样每个人不仅有自己的基准水平还会允许随时间的变化速度不同。如果这个模型的随机斜率方差显著不为零说明个体间的变化轨迹差异确实存在模型的拟合也更加细腻。4. 结果解读Stata输出表里到底写了什么4.1 看懂表里的主效应和交互效应很多同学跑完命令看到一大片数字就慌了。实际上Stata对重复测量方差分析结果的输出相当友好你只需要按顺序抓住三个区域。第一个区域是经典的方差分析表里面有Source、SS、df、MS、F和ProbF。你需要看的主要是time这一行的F值和p值。在完全被试内设计里被试间行和时间行的试验误差是分开的这也是重复测量方差分析与普通方差分析读表时的最大区别。Stata还有一行显示Error项时间效应的误差是误差项被试间效应的误差是id的误差项千万别看串了。第二个区域是紧跟在表下方的球形检验信息包括Mauchlys W、近似卡方值和p值以及Greenhouse-Geisser epsilon、Huynh-Feldt epsilon等。如果W检验p值大于0.05说明球形假设成立直接读上方未校正的time F值和p值就行。如果W检验p值很小表示球形假设被拒绝这时你需要往下读带Greenhouse-Geisser校正的那行p值。第三个区域是校正后的结果表。Stata会明确标注“Greenhouse-Geisser epsilon”和“Huynh-Feldt epsilon”同时提供两组校正后的p值。一般习惯以Greenhouse-Geisser为准因为它相对保守。Huynh-Feldt更宽松一点当epsilon接近1时两者差别很小。epsilon最低值是1/(k-1)k是时间点数当k4时最低epsilon是0.333如果你的epsilon小于0.5那说明球形性偏差已经比较严重了。4.2 效应量同样重要只汇报p值不汇报效应量是很多论文被审稿人打回的重灾区。重复测量方差分析的常用效应量指标是偏eta方写成部分η²。偏eta方的计算逻辑不复杂该效应平方和除以该效应平方和加对应误差平方和。Stata本身不会直接给出偏eta方但你可以根据方差分析表手动计算或者用estat esize命令尝试获取。手动计算也很简单display 120 / (120 80)比如time效应的平方和是120误差平方和是80那偏eta方就是0.6。这个数怎么判断大小呢大致上0.01左右算小效应0.06左右算中等效应0.14以上算大效应。注意这个参考值来自心理学行为研究的经验习惯不同学科门类略有出入但作为通用判断已经完全够用。如果你的导师要求比较规范论文里通常会写成这样“时间主效应显著F(2.34, 53.82) 6.03, p 0.003, 偏η² 0.208。”括号里的自由度不是整数就是因为使用了Greenhouse-Geisser校正后自由度被调整了这也是懂行读者的信号。5. 常见问题与避坑清单5.1 球形假设不满足怎么办这个问题几乎每个做纵向数据的人都躲不开。Mauchlys W检验p值小于0.05时最直接的应对就是读取Greenhouse-Geisser校正后的p值。如果报告了校正后结果仍然显著那结论基本站得住。如果校正后不显著说明原本的显著性可能是被夸大的这时你应该老实接受校正结果。さらに如果epsilon特别低、时间点数量又多我更建议避开这种纠结直接转向线性混合模型。混合模型通过指定协方差结构如非结构化、复合对称、一阶自回归来处理多次测量之间的相关性问题既不依赖球形假设又能保留所有观测。Stata里配合mixed命令加repeated()选项就能做到输出里还会给出协方差参数估计比方差分析信息量大得多。5.2 缺失值、异常值、数据类型要早检查重复测量数据最怕的就是时间点缺失。如果你用的是anova repeatedStata一般会把带有缺失值的观测整行剔除这意味着只要有缺失这个人的全部数据点都没法进入分析十分可惜。所以拿到数据第一件事就是检查缺失模式misstable summarize misstable patterns这两条命令能快速告诉你缺失的分布。如果缺失比例不高、模式为随机缺失可以先用多重插补再跑分析。如果缺失集中在某些时间点过去的经验是尽量少用整行删除优先考虑混合效应模型它能比较好地利用所有可用数据。还有一种坑是变量类型问题。有时候Excel里患者编号不小心被识别成数值变量跑到reshape时反复报错有时候时间变量被存成字符串命令报错说type mismatch。我习惯在建项目之初就把id变量转成字符串或数值保持统一的命名规范省市避免后面浪费大量时间排查。5.3 实测下来的几条操作心得第一条心得不要一上来就点菜单。Stata的菜单操作虽然看起来省事但可重复性太差换一台电脑、换一批数据你就不知道当时点了什么。所有操作尽量用do文件记录下来文件名、路径、每一步都标注清楚。第二条心得变量命名别用中文也尽量别用太长的大小写混合名。Stata对大小写敏感血压变量写成sbp就是sbp写成SBP就会报变量不存在。统一用首字母缩写加数字比如sbp0、sbp4、sbp8、sbp12一眼就能看懂。第三条心得如果样本量非常小比如总共10个人重复测量方差分析的检验效能很弱这时球形检验也容易不显著但别误以为条件满足就万事大吉。小样本下更要优先做描述统计和个体轨迹图数据质量问题往往在图上暴露得最快。6. 零基础小白入门Stata统计学/计量经济学学习路径6.1 阶段一统计基础必须补的几块砖零基础直接上手Stata最大的问题不是命令不熟而是统计概念缺失导致不知道自己跑出来的是什么。我建议先花两周夯实三个基础知识点描述统计与图表、正态分布与假设检验逻辑、t检验与方差分析的基本思想。这个阶段不必钻研数学推导只要明白p值是什么、第一类错误和第二类错误是怎么回事、方差分析为什么看F值就行。推荐《统计白话》这类用大白话讲统计的入门书一天看一两章配合网上的免费视频比啃大部头教材效率高得多。6.2 阶段二Stata软件与数据管理基本功下载安装Stata之后第一件事不是跑回归而是学会数据管理。数据管理占日常数据处理工作量的七成以上。你需要掌握的技能清单包括数据导入与清洗、变量生成与标签、子集选取和条件筛选、merge和append合并、reshape的长宽转换。以reshape为例这条命令对面板数据和重复测量数据来说如同命门。我的经验是找一份别人分享的练习数据把宽表和长表来回转换三次错了就重来半小时基本就能熟练。6.3 阶段三方差分析家族的系统学习方差分析不是孤立的一个命令而是一个家族。按照复杂程度排序我建议的学习路线是单因素方差分析、双因素方差分析、重复测量方差分析、协方差分析。每一类都要结合自己的学科数据做一次完整的分析流程从数据准备到结果解读再到图表展示。学完重复测量方差分析之后如果你有余力可以进一步了解一下其在meta分析中的应用场景比如剂量-反应meta分析中的数据整理常常涉及多次测量点这时Stata的网状meta分析相关的命令也会有帮助。不过那是专题性知识入门阶段不用强求先把自己课题的核心数据结构弄清楚更重要。6.4 阶段四回归分析与计量经济学入门方差分析本质上也是一种线性模型理解了方差分析之后进入回归分析会更顺滑。这个阶段重点是理解最小二乘法的直觉、回归系数的解释、虚拟变量的含义以及如何处理异方差和多重共线性问题。Stata里regress命令是永远的起点配合estat vif做共线性诊断是必修操作。经管类同学再往后走还会遇到固定效应和随机效应模型、面板数据分析这些都是重复测量思想的扩展版。你会发现当初在重复测量方差分析里学到的“个体稳定差异”概念到了面板数据里就是个体固定效应底层逻辑完全相通。6.5 阶段五面向实战的输出能力最后一个阶段其实是最容易被忽视的。你会发现真正拉开差距的不是会不会跑某个模型而是能不能把统计结果写成规范的论文段落。我建议大家准备一个自己的检查清单研究设计是什么类型、数据是否符合模型前提、有没有做多重比较校正、有没有报告效应量、图表是否清晰完整。另一个实用建议是把常用命令整理成自己的工具包。比如把重复测量方差分析、描述统计、数据清洗的do文件模板存好以后拿到新数据只改变量名就能直接跑。Stata里有一些第三方命令包比如ftool等工具可以根据需要随时安装使用但核心命令一定要做到不看帮助文档也能正确背诵。我在实际带教中总爱说一句话重复测量方差分析其实是你理解纵向数据的一把钥匙。你会发现钥匙本身能打开的门是有限的但它教会你的思维方式——把个体差异从随机误差中拆出来——几乎贯穿了你后面所有的计量和统计模型学习。所以哪怕你现在只为一个课题而学也值得静下心来把这套逻辑吃透这很可能是你往后数据分析路上最划得来的一笔时间投资。