王黎明,劉青松,張曉良,邱世芳
(重慶理工大學 理學院,重慶 400054)
研究某種疾病在人群中的流行率對生物醫學相關研究具有重要意義,一般情況下通過判斷受試者是否患有某種疾病對其進行分類從而獲得該疾病的患病率。然而,在現實生活中,金標準檢驗無誤判但價格昂貴,篩檢方法價格低但常有誤判。為了克服二者的不足,Tenenbein[1]提出了二重抽樣的方法,即從總體中隨機抽取N個個體利用篩檢方法進行檢驗,再從N個個體中隨機抽取n個個體接受金標準檢驗。這樣,這n個個體同時接受了2種檢驗,而N-n個個體只接受了篩檢檢驗,因此,Tang等[2]將采用二重抽樣方法獲得的數據稱為部分核實數據。
在部分核實數據的研究方面,已有大量成果,例如文獻[3-10],但是上述研究都是基于有金標準分類器存在的情況下進行的。實際生活中,金標準分類器并不一定存在或可用,而無金標準存在的部分核實數據在實際生活中廣泛存在。一個典型例子是Nedelman[11]利用二重抽樣數據研究了瘧疾的流行率,使用2種分類器(即初級顯微鏡者和高級顯微鏡者)對病人是否患有瘧疾進行分類。由于無論是初級還是高級顯微鏡都存在誤判,因而不存在金標準。表1列出了1個干旱季節和潮濕季節下2個年齡組(9~18歲,19~29歲)的瘧疾數據。感興趣的問題是:在不同的季節和不同的年齡組,瘧疾的流行率是否有顯著不同。為此,Qiu等[12]從齊性檢驗的角度進行了研究。然而,在醫學研究中,需要多少個體參與試驗也是一個重要的研究課題。Qiu等[13]從疾病流行率顯著性檢驗的角度對樣本量的確定進行了探究。Tang等[14]基于比例比,分別從假設檢驗與置信區間寬度的角度確定了近似樣本量公式。然而,在給定置信水平下,關于疾病流行率之差(比例差)的置信區間的寬度控制在指定范圍內的樣本量確定還沒有相關研究文獻。因此,本文中將從置信區間寬度的角度出發對此問題進行研究,提出幾種有效的樣本量的確定公式或有效算法。如Nedelman[11]所論述,對此類問題假定不存在假陽性誤判是合理的。因而,研究不存在假陽性誤判下基于流行率之差的區間寬度控制下的樣本量的確定問題。

表1 瘧疾數據
假設有來自第j組的Nj個個體,對每個個體先用初級分類器進行檢驗,Jj=1表示陽性,Jj=0表示陰性。接著從Nj個樣本中隨機選取nj個個體再用高級分類器進行檢驗,Sj=1表示陽性,Sj=0表示陰性(j=1,2),數據如表2所示。

表2 二重抽樣的數據類型
令Dj表示第j組中個體真實患病的情況,Dj=1表示個體患病,反之則沒有(j=1,2)。設πj=Pr(Dj=1)為患病率,ηj=Pr(Jj=1|Dj=1)和θj=Pr(Sj=1|Dj=1)分別表示初級和高級分類器的敏感度,即在患病條件下檢驗為陽性的條件概率。假設初級分類器和高級分類器滿足條件獨立性,即Pr(Jj,Sj|Dj)=Pr(Jj|Dj)Pr(Sj|Dj)。假設不存在假陽性誤判下,其概率結構(稱之為模型1)如表3所示。

表3 模型1的概率結構
設mj=(n11j,n10j,n01j,n00j,xj,yj),m={mj∶j=1,2},且π=(π1,π2),η=(η1,η2),θ=(θ1,θ2),則模型1下的似然函數為



其中C1=log(A1)。根據文獻[15],當n11jn00j≥n10jn01j時,參數πj,ηj,θj的極大似然估計為


在不假定條件獨立性下,Lie等[16]在假定Pr(Jj=1,Sj=0|Dj=1)=Pr(Sj=0|Dj=1)、Pr(Jj=0,Sj=1|Dj=1)=Pr(Jj=0|Dj=1)下,提出了模型2,如表4所示。

表4 模型2的概率結構
模型2的似然函數為:


其中C2是與參數無關的常數。根據文獻[15],可得參數πj,ηj,θj的極大似然估計分別為


假設核實驗證比例為pj,即nj=Njpj,(j=1,2),并在N2=N1r下考慮樣本量的確定問題。
基于以上得到的δ的估計和其估計的方差,易得到參數δ的基于Wald檢驗的置信水平為1-α的置信區間為:

其中,zα/2是標準正態分布的α/2分位數。為了將置信區間寬度控制在2ω以內,當ω≤δ時設否則由此可得置信區間寬度控制在指定寬度2ω所需的樣本量為:

對于模型1:

對于模型2:


則δ的置信水平為1-α的置信區間為:

通過計算得到,基于對數變換的置信區間寬度控制在2ω所需的樣本量為:

其中,模型1與模型2下,V分別由式(9)與(10)得到。

則δ的置信水平為1-α的置信區間為

通過計算可得基于Logit變換的置信區間寬度控制在2ω所需的樣本量為

其中,模型1與模型2下,V分別由式(9)與(10)得到。
由于基于Score和似然比檢驗統計量的δ的置信區間計算中需要參數π1,ηj,θj在δ=δ0下的限制性極大似然估計而這些估計沒有顯表達式,需要通過迭代算法得到,故考慮Zou等[17]提出的MOVER方法求得δ的置信區間。設δ的置信水平為1-α的置信區間上限δU和下限δL,根據中心極限定理易得:

同樣,通過中心極限定理可直接求出πj(j=1,2)的置信水平為1-α的雙邊置信區間的上限uj和下限lj分別為:

由此可得:

根據MOVER方法,下限δL的替換得到:

類似地,將上限δU的替換得到:

分別基于以下的Score檢驗統計量與似然比檢驗統計量求得πj(j=1,2)的置信水平為1-α的雙邊置信區間的上限uj和下限lj,代入上述δL和δU的式中,從而得到δ的置信水平為1-α的置信區間。
2.4.1 基于Score檢驗統計量的置信區間

上述方程組無顯式解,可通過迭代算法求得。
在模型2下:

基于Rao[18]的理論,在模型1下,檢驗假設H0j∶πj=π0j的Score檢驗統計量為

在模型2下:

根據文獻[19],當nj→∞時,在H0j下Tsc(π0j)漸近服從標準正態分布,利用迭代方法解如下方程可得到基于Score檢驗統計量的置信水平為1-α的πj置信區間上限與下限為

2.4.2 基于似然比檢驗統計量的置信區間
在模型1下,檢驗假設H0j∶πj=π0j的似然比檢驗統計量為

在模型2下,檢驗假設H0j∶πj=π0j的似然比檢驗統計量為

根據文獻[20],似然比檢驗統計量Tl(π0j)漸近服從自由度為1的卡方分布,通過解方程Tl(π0j)=可得到基于似然比檢驗統計量的置信水平為1-α的πj置信區間上下限。
2.4.3 樣本量的數值解法
由于基于以上檢驗統計量Tsc(π0j)和Tl(π0j)計算的πj置信區間均沒有顯表達式,采取如下的搜索算法來獲得近似樣本量的估計:
步驟1給定π1,ηj,θj,δ,r,pj,N1的值,產生K組隨機樣本mj=(n11j,n10j,n01j,n00j,xj,yj)。
令M(·)表示多項分布,在模型1下,(n11j,n10j,n01j,n00j)~M(nj;πjηjθj,πj(1-ηj)θj,πjηj(1-θj),πj(1-ηj)(1-θj)+(1-πj))。
在模型2下,(n11j,n10j,n01j,n00j)~M(pjNj;πj(ηj+θj-1),πj(1-ηj),πj(1-θj),(1-πj))。
兩種模型下(xj,yj)都服從二項分布B((1-pj)Nj,πjηj)。
步驟2基于步驟1產生的樣本,結合MOVER方法計算δ的置信區間。并通過K個置信區間寬度的平均值計算近似的區間寬度,記為2ω*(N1)。
步驟3若2ω*(N1)大于(小于)2ω,則增加(減少)N1的值,重復步驟1、2。
步驟4重復步驟3直到近似的半區間寬度ω*(N1)非常接近于給定的區間寬度ω,則N1=min{N1∶|ω*(N1)-ω|≤0.001}即為滿足條件的近似樣本量。
通過以上搜索方法可獲得基于Tsc(π0j)和Tl(π0j)的置信區間的樣本量Nsc和Nl。
對于模型1和模型2下所提出的樣本量的估計方法,通過計算在估計的樣本量下δ的置信區間的經驗覆蓋概率(ECP)和經驗覆蓋寬度(ECW)來評估本文中所提出方法的有效性。1)經驗覆蓋概率(ECP)

[δL(m(k)),δU(m(k))]為δ的第k個置信區間,:j=1,2},此時I{δ∈[δL(m(k)),δU(m(k))]}為示性函數。
2)經驗覆蓋寬度(ECW)

