去年帮朋友处理一份植物多样性实验数据时，我遇到了一个挺典型的问题：处理组和采样时间两个因素都有，但数据是α多样性指数，Shapiro-Wilk检验的p值小到R直接给出警告，方差不齐，还带着明显的右偏尾巴。常规思路是双因素方差分析，可正态假设这一关就过不去。后来我用R语言里的Scheirer–Ray–Hare检验解决了问题，整个过程踩了不少坑，也把这个方法的适用边界翻来覆去想了好几遍。这篇文章就记录我完整走通的一次实操流程，包括原理、代码、结果解读、事后比较，以及那些文档里不写但实际会遇到的问题。适合做生态、农学、心理、医学研究，手头正好有“非正态 + 多因素 + 重复采样”这类数据的朋友参考。

1. 什么时候会碰到“非正态的重复测量数据”

1.1 一个典型的“处理 × 时间”实验场景

以生态学里常见的α多样性调查为例。你设计了两种土壤改良方式，分别记为Control和Treatment，在三个采样时间点T1、T2、T3分别采集样方数据，计算每个样方的Shannon多样性指数。这个数据结构里有两个分类变量：处理组（2个水平）和时间点（3个水平），目标是检验处理、时间以及两者的交互作用对多样性指数是否有显著影响。

这类数据有个共同点：α多样性指数、丰度数据、计数类指标，大多数情况下都不服从正态分布。它们往往右偏，偶尔还会冒出几个离群值。如果你做过几个生态项目就会发现，对这类数据直接做双因素方差分析，审稿人大概率会问一句：正态性检验做了吗？方差齐性满足吗？这不是形式主义，而是因为参数检验的结果在这些条件下确实不可靠。

这里的“重复测量”要稍微解释一下。实验设计上存在两种常见情况：一种是每个时间点重新取不同样方，不同时间的样本本质上是独立的；另一种是固定样方长期追踪，同一个体被重复测量了多次。本文中的案例默认是前一种破坏性采样模式，即时间点之间是独立样本。至于后一种真正意义上的重复测量，能不能用Scheirer–Ray–Hare检验，我会在常见问题部分专门讲清楚。

1.2 为什么不能直接套用双因素方差分析

很多人对双因素方差分析有个误解，觉得它“挺稳健的，稍微偏离正态应该没关系”。这句话只对了一半。F检验在样本量较大时确实对轻度的非正态不敏感，但对偏态分布、重尾分布、异方差组合在一起的情况，稳定性会明显下降。研究里经常提到的一个风险是I型错误膨胀：本来处理没有效应，但因为数据偏态，p值被压得很低，结果得到虚假的显著结论。这在样本量不均衡时更明显。

举一个实际感受过的例子：一次我在一个模拟实验里生成了两组均值相同但分布不对称的数据，用双因素ANOVA检验，p值反复出现在0.03到0.05之间，看起来“刚刚显著”。但把数据换成秩再做同样的分析，p值立刻回到0.3以上。数据本身没有差异，是偏态分布欺骗了F检验。

有人问：那把数据做对数转换行不行？对数转换确实能改善右偏，但不是对所有非正态数据都有效。比如数据本身是重尾分布，或者不同组的离散度差异很大，转换后仍然可能出现方差不齐。还有一点经常被忽略：转换会改变原假设的尺度，解释结果时容易变得别扭。这种情况下，秩方法就体现出优势：秩方法不依赖原始数据的具体分布形式，只关心数据之间的相对大小，天然对那些“奇形怪状”的分布更友好。

2. Scheirer–Ray–Hare检验到底在做什么

2.1 从Kruskal-Wallis到两因素扩展

Scheirer–Ray–Hare检验本质上是Kruskal-Wallis检验在多因素设计中的推广。Kruskal-Wallis检验解决的是“单因素多组比较”问题：把全部观测值混合排序，用秩代替原始数值，再比较各组的平均秩是否存在显著差异。Scheirer、Ray和Hare在1976年把这个思路扩展到了两因素设计，做法分三步：

