四格表資料非條件精確檢驗的SAS實現(xiàn)

軍事醫(yī)學科學院生物醫(yī)學統(tǒng)計學咨詢中心(100850) 柳偉偉 胡良平 周詩國

對于四格表資料的分析，在大樣本時通常選用χ2檢驗，小樣本時一般選擇Fisher精確概率法。Fisher精確概率法屬于條件精確檢驗，它通過將四格表的行、列合計固定，從而去掉冗余參數(shù)，也就是未知的總體率π。與之相對的是非條件精確檢驗，該方法由Barnard提出〔1〕，它被認為比Fisher精確檢驗具有更高的檢驗效能〔2〕。在SAS軟件中，通過FREQ過程能夠很容易的實現(xiàn)χ2檢驗與Fisher精確檢驗，但是現(xiàn)有的SAS過程無法實現(xiàn)非條件精確檢驗。本文根據(jù)非條件精確檢驗的原理，編制出了非條件精確檢驗的SAS實現(xiàn)程序。

原理與方法

四格表的一般形式可見表1:

表1 四格表的一般形式

設(shè)A1和A2分別對應于兩個來自不同總體的獨立樣本，樣本量分別為n1和n2;B1和B2分別表示陽性結(jié)果與陰性結(jié)果;n代表總例數(shù)。想考察兩個總體率π1與π2之間的差異是否有統(tǒng)計學意義，原假設(shè)為H0:π1=π2=π。根據(jù)二項分布原理，在原假設(shè)成立的條件下，表1出現(xiàn)的概率為:

Fisher精確檢驗通過固定行、列合計，可以將上式中未知的π去掉，同時將表1出現(xiàn)的概率表達為超幾何概率的形式。然后在行、列合計固定的約束條件下，計算頻數(shù)變動時的各個四格表出現(xiàn)的超幾何概率，求出等于觀察值或比觀察值更極端的四格表出現(xiàn)的概率之和，就得到了檢驗的P值。該檢驗中需要計算概率的四格表的數(shù)目 k=min(n1，n2，m1，m2)+1。

在非條件精確檢驗中，僅僅假定兩個獨立樣本的例數(shù)n1與n2固定，陽性與陰性結(jié)果出現(xiàn)的頻數(shù)m1和m2是不固定的。對于一個給定的總體率π，在行合計固定而列合計變動的條件下，根據(jù)式(1)計算頻數(shù)變動時各個四格表出現(xiàn)的概率，求出等于觀察值或比觀察值更極端的四格表出現(xiàn)的概率之和，記作P(π)。此時需要計算概率的四格表的數(shù)目k=(n1+1)×(n2+1)。對于雙側(cè)檢驗，P(π)可以表示為:

式中T代表衡量四格表極端程度的統(tǒng)計量，T0代表觀察到的四格表所對應的統(tǒng)計量值。與Fisher精確檢驗類似，在非條件精確檢驗中，T的定義可以有多種方式〔4〕。其中計分統(tǒng)計量可表示如下〔5〕:

式中 p1=a/n1，p2=c/n2，p=(a+c)/n。如果 a=c=0 或 b=d=0，則 Zp=0。

由上可知，對于每一個特定的π，都可算得一個P(π)。現(xiàn)在的核心問題是:如何去掉未知總體率π。在非條件精確檢驗中，首先讓π在其可能的取值范圍〔0，1〕內(nèi)變動，針對每個π 計算P(π);然后取 P(π)的上確界作為檢驗的最終P值。這樣就達到了去除冗余參數(shù)π的目的，它能夠保證實際犯Ⅰ類錯誤的概率始終不會超過檢驗水準〔6〕。非條件精確檢驗的P值可表示如下:

在上式中，總體率π是在〔0，1〕內(nèi)變動的，然而，這個范圍內(nèi)的有些取值并不被觀察到的數(shù)據(jù)所支持，也就是說由觀察數(shù)據(jù)來判斷的話，總體率取這些值的可能性很小。基于此，Berger和Boos提出讓π在其1－γ可信區(qū)間內(nèi)變動，γ取較小的值，例如0.001或0.0001，該可信區(qū)間是通過觀察到的四格表數(shù)據(jù)計算

式中Cγ表示π的1－γ可信區(qū)間。

SAS程序?qū)崿F(xiàn)

