陳國勇 魯夢倩 劉仁權(quán)

·計(jì)算機(jī)應(yīng)用·

Friedman M檢驗(yàn)平均秩的多重比較在SAS軟件的實(shí)現(xiàn)

陳國勇1魯夢倩2劉仁權(quán)2

對于隨機(jī)區(qū)組設(shè)計(jì)資料不滿足方差分析的條件時(shí)，可使用基于秩的非參數(shù)檢驗(yàn)。統(tǒng)計(jì)軟件SPSS 18.0及其更高版本已經(jīng)提供了通過菜單操作實(shí)現(xiàn)Friedman M檢驗(yàn)平均秩的多重比較功能[1]，但SPSS提供的多重比較方法對P值進(jìn)行了校正，雖降低了犯Ⅰ類錯(cuò)誤的概率，但當(dāng)處理組數(shù)較多時(shí)，檢驗(yàn)結(jié)果偏保守?！吨袊l(wèi)生統(tǒng)計(jì)》雜志于2013年8月發(fā)表了一篇關(guān)于在SPSS中通過編程實(shí)現(xiàn)同樣功能的文章[1]，雖解決了這個(gè)問題，但部分計(jì)算還需要手工計(jì)算(如相同秩次的重復(fù)次數(shù)的計(jì)算)，且只能計(jì)算出多重比較的統(tǒng)計(jì)量q值，而對應(yīng)的P值還需要查q界值表[2]，不僅繁瑣，且無法得到精確概率。本研究使用SAS編程避免了所有手工計(jì)算，同時(shí)可獲得q值對應(yīng)的確切概率值，避免了以上不足。

q檢驗(yàn)原理

設(shè)隨機(jī)區(qū)組設(shè)計(jì)資料有g(shù)個(gè)處理組，n個(gè)區(qū)組，在不滿足方差分析的條件時(shí)，進(jìn)行Friedman M檢驗(yàn)。Friedman M檢驗(yàn)的結(jié)果為拒絕H0時(shí),可以認(rèn)為多個(gè)總體分布位置不全相同,需要進(jìn)一步分析具體哪兩組總體位置分布不同。此時(shí)，可先將原始數(shù)據(jù)在每個(gè)區(qū)組內(nèi)編秩，相同數(shù)據(jù)取平均秩,然后基于秩次做多個(gè)相關(guān)樣本兩兩比較的q檢驗(yàn)[3]。q檢驗(yàn)的原假設(shè)和備擇假設(shè)分別為H0：任意兩個(gè)總體分布位置相同，H1：任意兩個(gè)總體分布位置不同。規(guī)定檢驗(yàn)水準(zhǔn)為α，q統(tǒng)計(jì)量的計(jì)算公式是：

(1)

其中MS誤差=

(2)

Ri為第i個(gè)處理組的秩和，Rj為第j個(gè)處理組的秩和，n為區(qū)組個(gè)數(shù)，g為處理組個(gè)數(shù)，tk為各區(qū)組內(nèi)第k個(gè)相同秩的個(gè)數(shù)。q的分布形態(tài)依賴于自由度和樣本跨度，其中自由度v=(n-1)(g-1)，樣本跨度m為將g個(gè)樣本秩和排序后，Ri和Rj之間涵蓋的秩和的個(gè)數(shù)(包括Ri和Rj自身在內(nèi)，故2≤m≤g)。

通過式(1)可以計(jì)算任意兩組的q值和自由度以及樣本跨度，利用特定的SAS函數(shù)可進(jìn)一步計(jì)算出q值對應(yīng)的P值[4]，從而得到相應(yīng)的檢驗(yàn)結(jié)果。

實(shí)例分析

8名受試者在相同試驗(yàn)條件下接受4種不同頻率聲音的刺激，他們的反應(yīng)率(%)資料見表1[3]。問4種不同頻率聲音刺激的反應(yīng)率是否有差別？

q檢驗(yàn)的步驟和相應(yīng)的SAS程序

第一步：建立數(shù)據(jù)集，并對各處理組進(jìn)行正態(tài)性檢驗(yàn)。

data example;

input x group block @@;

/* x:反應(yīng)率，group:分組，block:區(qū)組*/

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;

正態(tài)性檢驗(yàn)結(jié)果：當(dāng)group=2時(shí)，Shapiro-Wilk統(tǒng)計(jì)量(W= 0.81044)對應(yīng)的P值為0.037，不滿足正態(tài)分布，使用非參數(shù)檢驗(yàn)。

第二步：隨機(jī)區(qū)組設(shè)計(jì)的Friedman M檢驗(yàn)。

