RPT對秩和檢驗的改進及Matlab實現

朱 凱 李 悅

秩和檢驗作為一種方便、有效的非參數檢驗方法，在實際中有著廣泛的應用，特別是在總體分布、參數均未知時，可以用來比較兩樣本的均值〔1〕。本文針對秩和檢驗對來自非正態(tài)總體兩樣本均值比較可能遇到準確性問題，介紹可以有效解決這一問題的方法——randomized permutation test(簡稱RPT)，并運用Matlab軟件編寫程序實現該方法。

RPT方法簡介

在實際應用秩和檢驗對兩樣本均數進行比較時，采用的方法為:首先求出樣本所有可能組合下的分布，根據現有樣本計算某一樣本的秩和，記為Robs，稱Robs為秩和的臨界值，然后根據該分布求出大于或等于Robs的頻率，即為p值。該方法稱為exact permutation test(簡稱EPT)。然而當樣本量較大時往往由于組合數太多(當n1=n2=10時組合數為C1020=184756)難以獲得秩和的確切分布，此時EPT方法難以實現。因此，在樣本量較大時，秩和檢驗將檢驗統(tǒng)計量近似服從正態(tài)分布，利用矩法估計求出其參數(下稱正態(tài)近似法)，這將產生一定的誤差。此外，當樣本存在相同秩次較多時，要求對檢驗統(tǒng)計量進行校正，這不僅增加計算難度且進一步增大誤差〔1－2〕。EPT的一個替代方法是:利用軟件產生隨機數，通過基于樣本的大量重復的隨機組合，得到某一樣本秩和的近似抽樣分布，即對所有可能的組合進行隨機抽樣后計算抽樣秩和大于臨界秩和的頻率，稱之為RPT〔3〕，該方法的關鍵步驟如下:

1. 建立假設，確定檢驗水準，H0:μ1= μ2，H1:μ1＞μ2，α =0.05(單側檢驗);

2.計算現有檢驗統(tǒng)計量Robs;

3.構造檢驗統(tǒng)計量R，在H0假設成立的條件下，利用軟件從兩樣本構成的總體中抽取與樣本容量相同的樣本，并計算其秩和，反復進行該步驟得到檢驗統(tǒng)計量的經驗抽樣分布;

實例應用與Matlab實現

為研究兩位化驗員讀得某種液體黏度讀數的差異，現對同一液體進行重復讀數，數據見表1，試判斷兩化驗員讀數是否有差異(α=0.05)〔1〕。

表1 A、B兩位化驗員讀數結果

解:H0:μ1= μ2，H1:μ1＞ μ2，其中 μ1，μ2分別為兩總體的均值。

由題中數據可知n1=10，n2=11，將第一組秩和定為檢驗統(tǒng)計量，Robs=121。用RPT模擬100000次，即從上述兩組的21個數值中有放回地隨機抽取10個數值，重復進行100000次，求這10個數的秩和RA大于Robs的頻率。筆者編寫了相應的Matlab程序實現該方法，程序及注釋如下:

a=［82 73 91 84 77 98 81 79 87 85］;

b=［80 76 92 86 74 96 83 79 80 75 79］;%數據錄入

n1=length(a);n2=length(b);n=n1+n2;%數據量統(tǒng)計

c=［a，b;ones(1，n)〕］;

［，，stats］=ranksum(a，b);%計算秩和

rank1=stats.ranksum;%計算Robs

m=100000;%模擬次數

t=0;rank0=zeros(1，m);%數據清零

for k=1:m

d=randperm(n);%生成1-21個隨機排列的整數

for k0=1:n1

rank0(k)=rank0(k)+c(2，d(k0));% 隨機抽取10個數并求其秩和

end

if rank0(k)＞rank1

t=t+1;%統(tǒng)計超過臨界秩和的數目

end

end

p=t/m;%超過臨界秩和的頻率

通過運行該程序，得到p=P(R≥Robs)=0.21811，p值均大于0.05，故接受H0，認為兩位化驗員所測得的數據無顯著差異。若使用EPT方法進行全排列，所有組合將達到C1021=352716種，且隨著樣本量的增加EPT方法運算次數還將呈幾何數增長，故EPT方法在實踐中是難以實現的。對上述程序進行適當修改后我們可以得出該問題的確切概率p'=P(R≥Robs)=0.2181239，RPT方法與之相比相對誤差只有0.0063626%。若使用傳統(tǒng)方法，求得的概率 p″=0.21891，此時相對誤差達到了0.36%，是RPT方法的近60倍。當P值較小時，傳統(tǒng)方法的相對誤差還將進一步增大。由此我們可知RPT方法較傳統(tǒng)方法減少了誤差，較EPT方法減少了運算次數，且該方法可根據實際情況適當調整抽樣次數，以達到增加精度或減少計算量的目的。

結 論

在樣本量較大、總體分布未知，沒有其他合適方法進行兩樣本均值比較時，可以使用秩和檢驗來進行統(tǒng)計推斷〔4－6〕。在應用RPT方法進行秩和檢驗時需注意如下幾個問題:

1．秩和檢驗作為一種非參數檢驗，由于不依賴資料的分布類型，故適用范圍廣泛，尤其在等級資料的分析中有較高的功效。

2．編秩時相同值要取平均秩次，否則將使秩和的臨界值發(fā)生錯誤，影響最終結果。

3．模擬誤差來源于Monte Carlo模擬抽樣。理論上，無限次的模擬將會完全消除模擬誤差，但顯然這是不可能也沒有必要的。因此確定模擬誤差足夠小并且計算可行的RPT抽樣次數是有效控制模擬誤差的必要步驟〔3〕。

4．由于程序運行中要用到隨機數，最終的結果有一定的誤差，因此建議反復運行程序，待結果相對穩(wěn)定時再下結論。

RPT作為EPT的一種近似方法，具有使用方便、誤差小、執(zhí)行效率高的優(yōu)點，對出現較多相同秩次時處理能力較強，是一種有效提高秩和檢驗效率的好方法。

1．盛驟，謝式千，潘承毅．概率論與數理統(tǒng)計(第二版)．北京:高等教育出版社，1989:118-121．

2．顏杰，李彩霞，等．完全隨機設計兩組t檢驗與秩和檢驗的功效比較．中國衛(wèi)生統(tǒng)計，2004，21(1):10-13．

3．丁元林，孔丹莉．多個樣本及其兩兩比較的秩和檢驗SAS程序．中國衛(wèi)生統(tǒng)計，2002，19(5):313-314．

4．荀鵬程，趙楊，柏建嶺，等．Permutation Test在假設檢驗中的．數理統(tǒng)計與管理，2006，26(5):616．

5．Cai JW，Shen Y．Permutation tests for comparing marginal survival functions with clustered failure time data．Statist．Med，2000，19:2963-2973．6．王試會，徐勇勇．隨機區(qū)組設計資料秩和檢驗的改進方法．中國衛(wèi)生統(tǒng)計，2003，20(4):231-232．