根據(jù)上文中所介紹的原理，本文針對雙側(cè)檢驗編制了非條件精確檢驗的SAS實現(xiàn)程序，對于單側(cè)檢驗，只要稍加改動就可以。由于Z2p=χ2，為了使程序更為簡潔，在編寫該程序的時候使用χ2統(tǒng)計量代替Zp統(tǒng)計量，通過χ2值衡量四格表極端程度。另外，我們所編制的程序計算出的P值是與式(5)相對應的，如果讀者想讓總體率π在〔0，1〕內(nèi)變動，根據(jù)式(4)來計算檢驗的P值，只需將程序中進行區(qū)間估計的內(nèi)容去掉就可以了。

實現(xiàn)四格表非條件精確檢驗的程序可以分為5個部分:第一部分用于估計總體率π的1－γ可信區(qū)間;第二部分用于計算觀察到的四格表所對應的χ2值;第三部分用于確定等于觀察值或比觀察值更極端的四格表;第四部分用于計算給定總體率π時各個四格表的概率，求出等于觀察值或比觀察值更極端的四格表出現(xiàn)的概率之和;第五部分用于確定P(π)的上確界，計算最終的P值。

具體程序如下:的〔7〕。此時非條件精確檢驗的P值為:

程序中的宏參數(shù)＆a、＆c分別代表兩組的陽性例數(shù)，也就是表1中的a和c;＆n1、＆n2分別代表兩組的樣本例數(shù)，＆n代表總例數(shù)，＆gamma代表顯著性水平;＆ite代表總體率的變動范圍被劃分成的小區(qū)間的個數(shù)，也就是π的取值個數(shù)減1。

實例應用

研究乙肝免疫球蛋白預防胎兒宮內(nèi)感染HBV的效果，將33例HbsAg陽性孕婦隨機分為預防注射組和非預防組，結(jié)果見表2。問兩組新生兒的HBV總體感染率有無差別?

表2 兩組新生兒HBV感染率的比較

運行非條件精確檢驗的SAS程序，程序?qū)懛?

該例中非條件精確檢驗的的P值為0.11778，此時對應的總體率π的取值為0.26727。對于本例，也可以算得Fisher精確檢驗的P值為0.1210，該值略大于非條件精確檢驗的結(jié)果。

討 論

為了驗證本文中SAS程序的正確性，把該程序應用于多個實例，將產(chǎn)生的結(jié)果與Statxact軟件的分析結(jié)果進行了反復比較，證明了本程序的運行結(jié)果是可靠的〔8〕。此外，我們也對文獻〔6〕中的例子進行了計算，結(jié)果與該文也是相同的。從運算速度來看，該程序也是比較令人滿意的:對于小樣本和中等規(guī)模樣本的資料，能夠很快得出運行結(jié)果;對于大樣本資料，耗時也不會很久。綜合來看，一般資料的非條件精確檢驗，都可以在數(shù)分鐘之內(nèi)運算完成。

筆者的程序不但可以作為SAS軟件四格表資料統(tǒng)計分析的補充，經(jīng)過適當改進，還可以用于不同方法之間的比較研究，以及實現(xiàn)其他一些定性資料的精確檢驗。

1．Barnard GA．Significance test for 2 ×2 tables．Biometrika，1947，34:123-138．

2．Lydersen S，F(xiàn)agerland MW，Laake P．Recommended tests for association in 2 ×2 tables．Statistics in Medicine，2009，28(7):1159-1175．

3．SAS Institute Inc．SAS/STAT 9．2 User's Guide．Cary，NC:SAS Institute Inc，2008:1675-1829．

4．Sahai H，Khurshid A．On analysis of epidemiological data involving a 2×2 contingency table:an overview of Fisher's exact test and Yates.correction for continuity．Journal of Biopharmaceutical Statistics，1995，5(1):43-70．

5．Suissa S，Shuster JJ．Exact unconditional sample sizes for the 2 ×2 binomial trial．Journal of the Royal Statistical Society，1985，148(4):317-327．

6．韓宏，王彤，何大衛(wèi)．完全隨機設(shè)計兩樣本率比較的非條件確切檢驗方法．中國衛(wèi)生統(tǒng)計，2005，22(4):207-209．

7．Berger RL，Boos DD．P values maximized over a confidence set for the nuisance parameter．Journal of the American Statistical Association，1994，89(427):1012-1016．

8．Statxact．Statxact 9 user's manual．Cytel Co．，2010:471-573．