在置信水平1-α=0.95下分別通過以下參數設置考察了各種樣本量估計方法的有效性:
1)δ的影響。為研究δ的變化對樣本量的影響,在2種模型下,考慮如下參數設置:δ=0.1(0.01)0.3,π1=0.3,η1=θ1=0.8,η2=0.85,θ2=0.75,ω=0.1以及(a):r=1.0,p1=p2=1/3,(b):r=2.0,p1=p2=1/3,(c):r=1.0,p1=1/3,p2=1/5,(d):r=2.0,p1=1/3,p2=1/5。圖1、3給出了估計的樣本量N1、置信區間的經驗覆蓋概率ECP(%)和經驗覆蓋寬度ECW。對于模型1和模型2,隨著δ的逐漸增大,基于所有方法計算的所需樣本量均逐漸增加,基于對數變換置信區間的樣本量Nlog和基于Logit變換置信區間的樣本量Nlogit相對其他3種方法較小,所有方法的ECP均接近于事先給定的置信水平1-α=0.95,ECW接近于2ω=0.2。
2)p1和p2的影響。為研究p1的變化對樣本量的影響,考慮如下參數設置:p1=0.1(0.1)0.9,π1=0.3,δ=0.2,η1=θ1=0.8,η2=0.85,θ2=0.75,ω=0.1,p2=0.3,(a):r=1.0,(b):r=2.0。為研究p2的變化對樣本量的影響,考慮如下參數設置:p2=0.1(0.1)0.9,π1=0.3,δ=0.2,η1=θ1=0.8,η2=0.85,θ2=0.75,ω=0.1,p1=0.3,(c):r=1.0,(d):r=2.0。圖2、4給出了估計的樣本量N1、置信區間的經驗覆蓋概率ECP(%)和經驗覆蓋寬度ECW。模擬研究結果表明,隨著p1和p2的增大,基于所有方法計算得到的樣本量均減少,所有方法的ECP均接近于事先給定的置信水平1-α=0.95,ECW接近于2ω=0.2。
圖1中,δ=0.1(0.01)0.3,π1=0.3,η1=θ1=0.8,η2=0.85,θ2=0.75,ω=0.1。

圖1 模型1:置信水平1-α=0.95下,估計的樣本量N1、置信區間的經驗覆蓋概率ECP(%)和經驗覆蓋寬度ECW隨δ變化的情況
圖2中,π1=0.3,δ=0.2,η1=θ1=0.8,η2=0.85,θ2=0.75,ω=0.1。

圖2 模型1:置信水平1-α=0.95下,估計的樣本量N1、置信區間的經驗覆蓋概率ECP(%)和經驗覆蓋寬度ECW隨p1、p2變化的情況
圖3中,δ=0.1(0.01)0.3,π1=0.3,η1=θ1=0.8,η2=0.85,θ2=0.75,ω=0.1。

圖3 模型2:置信水平1-α=0.95下,估計的樣本量N1、置信區間的經驗覆蓋概率ECP(%)和經驗覆蓋寬度ECW隨δ變化的情況
圖4中,π1=0.3,δ=0.2,η1=θ1=0.8,η2=0.85,θ2=0.75,ω=0.1。

圖4 模型2:置信水平1-α=0.95下,估計的樣本量N1、置信區間的經驗覆蓋概率ECP(%)和經驗覆蓋寬度ECW隨p1、p2變化的情況
采用所提出的方法對表1數據進行分析。表5列出了2個調查組下不同年齡組患病率之差δ、2組年齡下不同調查組患病率之差δ的參數估計,置信水平1-α=0.95下區間寬度控制在區間寬度ω=0.1下的樣本量估計,在估計樣本量下δ的置信區間的經驗覆蓋概率和經驗覆蓋寬度。結果表明:對于模型1,基于對數變換方法的ECW相對于其他4種方法較小,在模型2下則相反。值得注意的是,在模型1與模型2下,基于Wald方法置信區間得到的樣本量明顯大于其他4種方法的樣本量,但所有方法的ECP均接近于事先給定的置信水平1-α=0.95,ECW接近于2ω=0.2,現有樣本量滿足實驗需求。
基于無金標準的部分核實數據,從比例差的區間寬度的角度給出了比較2組患病率是否有顯著性差異所需要的樣本量。在2種不同二重抽樣模型下,給出了在給定置信水平下置信區間寬度控制在指定范圍內的樣本量估計公式。模擬結果表明:相同參數設置下模型1下各種方法所需樣本量大于模型2下相應方法所需的樣本量;在2種模型下,基于所提出的各種方法估計的樣本量下計算的δ的置信區間的經驗覆蓋概率接近給定的置信水平且經驗覆蓋寬度接近設定的寬度。由此可知,所提出的樣本量估計公式和數值算法有效,推薦在實際應用中使用。
附錄:Fisher信息陣
令Ij=(Ilkj)3×3表示第j組下2個模型的費希爾信息矩陣,其中,在模型1下,

在模型2下