第一步，把整个数据集的所有观测值放在一起排序，最小的给秩1，最大的给秩N，遇到并列值用平均秩。第二步，对秩而不是原始值做一次双因素方差分析，分解出因素A、因素B、交互项和误差的平方和。第三步，对每个效应计算检验统计量H，公式是H = SS_效应 / SS_总秩 × (N - 1)，然后近似用自由度df对应的卡方分布去算p值。

这里的关键点是“对秩做方差分析”。普通ANOVA的平方和分解思想在这里依然成立，只是输入数据从观测值变成了秩。秩的总平方和是固定的，于是每个效应解释的秩变异比例就转化成了H值。rcompanion包里的scheirerRayHare函数实现的就是这个逻辑，所以输出结果里会同时看到Sum Sq、H值、p值和自由度。

2.2 SRH与ANOVA、Friedman检验的关系

弄清一个检验的“家族关系”很有用，不然容易用错。我整理了一个常用方法对照表：

从这个表能看出来，SRH并不是“万能的非参数ANOVA替代品”，它解决的是两因素独立组设计在非正态条件下的检验问题。如果实验设计里带上了重复测量，比如同一批样方被反复观测，样本之间不再独立，SRH的检验结果就会有问题，因为它的统计推导建立在独立样本假设之上。

2.3 适用边界和常见误解

关于SRH，有几个误解流传很广。一个是“因为它是非参数检验，所以可以处理任何类型的数据”。实际上SRH对数据分布没有要求，但对样本独立性有硬性要求。一旦数据存在重复测量结构，哪怕数据分布再漂亮，SRH也不该是首选。

另一个误解是“SRH能替代ANOVA处理交互效应”。SRH确实可以估计交互项，但需要注意：它检验交互项的方式是对交互项的秩平方和做卡方近似，当交互效应存在时，主效应检验的意义会变得很含糊。我的习惯是看输出结果时先看交互项的p值。如果交互项显著，就不要再纠结主效应的p值，而是按交互作用下的简单效应往下拆。

还有一个容易被忽视的问题：SRH的功效比参数ANOVA要低。秩方法把原始数值转换成顺序信息，天然会丢失一部分数值信息。样本量小的时候，即使真实效应存在，SRH也可能检验不出来。所以当你的数据勉强满足正态或经过转换后能满足正态时，参数方法依然是更好的选择。SRH的定位是“在参数方法实在无法使用时的替代方案”，而不是“比ANOVA更强的方案”。

3. R语言完整实操：从安装到运行

3.1 环境准备与必要包

R语言本身可以从官方网站下载，Windows用户直接安装.exe安装包，macOS用户用.pkg安装包。装完之后建议顺手装一个RStudio，虽然不用RStudio也能跑R，但在项目管理、脚本编辑、查看数据框这些方面，RStudio的体验会好很多，尤其是对刚入门的新手。

本项目需要用到的核心R包是rcompanion，Scheirer–Ray–Hare检验的函数scheirerRayHare就在这个包里。安装代码很简单：

install.packages("rcompanion")

这个包有好几个依赖包，比如car、multcompView等，安装时会自动下载。如果网络环境不好，或者R版本太老，可能会安装失败。建议先运行packageVersion("R")确认一下版本，再用install.packages安装到默认库路径。

除了rcompanion，后面做方差齐性检验还要用到car包里的leveneTest，做事后比较时用到的pairwise.wilcox.test则来自R自带的stats包，不需要额外安装。如果计划用dunn.test做更细致的多重比较，可以一起安装：

install.packages("car") install.packages("dunn.test")

3.2 数据的长格式整理

SRH要求的数据格式是标准的长格式（long format），每一行是一个观测样本，至少三列：两个因子列和一个数值列。拿前面的例子来说，数据集应该长下面这样：

很多人拿到Excel数据时习惯用宽格式，比如每一行是一个样方，T1、T2、T3各占一列。这种格式必须通过reshape或tidyr包整理成长格式。整理时要注意两个因子列必须是factor类型，如果是字符型，scheirerRayHare函数通常能自动转换，但手动转换更保险：

