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.北京中醫藥大學