劉多偉，邱世芳，何 杰

(重慶理工大學(xué) 理學(xué)院， 重慶 400054)

在臨床實(shí)驗(yàn)和醫(yī)學(xué)研究中，對疾病的流行率進(jìn)行估計(jì)是一項(xiàng)重要的研究工作。為了估計(jì)疾病的流行率，需要對感興趣的總體是否患有某種疾病進(jìn)行診斷。為了節(jié)約成本或更快獲得診斷結(jié)果，常使用一種價(jià)格便宜、能迅速獲得結(jié)果但有誤判的篩檢方法進(jìn)行第一次診斷，而利用這種有誤判數(shù)據(jù)去估計(jì)流行率又會(huì)導(dǎo)致有偏的結(jié)果(Bross[1])。為此，沒有誤判的金標(biāo)準(zhǔn)常用于進(jìn)一步診斷。但是，金標(biāo)準(zhǔn)往往價(jià)格昂貴、對被檢驗(yàn)個(gè)體有傷害等不足，因而不適用于每個(gè)個(gè)體。為此，Tenenbein[2]提出了一種折中的抽樣方法——二重抽樣。所謂二重抽樣，是指從感興趣的總體中隨機(jī)抽取N個(gè)個(gè)體，使用篩檢方法對他們進(jìn)行診斷，然后從中再隨機(jī)抽取n個(gè)個(gè)體接受金標(biāo)準(zhǔn)診斷。由于這n個(gè)個(gè)體接受了兩種診斷，它反映了個(gè)體的真實(shí)特征并對篩檢方法的分類結(jié)果進(jìn)行了核實(shí)，因而這種二重抽樣方法獲得的數(shù)據(jù)又被稱之為部分核實(shí)數(shù)據(jù)[3]。

基于部分核實(shí)數(shù)據(jù)的統(tǒng)計(jì)研究已成為重要的研究課題。自Tenenbein[2]提出基于二分類數(shù)據(jù)估計(jì)二項(xiàng)比例的二重抽樣方法以來，已有大量文獻(xiàn)進(jìn)行了研究。如Tenenbein[4]將二重抽樣框架推廣到多項(xiàng)分布數(shù)據(jù)的研究中， Yiu等[5]基于潛在變量模型研究了二重抽樣數(shù)據(jù)的統(tǒng)計(jì)推斷。最近，Tang等[6-7]對于部分核實(shí)數(shù)據(jù)的統(tǒng)計(jì)推斷進(jìn)行了一系列的研究，給出了單組樣本下基于部分核實(shí)數(shù)據(jù)對疾病流行率的假設(shè)檢驗(yàn)、區(qū)間估計(jì)以及樣本量的統(tǒng)計(jì)推斷方法，以及兩組獨(dú)立樣本下基于疾病流行率差的統(tǒng)計(jì)推斷方法[8-9]。當(dāng)金標(biāo)準(zhǔn)不存在時(shí)，基于不完全無誤的“金標(biāo)準(zhǔn)”下的部分核實(shí)數(shù)據(jù)，Qiu等[10-11]給出了疾病流行率的假設(shè)檢驗(yàn)過程、區(qū)間估計(jì)以及樣本量的確定方法；邱世芳[12]研究了在給定的置信水平下，區(qū)間寬度控制在給定范圍內(nèi)的樣本量的估計(jì)公式。但是，疾病流行率可能與個(gè)體的年齡、性別或病人是否吸煙、喝酒等生活習(xí)慣有關(guān)。因此，為了避免這些不同因素帶來的混雜效應(yīng)，常將這些因素看作分層變量考慮分層設(shè)計(jì)下的統(tǒng)計(jì)推斷。對于不同的層，疾病流行率是否有顯著的不同是人們關(guān)心的問題。為此，本文將研究在金標(biāo)準(zhǔn)存在時(shí)基于分層設(shè)計(jì)下的部分核實(shí)數(shù)據(jù)對疾病流行率的齊性檢驗(yàn)過程。

1 模型與參數(shù)估計(jì)及檢驗(yàn)統(tǒng)計(jì)量

1.1 模型與參數(shù)估計(jì)