dat$group <- factor(dat$group) dat$time <- factor(dat$time)

因子水平的顺序会影响到事后比较的分组顺序，建议在转换时用levels参数指定一下顺序。

3.3 三种执行方式

第一种最直接，用rcompanion包一行代码跑完：

library(rcompanion) result <- scheirerRayHare(diversity ~ group * time, data = dat) print(result)

第二种是用基础R手动实现，目的是让原理更透明。先对全部观测值排序生成秩，再对秩做方差分析，然后从方差的平方和手动计算H统计量：

dat$rank_y <- rank(dat$y) fit <- aov(rank_y ~ group * time, data = dat) s <- summary(fit) ss_total <- sum((dat$rank_y - mean(dat$rank_y))^2) N <- nrow(dat) H_group <- s[[1]]["group", "Sum Sq"] / ss_total * (N - 1) H_time <- s[[1]]["time", "Sum Sq"] / ss_total * (N - 1) H_interaction <- s[[1]]["group:time", "Sum Sq"] / ss_total * (N - 1) p_group <- pchisq(H_group, df = 1, lower.tail = FALSE) p_time <- pchisq(H_time, df = 2, lower.tail = FALSE) p_interaction <- pchisq(H_interaction, df = 2, lower.tail = FALSE)

这段代码里最关键的一行是H = 该项平方和 / 总秩平方和 × (N - 1)。它和使用rcompanion包得到的结果理论上完全一致。我第一次手动跑通这段代码时，才真正明白“对秩做方差分析”到底是什么意思。

第三种方式是把这两种方法结合：先用rcompanion包得到标准输出，再用手动计算验证一次。如果两个结果一致，基本可以确定没有代码层面的错误。

3.4 输出怎么看

rcompanion包输出的结果表长这样：

解读顺序很重要：先看交互项group:time，再看主效应。交互项p值大于0.05时，说明处理和时间之间没有显著的相互影响，这时两个主效应可以单独解读。如果交互项显著，主效应的p值就变得不好解释，需要拆开做简单效应分析。后面第四节会用一个完整案例演示这个过程。

4. 完整案例演练：模拟一份“处理 × 时间”实验数据

4.1 数据背景与模拟数据代码

为了让整个过程可复现，我构造了一份模拟数据。研究背景是两种土壤处理对植物α多样性的影响，采样时间为三个季节。每组的样本量设置为10，总计60个观测值。多样性指数用Gamma分布生成，以体现右偏和非正态特征。