proc freq data=example;

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

run;

Friedman M檢驗(yàn)結(jié)果：當(dāng)g=4且n>5時(shí)，F(xiàn)riedman M檢驗(yàn)的統(tǒng)計(jì)量M的分布近似χ2分布，χ2=15.1519，對應(yīng)的P值為0.0017，可認(rèn)為4種頻率聲音刺激的反應(yīng)率總體分布位置不全相同，需使用多重比較進(jìn)行進(jìn)一步分析。

第三步：對反應(yīng)率x在各區(qū)組內(nèi)編秩，為多重比較作準(zhǔn)備。

proc rank data=example out=ex_rank;

by block;var x;ranks x_rank;

/*對x在區(qū)組內(nèi)編秩的結(jié)果存放在新變量x_rank中*/

run;

第四步：基于新變量x_rank計(jì)算多重比較的q統(tǒng)計(jì)量和對應(yīng)的P值。

data ex_rank1;

/*逐步累加求出4個(gè)處理組的秩和，放在數(shù)組R的4個(gè)元素中*/

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存放各個(gè)處理組的秩的平方和*/

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);

/*將各個(gè)區(qū)組內(nèi)編相同秩的個(gè)數(shù)放在數(shù)據(jù)集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;

/*只保留數(shù)據(jù)集中最后匯總的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用來計(jì)算公式(2)中的MS誤差*/

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

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

/*q12和p12分別用來計(jì)算頻率A和頻率B比較的q值和P值 */

run;

同理可以得到其他組多重比較的結(jié)果，最終結(jié)果見表2。

表3是之前的學(xué)者對相同的實(shí)例用SPSS計(jì)算出來的結(jié)果。其中第2列的P值是通過SPSS編程，然后查q界值表得到的概率，第3列是通過SPSS菜單操作得到的校正概率。

通過表2和表3的比較可以發(fā)現(xiàn)，SAS計(jì)算的結(jié)果和SPSS編程方式得到的結(jié)果是一致的，但SPSS編程方式還需要查表，只能得到一個(gè)范圍，而SAS計(jì)算結(jié)果更精確。SPSS菜單方式也可以得到精確的概率值，但它的檢驗(yàn)功效相對來說較低(頻率A與C比較以及頻率B與D比較用SAS計(jì)算出來都是有差異的，而用SPSS菜單方式得到的檢驗(yàn)結(jié)果顯示差異沒有統(tǒng)計(jì)學(xué)意義)。

討 論

本研究提供的Friedman M檢驗(yàn)平均秩的多重比較的SAS程序，在原始數(shù)據(jù)集的基礎(chǔ)上直接計(jì)算出Friedman M秩檢驗(yàn)的結(jié)果，以及多重比較的結(jié)果，中間不需要任何手工計(jì)算，在計(jì)算q統(tǒng)計(jì)量對應(yīng)的概率時(shí)，也不需要查q界值表。此方法不僅可以得到精確的概率值，還提高了統(tǒng)計(jì)人員在使用該方法時(shí)的工作效率。但以上程序只針對4個(gè)處理組，8個(gè)區(qū)組設(shè)計(jì)的資料，當(dāng)組別數(shù)量和區(qū)組數(shù)量發(fā)生變化時(shí)，需要在程序中做相應(yīng)的調(diào)整。有興趣的讀者也可將處理組數(shù)和區(qū)組數(shù)設(shè)置為宏變量，從而提高程序的適用性。

[1]申希平,祁海萍,劉小寧.Friedman M 檢驗(yàn)平均秩的多重比較在SPSS 軟件的實(shí)現(xiàn).中國衛(wèi)生統(tǒng)計(jì)，2013，30(4)：611-613.

[2]米術(shù)斌,張雷,段一娜,等.SPSS軟件進(jìn)行隨機(jī)區(qū)組設(shè)計(jì)非參數(shù)檢驗(yàn)的多重比較.現(xiàn)代預(yù)防醫(yī)學(xué)，2009，36(2)：217-219.

[3]孫振球主編.醫(yī)學(xué)統(tǒng)計(jì)學(xué).第3版.北京：人民衛(wèi)生出版社，2010，144-145.

[4]周詩國,胡良平.q臨界值、ψ值和λ值的含義及其計(jì)算.中國衛(wèi)生統(tǒng)計(jì)，2012，29(1)：27-30.

(責(zé)任編輯：劉 壯)

1.中國人民大學(xué)統(tǒng)計(jì)學(xué)院(100038)

2.北京中醫(yī)藥大學(xué)