假設(shè)從第j個(gè)總體中隨機(jī)抽取了Nj個(gè)個(gè)體，每個(gè)個(gè)體都接受篩檢檢驗(yàn)，檢驗(yàn)結(jié)果記為Fj。假設(shè)Fj=1表示該個(gè)體檢驗(yàn)為陽性，否則,Fj=0；不考慮篩檢的結(jié)果和其他因素，從這Nj個(gè)個(gè)體中再隨機(jī)抽取nj(nj

令πj=P(Dj=1)為在第j組中被診斷為患病的概率；篩檢方法的敏感度為ηj=P(Fj=1|Dj=1)；特異度為θj=P(Fj=0|Dj=0)。假設(shè)ηj=η,θj=θ對所有j=1,2,…,J，易知pj=πjη+(1-θ)·(1-πj)。概率結(jié)構(gòu)如表1所示。

表1 第j組的數(shù)據(jù)和概率結(jié)構(gòu)

本文研究如下的假設(shè)檢驗(yàn)：H0:π1=π2=…=πJ?H1:π1,…,πJ不全相等。令m={(n11j,n10j,n01j,n00j,xj,yj)′,j=1,…,J},π=(π1,…,πJ)′，則對數(shù)似然函數(shù)為：

l1=l1(m;π,η,θ)=

A1+n11+log(η)+n10+log(1-η)+

n00+log(θ)+n01+log(1-θ)+

(n01j+n00j)log(1-πj)+

xjlog(pj)+yjlog(1-pj)}

j=1,2,…,J

求解以上方程組，可以得到πj(j=1,…,J),η,θ的非限制性極大似然估計(jì)。其中：

j=1,2,…,J

(1)

(2)

在原假設(shè)H0:π1=π2=…=πJ下，不妨令π1=π2=…=πJ=π,則對數(shù)似然函數(shù)為：

l=l2(m;π,η,θ)=

A2+(n11++n10+)log(π)+

(n01++n00+)log(1-π)+n11+log(η)+

n10+log(1-η)+n00+log(θ)+

n01+log(1-θ)+x+log(p)+y+log(1-p)

其中：A2是常數(shù)；p=p(π,η,θ)=πη+(1-θ)(1-π)。根據(jù)Tang等[3]，易得π、η、θ的限制性極大似然估計(jì)分別為：

1.2 齊性檢驗(yàn)統(tǒng)計(jì)量

對于齊性檢驗(yàn)：H0:π1=π2=…=πJ? H1:π1,…,πJ不全相等。本文考慮了如下檢驗(yàn)統(tǒng)計(jì)量。

1.2.1加權(quán)最小二乘統(tǒng)計(jì)量

(3)

在原假設(shè)成立的情況之下，當(dāng)Nj,nj→∞(j=1,…,J),Twls～χ2(J-1)，其中χ2(J-1)表示自由度為J-1的卡方分布。

1.2.2加權(quán)最小二乘統(tǒng)計(jì)量的對數(shù)變換統(tǒng)計(jì)量

為了使加權(quán)最小二乘統(tǒng)計(jì)量Twls更接近正態(tài)分布，根據(jù)Fisher[14]，考慮對檢驗(yàn)統(tǒng)計(jì)量Twls進(jìn)行對數(shù)變換后的統(tǒng)計(jì)量：

Tlwls={log(Twls/(J-1))/2+

(4)

當(dāng)Nj,nj→∞(j=1,…,J),該檢驗(yàn)統(tǒng)計(jì)量漸近服從標(biāo)準(zhǔn)正態(tài)分布，即當(dāng)Tlwls≥z1-α，則拒絕原假設(shè)，其中z1-α為標(biāo)準(zhǔn)正態(tài)分布的上α分位數(shù)。

1.2.3基于對數(shù)變換的加權(quán)最小二乘檢驗(yàn)統(tǒng)計(jì)量

(5)

1.2.4基于logit變換的加權(quán)最小二乘檢驗(yàn)統(tǒng)計(jì)量

(6)

1.2.5基于雙對數(shù)變換的加權(quán)最小二乘檢驗(yàn)統(tǒng)計(jì)量

(7)

1.2.6Score檢驗(yàn)統(tǒng)計(jì)量

通過非限制性的對數(shù)似然函數(shù)，可得到關(guān)于向量π=(π1,…,πJ)′的一階偏導(dǎo)數(shù)，根據(jù)Score檢驗(yàn)的一般理論(Rao[15])，可以給出檢驗(yàn)H0:π1=π2=…=πJ的Score統(tǒng)計(jì)量為：

(8)

pj=πjη+(1-θ)(1-πj),j=1,…,J。

在原假設(shè)H0成立下，當(dāng)Nj,nj→∞(j=1,…,J),Tsc～χ2(J-1)。

1.2.7似然比檢驗(yàn)統(tǒng)計(jì)量

通過以上給出的對數(shù)似然函數(shù)l1(m;π,η,θ)和H0:π1=π2=…=πJ下的對數(shù)似然函數(shù)l2(m;π,η,θ)，可進(jìn)一步考慮如下的似然比檢驗(yàn)統(tǒng)計(jì)量：

Tl=2{l1(m;π,η,θ)-l2(m;π,η,θ)}

(9)

根據(jù)似然比檢驗(yàn)理論可知，在原假設(shè)成立時(shí)，當(dāng)Nj,nj→∞(j=1,…,J),Tl～χ2(J-1)。

1.3 漸近的檢驗(yàn)過程

設(shè)ti是檢驗(yàn)統(tǒng)計(jì)量的Ti(i=wls,lwls,log,logit,dlog,sc,l)的樣本觀測值，根據(jù)這些檢驗(yàn)統(tǒng)計(jì)量的漸近分布，在顯著性水平α下拒絕H0:π1=π2=…=πJ，當(dāng)

或者

tlwls≥z1-α

2 模擬研究

為了考察本文提出的檢驗(yàn)統(tǒng)計(jì)量的統(tǒng)計(jì)性質(zhì)，通過Monte Carlo模擬方法對基于以上檢驗(yàn)統(tǒng)計(jì)量的各種檢驗(yàn)過程進(jìn)行評價(jià)?？紤]了J=3的模型和如下的樣本量：平衡樣本量設(shè)計(jì)① 小樣本(n1,n2,n3,N1,N2,N3)=(30,30,30,50,50,50)；② 中等樣本(n1,n2,n3,N1,N2,N3)=(50,50,50,100,100,100);③ 大樣本(n1,n2,n3,N1,N2,N3)=(200,200,200,500,500,500)和(4)非平衡樣本(n1,n2,n3,N1,N2,N3)=(40,50,60,100,80,100)。為了考察每個(gè)檢驗(yàn)統(tǒng)計(jì)量犯第一類錯(cuò)誤的概率，在以上設(shè)置的各種樣本量下，考慮參數(shù)設(shè)置：π=0.1,0.3,0.5，η=0.6,0.7,θ=0.75,0.9，即考慮了3(π的值)×2(η的值)×2(θ的值)=12種參數(shù)組合。

