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

軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計(jì)學(xué)咨詢(xún)中心(100850) 柳偉偉 胡良平 周詩(shī)國(guó)

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

原理與方法

四格表的一般形式可見(jiàn)表1:

表1 四格表的一般形式

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

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

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

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

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

由上可知，對(duì)于每一個(gè)特定的π，都可算得一個(gè)P(π)?，F(xiàn)在的核心問(wèn)題是:如何去掉未知總體率π。在非條件精確檢驗(yàn)中，首先讓?duì)性谄淇赡艿娜≈捣秶?，1〕內(nèi)變動(dòng)，針對(duì)每個(gè)π 計(jì)算P(π);然后取 P(π)的上確界作為檢驗(yàn)的最終P值。這樣就達(dá)到了去除冗余參數(shù)π的目的，它能夠保證實(shí)際犯Ⅰ類(lèi)錯(cuò)誤的概率始終不會(huì)超過(guò)檢驗(yàn)水準(zhǔn)〔6〕。非條件精確檢驗(yàn)的P值可表示如下:

在上式中，總體率π是在〔0，1〕內(nèi)變動(dòng)的，然而，這個(gè)范圍內(nèi)的有些取值并不被觀(guān)察到的數(shù)據(jù)所支持，也就是說(shuō)由觀(guān)察數(shù)據(jù)來(lái)判斷的話(huà)，總體率取這些值的可能性很小?；诖耍珺erger和Boos提出讓?duì)性谄?－γ可信區(qū)間內(nèi)變動(dòng)，γ取較小的值，例如0.001或0.0001，該可信區(qū)間是通過(guò)觀(guān)察到的四格表數(shù)據(jù)計(jì)算

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

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

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

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

具體程序如下:的〔7〕。此時(shí)非條件精確檢驗(yàn)的P值為:

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

實(shí)例應(yīng)用

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

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

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

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

討 論

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

筆者的程序不但可以作為SAS軟件四格表資料統(tǒng)計(jì)分析的補(bǔ)充，經(jīng)過(guò)適當(dāng)改進(jìn)，還可以用于不同方法之間的比較研究，以及實(shí)現(xiàn)其他一些定性資料的精確檢驗(yà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)．完全隨機(jī)設(shè)計(jì)兩樣本率比較的非條件確切檢驗(yàn)方法．中國(guó)衛(wèi)生統(tǒng)計(jì)，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．