搞统计的人最常遇到的尴尬场景就是数据辛辛苦苦收回来画完箱线图一看歪歪扭扭的偏态分布再跑个Shapiro-Wilk检验p值小到怀疑人生。正常情况下双因素设计第一反应是two-way ANOVA可惜数据不给你这个面子。这时候如果转头去用Kruskal-Wallis又只能处理单个因素交互效应怎么办两个因素之间有没有协同或者拮抗总不能每个因素单独测一遍吧。这就要说到Scheirer–Ray–Hare检验了。这个检验在R语言里用起来其实不复杂核心思想就是把Kruskal-Wallis检验推广到双因素设计先用秩转换把非正态问题处理掉再做一次标准ANOVA。我在生态学、农学还有一部分心理学数据分析里反复用过这条路实测下来比硬着头皮变换数据或者直接忽视正态性假设靠谱得多。这篇博客就把整个流程完整拆开从统计原理讲到R代码实现最后再分享一些踩过的坑适合正在被非正态双因素数据折磨的研究生和从业者。1. 先搞懂Scheirer–Ray–Hare检验到底解决什么问题1.1 从Kruskal-Wallis到双因素检验的由来很多人在入门非参数检验时最先接触的是Wilcoxon秩和检验和Kruskal-Wallis检验。Kruskal-Wallis解决的是单因素多水平的问题被叫做非参数的单因素方差分析。它的思路很直接——把原始观测值全部混在一起排秩基于秩次构造检验统计量从而避开正态性假设。Scheirer、Ray和Hare三位统计学家在1976年的时候做了件很自然的事既然Kruskal-Wallis能处理单因素那把这种排秩的思路扩展到双因素设计行不行结果就是现在说的Scheirer–Ray–Hare检验一般简称SRH检验。它和Friedman检验并列是非参数双因素分析里的两大主力区别在于Friedman主要用于随机区组设计或重复测量设计而SRH适用于双因素完全随机设计而且能汇报交互效应这一点很关键。使用场景其实非常具体。比如我做过一个植物幼苗的实验设置两个种植密度水平、两个灌溉水平测每株幼苗的总生物量。数据分布明显右偏因为干旱胁迫下很多个体长得很差集中在低值区域还有一些长势特别好的离群值把均值拉得很高。这种数据你用参数双因素ANOVA残差正态性检验一跑一个不通过做log变换、Box-Cox变换效果时好时坏而且变换后的解释性也变差了。这种情况下SRH检验是非常合适的选择。1.2 什么情况该用它适用场景快速判断不是说非正态就一定无脑上SRH它有自己的适用范围和前提条件。我把判断标准整理一下。适用的典型情形因变量是连续变量或者有序分类变量且明显偏离正态分布各组方差不齐即使经过数据变换仍然无法满足参数检验条件研究设计是双因素完全随机设计两个自变量都是分类变量你关心交互效应而不仅仅是单个主效应需要注意的前提条件各组数据的分布形状需要大致相似。SRH检验本质上是检验各组位置参数是否相同如果有的组分布是左偏、有的组是右偏即使中位数相同检验也可能显著这种结果解释起来要格外小心观测必须独立不适用于重复测量或随机区组设计那个场景应该用Friedman检验数据中的“结”也就是相同观测值不宜过多。少量结用平均秩处理没大问题但如果大量数据集中在某几个值上检验结果会偏保守我特意画了个表把常见的几种分析方案放在一起对比方便快速选择检验方法数据类型正态性要求方差齐性要求能否检验交互效应典型场景双因素ANOVA连续变量要求要求能经典参数分析Friedman检验连续/有序不要求不要求一般不适用随机区组、重复测量Scheirer–Ray–Hare检验连续/有序不要求不要求能但功效偏弱双因素完全随机设计Aligned Rank Transform连续/有序不要求不要求能功效更优交互效应是研究重点时关于最后一行ART方法后面我会单独展开。如果研究重点是交互作用比如你想证明“处理A在水平B下才有用、在另一个水平下没用”那SRH检验的检验功效并不是最好更推荐对齐秩变换ART方法。但SRH作为经典的非参数双因素方案胜在计算简单、结果直观、新手容易掌握。2. 核心原理拆解秩转换背后的3步逻辑2.1 三个步骤说清检验的完整流程SRH检验的逻辑链条并不复杂核心就三步。我把每一步操作背后的统计目的也解释一下。第一步把所有观测值混合在一起进行排秩。这里用的是全局秩不是分组内秩。最小的观测值秩为1最大为NN是总样本量。如果出现相同的观测值取平均秩。这一步的目的是把原始数据映射到1到N的整数序列上消除原始数据的具体数值对结果的影响只保留大小顺序信息所以对分布形态和异常值天然不敏感。第二步以排秩结果为因变量对两个因素及其交互效应做一次常规的方差分析。这一步和普通的ANOVA在计算上完全一样只是因变量换成了秩。拟合模型是rank_value ~ factorA * factorB得到的是各效应因素A、因素B、交互项的平方和SS_A、SS_B、SS_AB以及误差平方和SS_error和对应的自由度。第三步是SRH检验最特别的地方。它不直接取ANOVA输出的F统计量而是用每个效应的平方和除以误差均方得到H统计量H_A SS_A / MS_error H_B SS_B / MS_error H_AB SS_AB / MS_error其中MS_error SS_error / df_error。每个H统计量近似服从自由度为对应效应自由度的卡方分布p值通过卡方分布计算。这里有个关键细节必须强调第三步为什么要用MS_error而不是用MS_A因为秩的总变异是固定的具体原因我在下一小节讲。先记住结论SRH检验的统计量是H值不是F值汇报结果的时候不要搞错。2.2 为什么统计量是SS效应/MS误差我一开始接触这个检验的时候最大的困惑就是为什么非要用效应平方和除以误差均方这跟普通ANOVA的F统计量有什么区别。后来查阅资料和实际跑数据后发现这正是SRH检验的精妙之处也是它能够和Kruskal-Wallis检验在数学上衔接的原因。关键在于秩的总平方和是固定不变的。假设总样本量为N将所有观测值排秩后秩的取值是1到N平均秩是(N1)/2。无论原始数据长什么样秩的总平方和都等于SS_total N(N² - 1) / 12举个例子N32时SS_total 32 × (1024 - 1) / 12 2728。这个值只跟样本量有关跟数据本身一点关系都没有。既然总变异固定那么误差平方和SS_error就可以理解为“总变异减去各效应解释掉的变异”。当处理效应显著时意味着处理解释了总变异里的较大部分剩下的误差均方MS_error会比较小。用SS_effect除以MS_error实际上是在测度“处理效应相对于剩余噪声”的倍数关系。这个比值越大说明处理效应越不可能只是随机波动听上去其实和F统计量的逻辑是相似的。更进一步说在单因素的情形下这个H值在数学上恰好等价于Kruskal-Wallis检验的H统计量。这一点很多资料都没有明确指出来但它验证了SRH检验和Kruskal-Wallis检验的理论同源性也解释了为什么H值服从卡方分布而不是F分布的关键原因——因为秩的分布不涉及额外的方差估计参数。2.3 单因素退化的验证SRH与Kruskal-Wallis结果一致理论说得再漂亮不如实际跑一次数据来得有说服力。我强烈建议任何第一次使用SRH检验的朋友先做一个单因素的验证测试确认你的手写代码没有逻辑错误。在R里单因素情况下手写的SRH检验流程和kruskal.test的结果应该非常接近。我写一小段验证代码你可以直接复制运行set.seed(123) group - factor(rep(c(A, B, C), each 15)) value - c( rgamma(15, shape 2, rate 0.8), rgamma(15, shape 3, rate 0.8), rgamma(15, shape 5, rate 0.8) ) # 标准 Kruskal-Wallis 检验 kruskal.test(value ~ group)$statistic # 手动实现 SRH 检验单因素版 rank_value - rank(value) model - aov(rank_value ~ group) summary_model - summary(model) ss_group - summary_model[[1]]$Sum Sq[1] ss_error - summary_model[[1]]$Sum Sq[2] df_error - summary_model[[1]]$Df[2] ms_error - ss_error / df_error H_srh - ss_group / ms_error H_srh我经常在实际数据分析里用这个方法做自查。如果样数据中结不太多两个数值会非常接近。如果有差异多半是平均秩处理方式的问题不太影响结论。这个验证技巧建议保留到自己的分析流程里能在早期发现代码错误。3. R语言实操从数据模拟到完整代码实现3.1 环境准备与模拟数据构造正式开始实操之前先把环境准备好。R和RStudio需要提前安装好如果你刚接触R语言直接从R官网下载对应操作系统的版本装好之后再安装RStudio日常写代码用RStudio体验会好很多。用到的包需要安装这几个install.packages(c(rcompanion, FSA, ggplot2, dunn.test))rcompanion是核心工具包提供了scheirerRayHare这个专门函数一行代码就能完成检验。FSA包里的dunnTest函数负责事后多重比较ggplot2用于可视化。为了完整演示整个流程我模拟一份植物生态学实验数据。研究假设是种植密度和灌溉水平对幼苗总生物量有影响。设计有两个因素密度分两个水平low、high灌溉分两个水平well、drought每个组合8个重复总共32个观测。为了制造非正态效果我用伽马分布生成数据让干旱胁迫组的分布明显右偏并且人为加入一个离群值set.seed(2024) n - 8 biomass - c( round(rgamma(n, shape 6, rate 0.9), 1), round(rgamma(n, shape 2.2, rate 0.8), 1), round(rgamma(n, shape 5.5, rate 1.1), 1), round(rgamma(n, shape 1.8, rate 0.7), 1) ) # 人为制造一个离群值 biomass[12] - biomass[12] * 2.5 density - factor(rep(c(low, high), each 16)) irrigation - factor(rep(rep(c(well, drought), each 8), 2)) data - data.frame(density, irrigation, biomass) str(data)跑完这段代码你能看到数据框是32行3列的结构。先做一次正态性检验让数据“自证清白”shapiro.test(data$biomass)输出的p值通常会远小于0.05这就说明标准差等参数方法的前提条件不满足需要用SRH检验来处理。3.2 手写实现每一步代码都有讲究自己手写实现SRH检验最大的好处是每一步都透明知道中间发生了什么。对于学习统计原理的人来说强烈建议至少完整写一遍。代码如下# 第一步全局排秩 data$rank_value - rank(data$biomass) # 第二步以秩为因变量做双因素ANOVA model - aov(rank_value ~ density * irrigation, data data) anova_result - summary(model) anova_result运行到这里你会得到一个标准的ANOVA表格。里面包含density、irrigation、density:irrigation三个效应的平方和以及Residuals的平方和和自由度。注意这里的F值和p值是ANOVA框架下的不能直接用还需要进一步计算。# 第三步提取中间结果计算H统计量 ss_density - anova_result[[1]]$Sum Sq[1] ss_irrigation - anova_result[[1]]$Sum Sq[2] ss_interaction - anova_result[[1]]$Sum Sq[3] ss_error - anova_result[[1]]$Sum Sq[4] df_density - anova_result[[1]]$Df[1] df_irrigation - anova_result[[1]]$Df[2] df_interaction - anova_result[[1]]$Df[3] df_error - anova_result[[1]]$Df[4] ms_error - ss_error / df_error H_density - ss_density / ms_error H_irrigation - ss_irrigation / ms_error H_interaction - ss_interaction / ms_error p_density - pchisq(H_density, df_density, lower.tail FALSE) p_irrigation - pchisq(H_irrigation, df_irrigation, lower.tail FALSE) p_interaction - pchisq(H_interaction, df_interaction, lower.tail FALSE) results - data.frame( Effect c(density, irrigation, density:irrigation), Df c(df_density, df_irrigation, df_interaction), H round(c(H_density, H_irrigation, H_interaction), 3), p.value round(c(p_density, p_irrigation, p_interaction), 4) ) results输出的结果大致长这样EffectDfHp.valuedensity12.1840.1394irrigation113.5870.0002density:irrigation10.1580.6908从这个模拟结果看灌溉水平对生物量的影响显著种植密度和交互效应都不显著。这和我们设计数据时的预期基本一致。3.3 一行代码方案rcompanion包快速出结果手写实现适合学习原理但实际分析中我更推荐用rcompanion包里的scheirerRayHare函数简洁不容易出错。用法如下library(rcompanion) result - scheirerRayHare(biomass ~ density * irrigation, data data) result输出表格里直接包含Df、Sum Sq、H和p.value这几列和我手写计算的结果完全一致。这里有一点值得注意scheirerRayHare函数默认使用I型平方和也就是依序添加效应。这意味着在平衡设计下结果没问题但如果数据不平衡效应的顺序会影响结果。如果数据不平衡我的建议是使用car包里的Anova函数配合type3参数来做排秩ANOVA再手动计算H值。不过话说回来既然都用到非参数检验了设计阶段最好控制好平衡性避免给自己挖坑。4. 结果解读、事后比较与论文汇报4.1 输出表格怎么读H值、自由度与p值跑完scheirerRayHare之后很多人对着输出表格发愣不知道该怎么往论文里写。我来把解读逻辑梳理一遍。先看p值。p值小于0.05说明该效应显著这是常规判断。接着看H值大小它类似于卡方检验的统计量H值越大代表效应越强。但注意H值不像效应量那样有标准化含义不能直接用H值比较不同因素的解释力大小。关键要看自由度。自由度对应检验该效应用了多少信息量。因素A有a个水平自由度就是a-1因素B有b个水平自由度是b-1交互项自由度是(a-1)(b-1)。自由度会直接影响p值的计算同一H值在不同自由度下显著性可能完全不同。在结果解读上我特别提醒一点当交互效应显著时主效应不能单独解读。比如密度和灌溉的交互p 0.05就说明“密度的效应依赖于灌溉水平”这时候报告“密度显著影响生物量”是不准确的应该进一步分析简单效应也就是固定一个因素看另一个因素在不同水平下的效应。论文里规范汇报的格式我总结如下“由于生物量数据显著偏离正态Shapiro-Wilk检验p 0.001采用Scheirer–Ray–Hare检验进行双因素非参数方差分析。结果显示灌溉水平对幼苗总生物量有显著影响H 13.587df 1p 0.001种植密度的主效应不显著H 2.184df 1p 0.139两因素交互效应不显著H 0.158df 1p 0.691。”4.2 显著之后怎么办Dunn检验做多重比较SRH检验本身只告诉你“有没有差异”不会告诉你“具体哪些组之间有差异”。当某个因素有显著效应且水平数多于2时需要做多重比较。最常用的是Dunn检验它本质上是逐对比较的秩和检验可以对p值做多重比较校正。以模拟数据为例如果irrigation显著但只有两个水平不需要做多重比较——两个水平下直接看图就知道谁高谁低。但如果灌溉水平是三个或四个就要跑Dunn检验library(FSA) dunnTest(biomass ~ irrigation, data data, method bonferroni)Dunn检验会输出两两比较的Z统计量和校正后的p值。method参数可以设置多重比较校正方法常用的有bonferroni、holm和bh。我个人偏好holm方法它在控制错误率和检验功效之间比较均衡。如果交互效应显著情况会复杂一些。这时需要把所有因素水平组合当作一个整体来做比较比如2×2设计就会有4个组两两比较一共6对这时可以用自定义的分组变量来做Dunn检验data$group - interaction(data$density, data$irrigation) dunnTest(biomass ~ group, data data, method bonferroni)这种做法的逻辑是交互显著意味着不能只看主效应要把每个单元格当作独立的处理来看。比较哪两个处理组合之间有差异才能定位交互效应的来源。4.3 可视化建议箱线图与交互效应图统计分析不配图等于没做。非参数检验的可视化原则和参数检验略有所不同因为汇报中位数比汇报均值更有代表性。我常用的第一张图是分组箱线图用中位数和四分位距来展示分布位置差异再加上抖动点显示原始数据分布library(ggplot2) ggplot(data, aes(x density, y biomass, fill irrigation)) geom_boxplot(outlier.color firebrick) geom_jitter(width 0.15, alpha 0.5, size 1.5) stat_summary(fun median, geom point, aes(group irrigation), position position_dodge(0.75), color black, size 3) labs(x 种植密度, y 总生物量g, fill 灌溉处理) theme_classic(base_size 14)这张图能直观展示各组的中位数、离散程度和离群值情况。配合SRH检验结果读者一眼就能理解“灌溉显著”这个结论的具体表现。第二张推荐交互效应图用来直观观察两个因素之间是否存在交互。R基础绘图就能实现interaction.plot( x.factor data$density, trace.factor data$irrigation, response data$biomass, fun median, type b, col c(black, gray), trace.label 灌溉, xlab 种植密度, ylab 生物量中位数g )如果两条线基本平行说明交互不显著如果明显交叉或者角度差异很大说明存在交互作用。这个图比任何统计检验都直观我在审稿的时候最喜欢看图来判断结论是否可信。4.4 论文中怎么写标准报告模板写论文时关于SRH检验的完整汇报至少应该包含以下内容第一说明为什么用这个检验。数据不符合正态性假设或者方差不齐这是选择SRH检验的理由。第二报告检验的关键统计量。每个效应的H值、自由度、p值都需要写清楚表格是最高效的形式。第三汇报描述性统计量。因为是非参数检验推荐报告各组中位数和四分位距而不是均值和标准差。第四如果做了事后比较报告校正方法、比较组的统计量和校正后的p值。下面是一个可以借鉴的结果汇报表格效应H值自由度p值密度2.18410.139灌溉13.58710.001密度×灌溉0.15810.691被小写“Table”或“表”标注后放入论文结果部分就能让审稿人一眼看懂分析逻辑。注意在表格下方添加注释说明使用的是Scheirer–Ray–Hare检验p值由卡方近似得到。5. 常见报错与避坑指南5.1 典型报错和对应处理实际使用过程中我遇到过几种典型报错每次处理方式都大同小异。汇总成一张速查表报错信息可能原因解决方案variable lengths differ向量长度不一致或数据有缺失值用complete.cases()检查数据确保三个向量等长NA/NaN/Inf in foreign function call数据包含缺失值无法排秩在分析前用na.omit()剔除缺失值system is computationally singular某组只有1个观测或观测数过少检查组容量建议每组至少5个观测出现p值全为1样本太少卡方近似极差考虑精确检验或换用置换检验包安装失败R版本过旧或依赖缺失更新R到最新版本重新安装依赖包最常见的是第一个问题。我习惯在跑分析之前先检查数据结构str(data) sum(is.na(data$biomass))这两行代码能避免大半低级错误。再补充一个容易忽略的点scheirerRayHare函数要求输入的数据格式正确公式里因素变量必须是factor类型。如果因素变量是字符型或数值型分析结果可能和预期不同。可以在跑检验前强制转换data$density - as.factor(data$density) data$irrigation - as.factor(data$irrigation)5.2 六个实际数据分析中容易踩的坑第一个坑交互效应显著时直接解读主效应结果。这个我在前面说过但在实际分析中太多人犯这个错误。交互显著意味着某个因素的效果依赖另一个因素的水平此时主效应的H值和p值已经不能直接回答研究问题了正确姿势是先做简单效应分析。第二个坑忽视数据平衡性。SRH检验在平衡数据下表现稳定但数据严重不平衡时I型平方和的顺序依赖问题会被放大。如果各组样本量差异很大建议改用art.con包里的ART-C方法或者使用car包的Anova函数配合type3来重新计算。第三个坑样本量过小。我试过每组只有3个重复的数据跑出来的p值波动非常大不同随机种子下结论都可能反转。SRH检验基于卡方近似小样本下近似效果不好。我的经验是每单元格至少5到8个观测稳妥起见的话最好每组10个以上。第四个坑处理重复测量数据。SRH检验要求观测独立如果你的实验是同一批个体在不同时间点的重复测量那就是随机区组或者重复测量设计应该用Friedman检验而不是SRH检验。这个错误在纵向研究里经常出现。第五个坑报告中缺失效应量。很多人汇报非参数检验结果只给H和p值不给效应量。审稿人通常希望看到效应的大小衡量。SRH检验可以算eta-squared就是效应平方和除以总平方和eta_sq_density - ss_density / (ss_density ss_irrigation ss_interaction ss_error)当然更规范的做法是用秩方差解释比例。这个数值可以用来比较不同因素的相对重要性。第六个坑把所有检验结果完全交给函数。scheirerRayHare确实方便但你不理解内部逻辑遇到特殊数据时很容易被函数输出误导。比如结处理方式、因子水平顺序、NA值处理这些细节都可能影响结果。我强烈建议在具体数据上跑手写版验证一遍确认结果一致后再使用包函数出正式结论。关于ART方法的补充如果你的研究核心就是交互效应SRH检验相对偏保守检验功效不足。我之前做过一个模拟实验人为注入强交互效应后SRH检验的检出率明显低于ART-C方法。所以如果交互效应是研究的重中之重尽量用art.con包里的ART-C替代。它通过对齐和重排流程能在非参数框架下更好地保留下交互效应信息。最后再分享一个实际操作中的小技巧SRH检验跑完之后把排秩结果拟合的ANOVA表格也保存下来后续做敏感性分析可能用得到。有时候数据中有几个极端离群值影响结论你可以对比删除离群值前后的H值和p值变化判断结论的稳健性。这种做法在审稿时经常被问到提前准备好答案能让回复从容很多。做数据分析这些年我最大的体会是统计检验永远只是工具理解工具背后的假设和边界条件比记住几个函数接口重要得多。SRH检验在非正态双因素数据分析里确实是救场利器但它也不是万能钥匙。如果数据条件允许优先考虑参数方法如果非正态数据迫不得已SRH检验是一个可靠的选择如果交互效应是核心研究问题ART方法值得再多花一点时间学习。