為了考察每個(gè)檢驗(yàn)統(tǒng)計(jì)量的經(jīng)驗(yàn)功效，考慮參數(shù)設(shè)置：η=0.6,0.7,θ=0.75,0.9，(π1,π2,π3)=(0.1,0.15,0.2),(0.2,0.35,0.35),(0.3,0.55,0.5)，即考慮了3×2×2=12種參數(shù)組合。對于每個(gè)樣本量和每種參數(shù)組合，隨機(jī)產(chǎn)生5 000組隨機(jī)樣本{(n11j,n10j,n01j,n00j,xj,yj):j=1,2,3}分別計(jì)算本文提出的7種檢驗(yàn)統(tǒng)計(jì)量的值，其中，{(n11j,n10j,n01j,n00j):j=1,2,3}從多項(xiàng)分布(nj:πjη,πj(1-η),(1-θ)(1-πj),θ(1-πj))中產(chǎn)生，{(xj,yj):j=1,2,3}從二項(xiàng)分布(Nj-nj:pj,1-pj)中產(chǎn)生，進(jìn)而估計(jì)每個(gè)檢驗(yàn)統(tǒng)計(jì)量下犯第一類錯(cuò)誤的概率和檢驗(yàn)功效。其中，犯第一類錯(cuò)誤的概率估計(jì)公式為：被檢驗(yàn)統(tǒng)計(jì)量T拒絕的次數(shù)/5 000(H0成立下)，經(jīng)驗(yàn)功效的估計(jì)公式為：被檢驗(yàn)統(tǒng)計(jì)量T拒絕的次數(shù)/5 000(H1成立下)。顯著性水平α=0.05，犯第一類錯(cuò)誤的模擬結(jié)果如表2～5所示。

由于篇幅的限制，本文只列出了中等樣本(平衡設(shè)計(jì)和非平衡設(shè)計(jì))下的經(jīng)驗(yàn)功效的模擬結(jié)果，如表6和表7所示。

表2 小樣本(n1,n2,n3,N1,N2,N3)=(30,30,30,50,50,50)下各種檢驗(yàn)統(tǒng)計(jì)量犯第一類錯(cuò)誤的概率(%)(顯著性水平為α=0.05)

表3 中等樣本(n1,n2,n3,N1,N2,N3)=(50,50,50,100,100,100)下各種檢驗(yàn)統(tǒng)計(jì)量犯第一類錯(cuò)誤的概率(%)(顯著性水平為α=0.05)

表4 大樣本(n1,n2,n3,N1,N2,N3)=(200,200,200,500,500,500)下各種檢驗(yàn)統(tǒng)計(jì)量犯第一類錯(cuò)誤的概率(%)(顯著性水平為α=0.05)

表5 不平衡樣本(n1,n2,n3,N1,N2,N3)=(40,50,60,100,80,100)下各種檢驗(yàn)統(tǒng)計(jì)量犯第一類錯(cuò)誤的概率(%)(顯著性水平為α=0.05)

表6 中等樣本(n1,n2,n3,N1,N2,N3)=(50,50,50,100,100,100)下各種檢驗(yàn)統(tǒng)計(jì)量的檢驗(yàn)功效(%)(顯著性水平為α=0.05)