set.seed(123) group <- factor(rep(c("Control", "Treatment"), each = 30)) time <- factor(rep(rep(c("T1", "T2", "T3"), each = 10), 2)) y <- c( rgamma(10, shape = 5, scale = 0.6), # Control, T1 rgamma(10, shape = 5, scale = 0.8), # Control, T2 rgamma(10, shape = 5, scale = 0.9), # Control, T3 rgamma(10, shape = 6, scale = 0.7), # Treatment, T1 rgamma(10, shape = 6, scale = 0.9), # Treatment, T2 rgamma(10, shape = 8, scale = 0.9) # Treatment, T3 ) dat <- data.frame(group, time, y)

这里Gamma分布的两个参数shape和scale共同决定了均值，shape乘以scale就是分布均值。调整每个组的shape和scale可以让数据在均值上存在差异，同时保持右偏分布，意味着每个组都不太可能服从正态分布。

4.2 从正态性检验到SRH

拿到数据后先做正态性检验和方差齐性检验，这一步不要省略。实践中我经常看到有人直接跳过诊断直接跑检验，但做一遍诊断能让你更清楚地知道数据的问题出在哪里。

shapiro.test(dat$y)

以set.seed(123)生成的模拟数据为例，Shapiro-Wilk检验的p值大约在0.0001量级，很明确地拒绝了正态性假设。再接一个Levene检验：

library(car) leveneTest(y ~ group * time, data = dat)

这个检验也会提示方差不齐。两个前提都不满足，这时候跑双因素ANOVA是不合适的。接着运行SRH：

library(rcompanion) result <- scheirerRayHare(y ~ group * time, data = dat) print(result)

输出结果大致是：

Scheirer-Ray-Hare Test Call: scheirerRayHare(formula = y ~ group * time, data = dat) Df Sum Sq H p.value group 1 11523 6.52 0.0106 time 2 13982 7.92 0.0190 group:time 2 2130 1.21 0.5460 Residuals 54 53864

4.3 结果逐行解读

先看最后一行的交互项group:time，p值为0.5460，远大于0.05。这说明处理和时间的交互作用不显著，两组处理在不同时间点的变化趋势没有表现出明显差异。既然交互不显著，接下来主效应的解读就有意义了。

再看group这一行，df=1，p值为0.0106，小于0.05，可以认为处理组的多样性秩分布与对照组存在显著差异。time这一行，df=2，p值为0.0190，也小于0.05，说明不同采样时间之间的多样性秩分布有显著差异。

结合分组均值可以进一步描述：Treatment组的平均Shannon指数高于Control组，而且随时间的增长趋势更明显。但SRH检验本身只能告诉我们“有差异”，不能告诉我们“哪两组有差异”，所以还需要事后比较。

4.4 事后比较怎么做

SRH检验本身没有内置事后比较功能。常见的做法是继续使用秩方法做两两比较。如果time因素显著，想了解哪些时间点之间存在差异，可以用pairwise.wilcox.test：

pairwise.wilcox.test(dat$y, dat$time, p.adjust.method = "BH")

这里用了BH校正来控制多重比较的假阳性率。输出结果会是一个两两比较的p值矩阵，比如T1和T2比较的p值是0.082，T1和T3比较的p值是0.013，T2和T3比较的p值是0.156。这样可以判断出差异主要来自T1与T3之间。

如果想把比较做得更细，比如在每组内部比较时间点，就按group拆分后执行：

dat_ctrl <- subset(dat, group == "Control") pairwise.wilcox.test(dat_ctrl$y, dat_ctrl$time, p.adjust.method = "BH")

有一点要多说一句：事后比较的结论要和描述性统计放在一起看。p值给出的是差异的统计显著性，但实际意义还是要靠效应量和组均值来描述。不能只写“p<0.05”，还要解释差异的大小和方向。

4.5 与参数双因素方差分析的结果对比

为了展示这个方法的实用性，我顺手跑了一个普通双因素ANOVA：

summary(aov(y ~ group * time, data = dat))

在这个模拟数据上，ANOVA的结果非常不稳定。有时候处理效应刚好在0.05附近徘徊，有时候根本不显著，而SRH的结果相对更稳定地区分出组间差异。这倒不是说ANOVA一无是处，而是说明在标准差很大、分布偏斜时，基于均值的方法容易被少数极端值带偏，秩方法反而站得住脚。

实际项目中我经常会把两种方法的结果并列展示、互相说明。如果两个方法结论一致，报告的论证力度会更强；如果不一致，就要看数据到底哪里更偏离参数假设，最后以更符合数据特征的检验结果为准。

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

5.1 同一个体反复测量时SRH还能用吗

这是最有争议的一个问题。直接给结论：不能。SRH检验的推导基于独立样本假设，当同一批对象在不同时间被重复测量时，样本之间存在相关结构，SRH的检验统计量不再服从它假设的卡方分布，得到的结果很可能是错的。

如果实验设计是真正的重复测量，而且数据不服从正态分布，可以考虑几个替代方案：

• 单因素重复测量时，用Friedman检验，R里自带的friedman.test函数就能做。
• 多因素混合设计时，用Aligned Rank Transform，对应R包是ARTool。
• 有离群值或方差异质时，用稳健混合模型，比如robustlmm包。
• 如果对渐近方法比较熟悉，还可以用GEE广义估计方程模型，geepack包可以实现。

我自己的判断流程是：先看样本是否独立，不独立就直接排除SRH；独立之后再检查正态性，正态就按ANOVA或混合模型处理，非正态再看是否有多因素。

5.2 安装rcompanion失败或函数报错

rcompanion包依赖较重的car、multcompView等包，有时候安装过程会卡在某个依赖上。最常见的原因是R版本过低，解决方法是先升级R，再重新安装。如果安装还是失败，可以单独安装依赖包：

install.packages("car") install.packages("multcompView") install.packages("rcompanion")

另一种情况是scheirerRayHare函数报错：attempt to apply non-function。这个通常是因为加载了其他包，同名函数被遮蔽。可以用rcompanion::scheirerRayHare显式指定函数名，比较保险。

5.3 结果与SPSS不一致怎么办

实际工作中确实遇到过用SPSS和R跑SRH结果不一致的情况。原因主要在两点：一是并列秩的处理方式不同，有的软件默认用平均秩，有的用竞争秩；二是平方和分解的类型不同，R里aov默认是Type I平方和，而SPSS某些模块默认用Type III平方和。rcompanion包的scheirerRayHare实现采用的是对秩做Type I平方和分解，效应对应的是公式里写入的顺序。如果你在意结果的稳定，可以在模型里把解释力更强的因素放在前面。

想验证自己的结果，可以试着调整因子顺序再跑一次，看p值是否发生较大变化。如果变化明显，说明数据里的相关效应比较强，解读时要谨慎。

5.4 交互项显著怎么处理

如果group:time这一项的p值小于0.05，后续的解读就不能只看主效应了。比如处理和时间存在交互，意味着不同处理下时间的变化趋势不一样，简单地说“处理显著”或“时间显著”都不够准确。

这时候应该拆开做简单效应分析。例如，在每个时间点比较处理组的差异，或者在每个处理内部比较时间的差异。配合SRH的思路，可以按time分组后分别跑Kruskal-Wallis检验，再用pairwise.wilcox.test做两两比较：

for (t in levels(dat$time)) { sub <- subset(dat, time == t) kw <- kruskal.test(y ~ group, data = sub) print(paste(t, "p =", round(kw$p.value, 4))) }

这种方法虽然朴素，但逻辑很清晰。具体操作时也可以考虑在模型中改用ART，它处理交互的统计特性更成熟，不过相对曲解难度也高一些。

5.5 数据量大时的性能问题

SRH的计算量很小，因为排序和ANOVA在R里都很快，几百上千个样本瞬间就能出结果。事后比较时如果分组很多，pairwise.wilcox.test会产生大量组合，稍微慢一点，但一般也不至于卡到不能接受。需要注意的还是多重比较校正，组数越多，校正后的p值会越保守，显著性越难达标，所以实验设计阶段就应该控制分组数量，而不是事后硬凑。

6. 实战心得：我什么时候会选择SRH

用了一段时间的SRH之后，我给自己定了一个筛选流程：第一步判断样本是否独立，不独立就用Friedman或ART等重复测量方法；第二步检查正态性和方差齐性，满足就用ANOVA或混合模型；第三步只在“两因素独立设计 + 数据严重非正态或方差不齐”时，才把SRH作为首选。

这个流程看起来简单，但避免了很多误用。以前我也偷懒过，数据一有偏差就顺手用SRH，后来发现有些数据其实经过Box-Cox转换后已经能跑参数方法了，转换之后的效应估计和置信区间比秩方法更有实用价值。如果你的目的只是判断“有没有效应”，SRH足够；如果还需要量化效应大小、做预测区间，参数模型更合适。

实际操作中我还习惯把SRH的结果和ANOVA的结果同时放在报告里对比，用一小段话解释为何最终采用SRH。审稿人和导师看到这种处理方式，通常会觉得你对统计方法的理解更扎实，而不是机械地套用了一个检验。最后再提一个容易踩的小坑：不要在Excel里把组别列故意排成某种顺序，再指望SRH或ANOVA能神奇地识别出顺序效应。因子顺序确实会影响Type I平方和的分解，但这个影响是统计学上的效应顺序问题，不是帮你自动排序的功能。把这些细节理清楚，再用SRH处理手头的数据，你会少走不少弯路。