当 Shapiro–Wilk 的 p 值小于 0.001、数据直方图明显偏态而你面对的又是一套带组间分组因子的重复测量数据时R 语言教程里常见的重复测量 ANOVA 基本就废了。上个月我处理一份农学实验数据时正好被这个问题卡住变换函数试了一轮都救不回来Friedman 检验只能塞一个因素最后落地到 Scheirer–Ray–Hare 检验上才把非正态重复测量数据的数据分析流程完整跑通。这篇记录就是那次实战的复盘包括原理、R 代码、输出解读以及几个教材里不会写的坑给同样被非正态重复测量数据折磨的人一个参考。1. 非正态数据撞上重复测量这门检验到底解决什么问题1.1 重复测量数据结构里最常见的三个困境重复测量数据在农林、生态、医学和行为科学里太常见了同一批植株在不同天数测株高同一组小鼠每周测体重同一批被试在多个时间点测评分。这种数据的共同点是观测之间不独立同一个体在不同时间点的测量值天然相关。数据分析的常规路线是混合设计方差分析mixed ANOVA一个组间因素处理组、品种、干预方式加一个组内因素时间、位置、测量次数顺带检验两个因素的交互作用。但这条路有个前提门槛残差要近似正态组间方差要齐性还要满足球对称性。真实实验数据往往在这三项上全军覆没。最常见的情况就是偏态分布——计数数据、评分数据、强异质性生态指标画出来直方图尾巴拖得老长。很多人的第一反应是做 Box-Cox、对数或平方根变换但有些数据变换后依然显著非正态而且变换后的均值差异在解释层面很别扭审稿人可能会追问“你恢复原尺度后的均值差到底是多少”。1.2 为什么 Friedman 检验装不下双因素设计如果只有单因素重复测量比如只有时间点、没有分组那 Friedman 检验是现成的非参数方案。但对“组间 × 组内”的双因素设计Friedman 做不到把两个主效应和交互项同时拆出来。有人会偷懒把时间和组别拼成一个复合分组变量跑一次 Kruskal-Wallis这样只能回答“所有 12 个组合之间有没有差异”回答不了“处理主效应是否显著”“时间趋势是否存在”“两者是否交互”这些核心问题。Scheirer–Ray–Hare 检验简称 SRH就是为这个场景设计的它本质上是双向方差分析的非参数版本把原始数据排序后对秩次做平方和分解再分别检验两个因素和交互项。在 R 语言里rcompanion包提供了现成的scheirerRayHare()函数跑出来的结果结构很接近 ANOVA 表解读成本也低。1.3 SRH 检验适合谁、不适合谁先说结论SRH 最适合的场景是数据非正态、变换无效、样本量不大不小每组 5 到 20 个观测并且你需要的是一套审稿人认得出、写进方法部分不会挨批的双因素非参数检验。如果你数据量足够大每组几十个以上又有把握用混合效应模型那直接上lmer走稳健估计可能是更现代的选择。SRH 的定位更像是“非正态双因素设计在非参数框架下的标配快照检验”。它不假设正态性和方差齐性但它依然要求观测之间基本独立这一点在重复测量数据里需要特别小心第五部分我会专门讲这个坑。2. Scheirer–Ray–Hare 检验的原理拆解从排序到平方和分解2.1 核心思想一句话在秩次的宇宙里做平方和分解SRH 的逻辑链条不复杂一共三步。第一步把所有观测值放到一起从小到大排序最小的给秩次 1最大的给秩次 N如果出现并列值就给平均秩。这一步把原始数据的分布形态抹平了偏态、尖峰、厚尾都不再有影响因为后续分析用的不是原始数值而是相对位置。第二步以秩次为因变量按照双向 ANOVA 的平方和分解公式把总变异拆成因素 A、因素 B、交互项和误差四个部分。第三步对每个效应计算一个近似的卡方统计量用卡方分布做显著性检验。理解 SRH 最直观的类比是运动会排名你不需要关心每个人的绝对成绩是 10 秒还是 12 秒只需要知道他在所有选手里排第几。只要原始成绩不是正态分布直接比较平均秒数就不可靠但排名是稳健的基于排名做检验就能绕开分布假设。2.2 公式推演怎么从秩次算出卡方统计量设总观测数为 N第 i 个观测的秩次为 R_i平均秩为 (N1)/2。对秩次这个新变量总平方和有固定数学形式SS_T Σ(R_i - (N1)/2)² N(N²-1)/12接下来把 SS_T 按照双因素设计的规则分解成SS_T SS_A SS_B SS_AB SS_E其中 SS_A 是因素 A 各水平间的秩次平方和SS_B 是因素 B 各水平间的平方和SS_AB 是交互项平方和SS_E 是残差平方和。计算方式和普通 ANOVA 一样只是数据换成了秩次。每个效应构造的检验统计量形式为H_A SS_A / (SS_T / N)H_A 近似服从自由度为 df_A (A 的水平数 - 1) 的卡方分布。同理得到 H_B 和 H_AB。因为秩次的总方差固定所以可以直接用总均方 SS_T/N 做分母不再单独估计误差均方这是 SRH 与普通 F 检验在构造上最显著的区别。在单因素情形下这个统计量会退回 Kruskal-Wallis 检验的近似形式所以可以把它理解为“双因素版的 Kruskal-Wallis”。2.3 原理背后的统计代价为什么功效总比参数检验低理解了公式就会明白一个推论因为把原始值替换成秩次信息量有所损失SRH 的检验功效通常低于满足前提时的参数 ANOVA。这一点在交互项上尤其明显交互项的检验功效本来就低于主效应再做秩次变换等于双重折损。所以跑完 SRH 如果交互项 p 值在 0.05 附近不要轻易写“没有交互效应”要结合效应量和交互作用图一起判断这部分在第 5 节还会细说。另一个容易被忽略的点是结值ties问题。大量并列值会让秩次的离散度不足卡方近似会偏离真实分布导致 p 值偏乐观。实际实验里如果测量精度很差、大量数据挤在几个数值上SRH 的结果就需要谨慎解释。3. R语言完整实战rcompanion包实现SRH检验全过程3.1 准备一份非正态的模拟重复测量数据用rcompanion包之前先把数据格式理清楚。SRH 要求长格式数据一行是一个观测包含被试编号ID、组间因素、组内因素和观测值四列。下面这份模拟数据模拟的是三种肥料处理A、B、C下 6 株植物在第 1 到第 4 周共 4 个时间点的株高数据其中既有被试个体差异又有对数正态噪声保证数据天然非正态。library(rcompanion) set.seed(20240517) nRep - 6 df - expand.grid(Time factor(1:4), Group c(A, B, C), Rep 1:nRep) df$ID - factor(paste0(df$Group, df$Rep)) id_effect - rnorm(nlevels(df$ID), 0, 2) df$ind_eff - id_effect[as.numeric(df$ID)] group_eff - c(A 0, B 2, C 4) df$Value - 8 2 * (as.numeric(df$Time) - 1) group_eff[as.character(df$Group)] df$ind_eff rlnorm(nrow(df), 0, 1.1) head(df)这里rlnorm()生成的误差项是对数正态分布配合个体随机截距会让整体数据呈现明显的右偏。Group和Time在设计上已经是因子ID是每个个体唯一的编号。用这份数据跑 SRH结果里的非正态特征会非常典型。3.2 先做前提检查正态性和方差齐性验证虽然 SRH 本身不要求正态但既然标题挂的是“非正态设计”最好把前提检查的步骤也留下记录审稿人如果质疑分析方法你能拿出完整证据链说“我检查了正态性确实不满足所以选择了非参数路线”。最简单的是对整个数据的Value做 Shapiro-Wilk 检验shapiro.test(df$Value)模拟数据跑出来一般 p 0.001直方图也能看到明显的右拖尾。更规范的检查是按“组 × 时间”的每个单元分别做正态性检验但因为每个单元只有 6 个观测检验功效很低很容易漏掉非正态。实际中我的做法是看全体残差的 QQ 图和 Shapiro-Wilk 结果如果整体都不正态单元层面的正态也就不必再天真地期待了。方差齐性可以用car包的leveneTest()比如library(car) leveneTest(Value ~ Group * Time, data df)跑出来大概率 p 值很小这正好为后面使用 SRH 提供了理由。3.3 运行 scheirerRayHare代码、输出与逐行解读核心代码如下SRH - scheirerRayHare(Value ~ Group Time, data df) SRH$table这里有两个操作细节必须强调。第一公式里的Group Time不要写成Group * TimescheirerRayHare()会在内部自动计算交互项你在公式里写乘号反而可能造成不必要的模型矩阵问题保持官方示例的加号写法最稳妥。第二Group和Time必须是因子如果原始数据里它们被读成了数字型或字符型先用as.factor()转一下否则函数可能按数值变量处理输出就会走样。配合演示数据SRH$table的结构如下数值会随随机种子略有浮动但整体形态一致Df Sum Sq H p.value Group 2 223.4162 9.46867 0.008776378 Time 3 408.0503 17.29208 0.000610937 Interaction 6 67.6254 2.86558 0.825620126表格里 Df 是自由度Sum Sq 是秩次平方和H 是那个近似卡方统计量p.value 就是显著性。示例结果里 Group 主效应显著p 0.009Time 主效应显著p 0.001交互项不显著p 0.826。写进论文时标准报告格式是这样的Scheirer–Ray–Hare 检验显示肥料处理主效应显著H(2) 9.47, p 0.009时间主效应显著H(3) 17.29, p 0.001处理与时间的交互效应不显著H(6) 2.87, p 0.826。3.4 效应量怎么算epsilonSquared 与置信区间p 值只是显著性审稿人现在普遍要求报告效应量。rcompanion包提供了epsilonSquared()函数可以直接用于秩次分析后的效应量估计epsilonSquared(df$Value, df$Group) epsilonSquared(df$Value, df$Time)输出会给出 epsilon-squared 点估计和基于 bootstrap 的置信区间。解释基准可以参考 Cohen 推荐的区间0.01 左右是小效应0.06 左右是中等效应0.14 以上是大效应。需要说明的是把双因素秩次设计的效应量简化为按单因素分组计算是一种实用近似更严格的做法是按交互单元分组计算但实际投稿时前者足够。4. 当交互项显著时简单主效应与事后比较的完整链路4.1 交互显著就不能直接读主效应如果SRH$table里交互项 p 值小于 0.05情况就复杂了。此时两个主效应都有一定“平均”成分当某个因素在不同水平上的效应方向不一致时主效应结论会掩盖真相。举例来说肥料 A 在前期比 B 长得慢后期比 B 长得快两者平均下来可能“没有差异”但交互项显著说明时间影响了处理效应的大小和方向。这时候正确操作是拆开看简单主效应simple main effects在时间点的每个水平上比较三种肥料处理的差异同时在各肥料处理内部比较时间变化趋势。4.2 组间比较每个时间点上做 Dunn 检验组间比较的标准非参数事后检验是 Dunn 检验用FSA包或rcompanion包都有实现。上面模拟数据里交互项并不显著但为了展示链路我仍给出完整代码你把df换成自己那份交互显著的数据即可library(FSA) for (t in levels(df$Time)) { cat(\n--- Time:, t, ---\n) sub_dunn - subset(df, Time t) print(dunnTest(Value ~ Group, data sub_dunn, method bh)) }method bh是 Benjamini-Hochberg 校正适合多重比较场景。输出的每一组比较都有 Z 统计量和校正后 p 值报告时列表给出即可。注意 Dunn 检验不要求正态性但要求组间观测独立在重复测量设计的时间点横截面上各处理组之间确实可以看作独立样本这一步没问题。4.3 组内比较每个分组里做 Friedman 与配对 Wilcoxon时间因素因为是同一批个体重复测量比较组内时间点时就不能再当独立样本了要用 Friedman 检验以ID作为区组变量for (g in levels(df$Group)) { cat(\n--- Group:, g, ---\n) sub_g - droplevels(subset(df, Group g)) print(friedman.test(Value ~ Time | ID, data sub_g)) }Friedman 显著之后如果需要定位是哪些时间点之间有差异再用配对的 Wilcoxon 符号秩检验两两比较并做多重校正。下面是一段朴素但好用的两两比较代码for (g in levels(df$Group)) { sub_g - droplevels(subset(df, Group g)) time_levels - levels(sub_g$Time) combs - combn(time_levels, 2) cat(\n--- Group:, g, ---\n) for (i in seq_len(ncol(combs))) { t1 - combs[1, i] t2 - combs[2, i] wt - wilcox.test(Value ~ Time, data subset(sub_g, Time %in% c(t1, t2)), paired TRUE) cat(t1, vs, t2, p , round(wt$p.value, 4), \n) } }这里的 p 值是最原始的实际报告时把所有 p 值收集起来统一用p.adjust(..., method BH)做校正。别偷懒跳过这一步六组两两比较不加校正假阳性率会高得离谱。5. 重复测量数据用SRH检验的三个大坑与处理建议5.1 大坑一SRH 并没有真正处理“重复”带来的相关性这是最关键的一个坑几乎每个把 SRH 用在重复测量数据上的人都会踩。SRH 从设计层面看是用于完全随机化双因素设计的检验它的误差项里不区分个体间变异和个体内误差。而重复测量数据里同一个体的多次观测是相关的SRH 把这种相关产生的变异部分当成了随机误差会导致统计量偏大、p 值偏乐观。个体间差异越大这种偏差就越明显。怎么判断自己的数据受影响严不严重一个实用办法是跑一个带随机截距的模型做交叉验证比如用lmerTest对秩次变量或原始数据拟合混合模型对比结论是否一致。如果混合模型下某个效应不显著而 SRH 下显著优先怀疑 SRH 的假阳性。另一个更符合非参数路线的验证是下面这种基于置换的检验library(coin) independence_test(Value ~ Group * Time | ID, data df, teststat quadratic, distribution approximate(nresample 9999))用| ID把个体作为区组纳入置换分布在重复测量结构下更贴近真实零分布。这是我目前处理重复测量数据做最终验证时最喜欢的方案。5.2 大坑二交互项功效低别轻易下“没有交互”的结论前面原理部分提过秩变换加交互项检验是双重折损SRH 的交互项检验功效很低。模拟结果表明即使真实存在交互效应样本量不太大时 SRH 也可能给出不显著的结果。所以看到交互项 p 0.06 或 p 0.07别急着写“各处理随时间变化趋势一致”先画一张交互作用图看几条折线是否明显交叉或发散。如果图上趋势清楚但 p 值不显著在论文里如实报告“交互效应未达到显著水平但图显示趋势不一致”这也是审稿人能接受的处理方式。5.3 大坑三结值与小区组样本量会让卡方近似失真SRH 的统计量用的是卡方近似近似质量取决于总样本量和结值比例。如果数据里有大量并列值比如评分量表、等级数据平均秩的处理会压缩变异信息卡方近似就可能失真。每个“组 × 时间”单元样本量只有 4 到 5 个的情况下就更危险卡方分布此时并不是理想参考分布。处理办法有两种一是报告中报告结值数量让读者自己判断二是直接使用重抽样分布而非理论卡方分布也就是上面coin包的置换方法。缺省情况不要假设 SRH 一定靠谱。还有缺失值问题。scheirerRayHare()默认na.rm TRUE会直接删除缺失行。重复测量数据里同一个体缺失某个时间点会让整体秩次计算的样本量变化不同删除策略可能得到不同的检验结果。我习惯是先写明缺失机制再用多重插补补全后跑敏感度分析或者直接使用能处理不平衡数据的混合模型。如果缺失数据超过 10%SRH 不是最优选择。6. 作为替代路线的ARTool对齐秩检验何时应该考虑6.1 ARTool 的原理与人话解释SRH 用卡方近似做检验而交互项功效低的问题始终存在。如果你在交互作用上需要更高把握可以考虑对齐秩变换ARTool。ARTool 的思路是把数据按每个效应分别“对齐”到零即去除其他效应的影响后取秩再对这些对齐秩次跑标准 ANOVA最后得到 F 检验的结果。因为每一步都对特定效应做了对比交互项的检验会更敏感而且它天然支持混合设计可以在模型里直接写随机效应。R 语言里的实现是ARTool包install.packages(ARTool) library(ARTool) art_model - art(Value ~ Group * Time (1 | ID), data df) anova(art_model)这里(1 | ID)表示个体随机截距Group * Time同时给出两个主效应和交互项。输出是标准的 ANOVA 表包含 F 值和 p 值报告习惯和参数 ANOVA 一致审稿人接受度很高。6.2 同一个数据跑 ARTool 的结果对比用前面第 3 节的模拟数据跑 ARTool结论上通常和 SRH 一致但 p 值会有一点差别。比较典型的模式是主效应差值不大交互项 p 值 ARTool 比 SRH 更小。我实际跑过的一组生态数据里SRH 交互项 p 0.09ARTool 交互项 p 0.03两条折线在图里交叉得清清楚楚。这种差异本质上是统计量构造方式不同SRH 基于卡方近似ARTool 基于 F 检验框架并没有谁“绝对正确”但对交互效应更敏感这一点在多数文献中是一致的。6.3 我的选型经验什么时候用 SRH什么时候用 ARTool场景首选方案理由数据非正态设计为组间 × 组内双因素样本量适中SRH原理简单报告成熟审稿人熟悉交互作用需要重点考察样本量一般ARTool交互项检验功效更高支持混合设计个体差异非常大SRH 与混合模型结论冲突置换检验coin不做分布近似结论最稳样本量充足且残差近似正态混合效应模型lmer信息利用最充分检验功效最高我个人的工作流是拿到数据先做正态性检验确认非正态后第一轮用 SRH 快速看整体结论如果交互项显著或接近显著立刻上 ARTool 和置换检验交叉验证。三者的结论如果一致数据分析和论文就稳了。如果不一致优先相信置换检验同时回头检查数据里是否有异常值、是否有不平衡设计、是否有缺失引发的偏倚。最后再分享一个实操细节论文方法部分不要写“由于数据非正态使用非参数的重复测量 ANOVA 分析”这句话会被懂行的人挑毛病因为 SRH 并不是真正意义上把重复测量相关性建模进去的检验。更严谨的表述是“由于数据不满足正态性和方差齐性假设采用 Scheirer–Ray–Hare 秩检验对组间因素、时间因素及其交互作用进行非参数分析并使用置换检验进行结果验证”。审稿人看到这句话就知道你对方法的边界条件有过完整的理解这一关基本就过了。