表7 不平衡樣本(n1,n2,n3,N1,N2,N3)=(40,50,60,100,80,100)下各種檢驗(yàn)統(tǒng)計(jì)量的檢驗(yàn)功效(%)(顯著性水平為α=0.05)

根據(jù)模擬結(jié)果可知，① 即使在小樣本(如(n1,n2,n3,N1,N2,N3)=(30,30,30,50,50,50))下，似然比檢驗(yàn)統(tǒng)計(jì)量都表現(xiàn)很好，其犯第一類錯(cuò)誤的概率都很接近給定的顯著性水平且具有較高的功效；② 除了π很小(如π=0.1)外，即使在樣本量較小的情況下，score統(tǒng)計(jì)量都能很好地控制其犯第一類錯(cuò)誤的概率且具有較高的檢驗(yàn)功效；③ 隨著樣本量的增大，各個(gè)統(tǒng)計(jì)量犯第一類錯(cuò)誤的概率都越來越接近名義水平，特別地，當(dāng)樣本量較大(如(n1,n2,n3,N1,N2,N3)=(50,50,50,100,100,100))時(shí)，Tlog、Tlogit、Tdlog也能夠很好地將犯第一類錯(cuò)誤的概率都控制在名義水平左右。綜上所述，從各個(gè)檢驗(yàn)統(tǒng)計(jì)量犯第一類錯(cuò)誤的概率和檢驗(yàn)功效角度研究，可以看出統(tǒng)計(jì)量Tsc、Tl總是具有較好的統(tǒng)計(jì)性質(zhì)，因而實(shí)際應(yīng)用中推薦使用這兩個(gè)統(tǒng)計(jì)量。當(dāng)樣本量較大時(shí)，Tlog、Tlogit、Tdlog也值得推薦。

3 實(shí)際數(shù)據(jù)分析

為了評價(jià)本文所提出的檢驗(yàn)統(tǒng)計(jì)量的統(tǒng)計(jì)性質(zhì)，本文利用Nedelman[16]中的瘧疾數(shù)據(jù)來說明本文提出方法的有效性。世界衛(wèi)生組織和尼日利亞加爾基當(dāng)?shù)卣疄榱搜芯慨?dāng)?shù)丿懠驳幕疾÷?，分別考慮了調(diào)查的方法、村莊類型、年齡、性別等不同因素對當(dāng)?shù)丿懠不疾÷实难芯?。為了研究不同年齡對瘧疾患病率的影響，本文選取了3個(gè)年齡段(成年人數(shù)據(jù))：即19～28歲、29～43歲、≥44歲的調(diào)查數(shù)據(jù)進(jìn)行分析，數(shù)據(jù)如表8所示。

表8 尼日利亞加爾基瘧疾數(shù)據(jù)

注：+表示陽性，-表示陰性。

4 結(jié)論

為了研究基于分層設(shè)計(jì)下的部分核實(shí)數(shù)據(jù)的疾病患病率是否相等的問題，即對于如下的齊性檢驗(yàn)：H0:π1=π2=…=πJ? H1:π1,…,πJ不全相等，本文提出7種檢驗(yàn)統(tǒng)計(jì)量Twls、Tlwls、Tlog、Tlogit、Tdlog、Tsc、Tl，并通過隨機(jī)模擬研究了各種檢驗(yàn)統(tǒng)計(jì)量的統(tǒng)計(jì)性質(zhì)，即研究了基于這些檢驗(yàn)統(tǒng)計(jì)量的漸近檢驗(yàn)過程的犯第一類錯(cuò)誤的概率和經(jīng)驗(yàn)功效。模擬結(jié)果表明，Tsc、Tl統(tǒng)計(jì)量不論犯第一類錯(cuò)誤還是經(jīng)驗(yàn)功效都表明基于它們的齊性檢驗(yàn)過程都有良好的表現(xiàn)，即它們犯第一類錯(cuò)誤的概率接近名義水平，同時(shí)具有較高的功效，因此，在實(shí)際應(yīng)用中推薦使用。當(dāng)樣本量較大時(shí)，Tlog、Tlogit、Tdlog也能很好地控制犯第一類錯(cuò)誤的概率，因而被推薦使用于實(shí)際應(yīng)用中。

附錄：Score檢驗(yàn)統(tǒng)計(jì)量的推導(dǎo)

由對數(shù)似然函數(shù)l1,可得到Fisher信息陣為：

I(J+1)1=I1(J+1)

則I11(ψ)是J×J矩陣，I12(ψ)是J×2矩陣，I21(ψ)是2×J矩陣，I22(ψ)是2×2矩陣。

在原假設(shè)H0成立之下的Score檢驗(yàn)統(tǒng)計(jì)量為