基于有金標準下的部分核實數據對疾病流行率的齊性檢驗

劉多偉，邱世芳，何 杰

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

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

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

1 模型與參數估計及檢驗統計量

1.1 模型與參數估計

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

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

表1 第j組的數據和概率結構

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

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=1,2,…,J

(1)

(2)

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

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是常數；p=p(π,η,θ)=πη+(1-θ)(1-π)。根據Tang等[3]，易得π、η、θ的限制性極大似然估計分別為：

1.2 齊性檢驗統計量

對于齊性檢驗：H0:π1=π2=…=πJ? H1:π1,…,πJ不全相等。本文考慮了如下檢驗統計量。

1.2.1加權最小二乘統計量

(3)

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

1.2.2加權最小二乘統計量的對數變換統計量

為了使加權最小二乘統計量Twls更接近正態分布，根據Fisher[14]，考慮對檢驗統計量Twls進行對數變換后的統計量：

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

(4)

當Nj,nj→∞(j=1,…,J),該檢驗統計量漸近服從標準正態分布，即當Tlwls≥z1-α，則拒絕原假設，其中z1-α為標準正態分布的上α分位數。

1.2.3基于對數變換的加權最小二乘檢驗統計量

(5)

1.2.4基于logit變換的加權最小二乘檢驗統計量

(6)

1.2.5基于雙對數變換的加權最小二乘檢驗統計量

(7)

1.2.6Score檢驗統計量

通過非限制性的對數似然函數，可得到關于向量π=(π1,…,πJ)′的一階偏導數，根據Score檢驗的一般理論(Rao[15])，可以給出檢驗H0:π1=π2=…=πJ的Score統計量為：

(8)

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

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

1.2.7似然比檢驗統計量

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

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

(9)

根據似然比檢驗理論可知，在原假設成立時，當Nj,nj→∞(j=1,…,J),Tl～χ2(J-1)。

1.3 漸近的檢驗過程

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

或者

tlwls≥z1-α

2 模擬研究

為了考察本文提出的檢驗統計量的統計性質，通過Monte Carlo模擬方法對基于以上檢驗統計量的各種檢驗過程進行評價。考慮了J=3的模型和如下的樣本量：平衡樣本量設計① 小樣本(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)。為了考察每個檢驗統計量犯第一類錯誤的概率，在以上設置的各種樣本量下，考慮參數設置：π=0.1,0.3,0.5，η=0.6,0.7,θ=0.75,0.9，即考慮了3(π的值)×2(η的值)×2(θ的值)=12種參數組合。

為了考察每個檢驗統計量的經驗功效，考慮參數設置：η=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種參數組合。對于每個樣本量和每種參數組合，隨機產生5 000組隨機樣本{(n11j,n10j,n01j,n00j,xj,yj):j=1,2,3}分別計算本文提出的7種檢驗統計量的值，其中，{(n11j,n10j,n01j,n00j):j=1,2,3}從多項分布(nj:πjη,πj(1-η),(1-θ)(1-πj),θ(1-πj))中產生，{(xj,yj):j=1,2,3}從二項分布(Nj-nj:pj,1-pj)中產生，進而估計每個檢驗統計量下犯第一類錯誤的概率和檢驗功效。其中，犯第一類錯誤的概率估計公式為：被檢驗統計量T拒絕的次數/5 000(H0成立下)，經驗功效的估計公式為：被檢驗統計量T拒絕的次數/5 000(H1成立下)。顯著性水平α=0.05，犯第一類錯誤的模擬結果如表2～5所示。

由于篇幅的限制，本文只列出了中等樣本(平衡設計和非平衡設計)下的經驗功效的模擬結果，如表6和表7所示。

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

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

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

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

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

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

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

3 實際數據分析

為了評價本文所提出的檢驗統計量的統計性質，本文利用Nedelman[16]中的瘧疾數據來說明本文提出方法的有效性。世界衛生組織和尼日利亞加爾基當地政府為了研究當地瘧疾的患病率，分別考慮了調查的方法、村莊類型、年齡、性別等不同因素對當地瘧疾患病率的研究。為了研究不同年齡對瘧疾患病率的影響，本文選取了3個年齡段(成年人數據)：即19～28歲、29～43歲、≥44歲的調查數據進行分析，數據如表8所示。

表8 尼日利亞加爾基瘧疾數據

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

4 結論

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

附錄：Score檢驗統計量的推導

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

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

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

在原假設H0成立之下的Score檢驗統計量為