Friedman M检验平均秩的多重比较在SAS软件的实现

陈国勇 鲁梦倩 刘仁权

·计算机应用·

Friedman M检验平均秩的多重比较在SAS软件的实现

陈国勇1鲁梦倩2刘仁权2

对于随机区组设计资料不满足方差分析的条件时，可使用基于秩的非参数检验。统计软件SPSS 18.0及其更高版本已经提供了通过菜单操作实现Friedman M检验平均秩的多重比较功能[1]，但SPSS提供的多重比较方法对P值进行了校正，虽降低了犯Ⅰ类错误的概率，但当处理组数较多时，检验结果偏保守。《中国卫生统计》杂志于2013年8月发表了一篇关于在SPSS中通过编程实现同样功能的文章[1]，虽解决了这个问题，但部分计算还需要手工计算(如相同秩次的重复次数的计算)，且只能计算出多重比较的统计量q值，而对应的P值还需要查q界值表[2]，不仅繁琐，且无法得到精确概率。本研究使用SAS编程避免了所有手工计算，同时可获得q值对应的确切概率值，避免了以上不足。

q检验原理

设随机区组设计资料有g个处理组，n个区组，在不满足方差分析的条件时，进行Friedman M检验。Friedman M检验的结果为拒绝H0时,可以认为多个总体分布位置不全相同,需要进一步分析具体哪两组总体位置分布不同。此时，可先将原始数据在每个区组内编秩，相同数据取平均秩,然后基于秩次做多个相关样本两两比较的q检验[3]。q检验的原假设和备择假设分别为H0：任意两个总体分布位置相同，H1：任意两个总体分布位置不同。规定检验水准为α，q统计量的计算公式是：

(1)

其中MS误差=

(2)

Ri为第i个处理组的秩和，Rj为第j个处理组的秩和，n为区组个数，g为处理组个数，tk为各区组内第k个相同秩的个数。q的分布形态依赖于自由度和样本跨度，其中自由度v=(n-1)(g-1)，样本跨度m为将g个样本秩和排序后，Ri和Rj之间涵盖的秩和的个数(包括Ri和Rj自身在内，故2≤m≤g)。

通过式(1)可以计算任意两组的q值和自由度以及样本跨度，利用特定的SAS函数可进一步计算出q值对应的P值[4]，从而得到相应的检验结果。

实例分析

8名受试者在相同试验条件下接受4种不同频率声音的刺激，他们的反应率(%)资料见表1[3]。问4种不同频率声音刺激的反应率是否有差别？

q检验的步骤和相应的SAS程序

第一步：建立数据集，并对各处理组进行正态性检验。

data example;

input x group block @@;

/* x:反应率，group:分组，block:区组*/

cards;

8.4 1 1 9.6 2 1 9.8 3 1 11.7 4 1

11.6 1 2 12.7 2 2 11.8 3 2 12.0 4 2

9.4 1 3 9.1 2 3 10.4 3 3 9.8 4 3

⋮

7.8 1 8 8.2 2 8 8.5 3 8 10.8 4 8

;

proc univariate normal data=example;

var x;class group;run;

正态性检验结果：当group=2时，Shapiro-Wilk统计量(W= 0.81044)对应的P值为0.037，不满足正态分布，使用非参数检验。

第二步：随机区组设计的Friedman M检验。

proc freq data=example;

tables block*group*x/scores=rank cmh2;

run;

Friedman M检验结果：当g=4且n>5时，Friedman M检验的统计量M的分布近似χ2分布，χ2=15.1519，对应的P值为0.0017，可认为4种频率声音刺激的反应率总体分布位置不全相同，需使用多重比较进行进一步分析。

第三步：对反应率x在各区组内编秩，为多重比较作准备。

proc rank data=example out=ex_rank;

by block;var x;ranks x_rank;

/*对x在区组内编秩的结果存放在新变量x_rank中*/

run;

第四步：基于新变量x_rank计算多重比较的q统计量和对应的P值。

data ex_rank1;

/*逐步累加求出4个处理组的秩和，放在数组R的4个元素中*/

setex_rank;

array R[4](0,0,0,0);

retain R;sum_Ri=0;

do i=1 to 4;

if group=i then R[i]=R[i]+x_rank;

sum_Ri=sum_Ri+R[i]*R[i];

/*sum_Ri存放各个处理组的秩的平方和*/

end;run;

data ex_rank2;/*去掉不用的变量和观测*/

set ex_rank1;drop x i;

where group=4 and block=8;run;

ods output table=again(where=(x_N>1) keep=x_N);

/*将各个区组内编相同秩的个数放在数据集again中*/

proc tabulate data=ex_rank1;

class block x_rank;var x;

table block*x_rank,(x),(n);run;

data again_1;

set again;

retain sum_tj 0;

tj=x_N**3-x_N;

sum_tj=sum_tj+tj;

drop x_N;run;

data again_2;

/*只保留数据集中最后汇总的sum_tj*/

set again_1 end=last;

if last=1;drop tj;run;

data last;

merge ex_rank2 again_2;

ms=(block*group*(group+1)*(2*group+1)/6-sum_Ri/block-sum_tj/12)/((block-1)*(group-1)); /*ms用来计算公式(2)中的MS误差*/

q12=abs(R1-R2)/(block*ms)**0.5;

P12=1-probmc(“range”,q12,.,(block-1)*(group-1),2);

/*q12和p12分别用来计算频率A和频率B比较的q值和P值 */

run;

同理可以得到其他组多重比较的结果，最终结果见表2。

表3是之前的学者对相同的实例用SPSS计算出来的结果。其中第2列的P值是通过SPSS编程，然后查q界值表得到的概率，第3列是通过SPSS菜单操作得到的校正概率。

通过表2和表3的比较可以发现，SAS计算的结果和SPSS编程方式得到的结果是一致的，但SPSS编程方式还需要查表，只能得到一个范围，而SAS计算结果更精确。SPSS菜单方式也可以得到精确的概率值，但它的检验功效相对来说较低(频率A与C比较以及频率B与D比较用SAS计算出来都是有差异的，而用SPSS菜单方式得到的检验结果显示差异没有统计学意义)。

讨 论

本研究提供的Friedman M检验平均秩的多重比较的SAS程序，在原始数据集的基础上直接计算出Friedman M秩检验的结果，以及多重比较的结果，中间不需要任何手工计算，在计算q统计量对应的概率时，也不需要查q界值表。此方法不仅可以得到精确的概率值，还提高了统计人员在使用该方法时的工作效率。但以上程序只针对4个处理组，8个区组设计的资料，当组别数量和区组数量发生变化时，需要在程序中做相应的调整。有兴趣的读者也可将处理组数和区组数设置为宏变量，从而提高程序的适用性。

[1]申希平,祁海萍,刘小宁.Friedman M 检验平均秩的多重比较在SPSS 软件的实现.中国卫生统计，2013，30(4)：611-613.

[2]米术斌,张雷,段一娜,等.SPSS软件进行随机区组设计非参数检验的多重比较.现代预防医学，2009，36(2)：217-219.

[3]孙振球主编.医学统计学.第3版.北京：人民卫生出版社，2010，144-145.

[4]周诗国,胡良平.q临界值、ψ值和λ值的含义及其计算.中国卫生统计，2012，29(1)：27-30.

(责任编辑：刘 壮)

1.中国人民大学统计学院(100038)

2.北京中医药大学