四格表資料非條件精確檢驗的SAS實現

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

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

原理與方法

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

表1 四格表的一般形式

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

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

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

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

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

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

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

式中Cγ表示π的1－γ可信區間。

SAS程序實現

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

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

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

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

實例應用

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

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

運行非條件精確檢驗的SAS程序，程序寫法為:

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

討 論

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

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

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

2．Lydersen S，Fagerland 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．韓宏，王彤，何大衛．完全隨機設計兩樣本率比較的非條件確切檢驗方法．中國衛生統計，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．