邱世芳，曾小松

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

1 背景

在生物醫(yī)學(xué)和流行病學(xué)研究中，為了分析一種疾病的流行率，人們通常使用價格便宜的篩檢方法對感興趣的個體是否感染這種疾病進行第1次臨床診斷。但是，這種篩檢方法經(jīng)常導(dǎo)致因不可忽略的診斷誤差而出現(xiàn)存在誤判的分類數(shù)據(jù)；另一方面，雖然可以選擇一種完全準(zhǔn)確的金標(biāo)準(zhǔn)檢驗，但是這種檢驗具有價格昂貴、耗時長而且對檢驗個體有副作用等不足。為了克服這兩種方法的不足，常采用二重抽樣方法獲得部分核實數(shù)據(jù)，即從感興趣的總體中隨機抽取N個個體，對每個個體都使用篩檢方法進行診斷檢測，然后再從中隨機抽取n(n

在完全無誤判的金標(biāo)準(zhǔn)存在時，已有文獻對基于部分核實數(shù)據(jù)的統(tǒng)計推斷問題進行研究。例如：Boese[4]對二重抽樣設(shè)計下有假陽性錯誤分類時，提出了關(guān)于二項分布參數(shù)的幾種基于似然的近似區(qū)間估計；Morvan[5]對篩檢過程的準(zhǔn)確性以及兩階段調(diào)查中疾病流行率的修正估計的準(zhǔn)確性提出了兩種評價方法；對于單個樣本情形，Tang等[3]基于疾病的流行率提出了有實際意義的統(tǒng)計假設(shè)，并對此假設(shè)考慮了Score檢驗統(tǒng)計量、似然比檢驗統(tǒng)計量以及兩種Wald型檢驗統(tǒng)計量，從大樣本的角度考慮了漸近的檢驗過程以及從小樣本的角度考慮了近似非條件的檢驗過程；Tang等[6]基于上述檢驗統(tǒng)計量提出了關(guān)于疾病流行率的12種區(qū)間估計方法；Qiu等[7]分別從顯著性檢驗的角度和置信區(qū)間的角度研究了疾病流行率顯著性檢驗的樣本量的計算公式。對于兩組獨立樣本情形，Tang等[8]基于疾病流行率之差考慮了兩組獨立樣本下疾病流行率的顯著性差異的假設(shè)檢驗問題，分析了多種檢驗統(tǒng)計量(如Score或者Score型檢驗)的大樣本的漸近檢驗過程和小樣本的近似非條件的檢驗過程，并從檢驗功效的角度研究了樣本量的確定問題； Qiu等[9]基于比例差對不同組的流行病治愈率的等價性評價進行了研究，從置信區(qū)間的角度提出了等價性檢驗中基于比例差的置信區(qū)間的估計方法。

以上內(nèi)容都是基于有完全無誤判金標(biāo)準(zhǔn)存在時部分核實數(shù)據(jù)的研究。事實上，完全無誤判的金標(biāo)準(zhǔn)通常是不存在的，接受核實驗證的個體通常接受的是不完全無誤判金標(biāo)準(zhǔn)檢驗。當(dāng)這種不完全金標(biāo)準(zhǔn)檢驗比篩檢有更高的敏感度和特異度時，這樣的二重抽樣設(shè)計仍然可以降低統(tǒng)計推斷的偏差，因而基于不完全無誤判金標(biāo)準(zhǔn)下的核實數(shù)據(jù)的統(tǒng)計推斷更具有現(xiàn)實意義和實際研究價值。一個典型的例子是Nedelman[10]利用二重抽樣方法對瘧疾病的流行率進行了研究，根據(jù)年齡的不同將總體分為了7個年齡組。首先初級的顯微鏡醫(yī)生對每一個年齡組的個體進行檢驗，然后從中隨機抽取一部分個體再接受高級顯微鏡醫(yī)生的檢驗，從而得到部分核實數(shù)據(jù)。在此研究中，高級顯微鏡醫(yī)生雖然比初級醫(yī)生有更高的敏感度，但是他們的檢驗仍然有誤判可能性?；谶@樣的部分核實數(shù)據(jù)，Qiu[11]在沒有假陽性假定下對疾病流行率的區(qū)間估計進行了研究，考慮了兩種模型下的多種有效的區(qū)間估計方法。但是，實驗樣本量的確定問題是實際醫(yī)務(wù)工作者最為關(guān)心的問題之一，而從區(qū)間估計的角度對該問題的研究還沒有文獻涉及，因此本文對不完全無誤判金標(biāo)準(zhǔn)下的二重抽樣設(shè)計中所需要的樣本量進行研究。

2 模型與參數(shù)估計

假設(shè)從感興趣的總體中隨機抽取N個個體并用篩檢方法進行檢驗，檢驗結(jié)果用變量J表示，J=1表示陽性，J=0表示陰性。然后從抽取的N個個體中再隨機抽取n個個體接受不完全無誤判金標(biāo)準(zhǔn)檢驗，檢驗結(jié)果用變量S表示，S=1表示陽性，反之則表示陰性。二重抽樣的數(shù)據(jù)結(jié)構(gòu)如表1所示。

表1 二重抽樣的數(shù)據(jù)結(jié)構(gòu)

假設(shè)D表示個體患病狀態(tài)變量，D=1表示患有疾病，反之表示正常。令π=Pr(D=1)，η=Pr(J=1|D=1)和θ=Pr(S=1|D=1)，即η和θ分別表示篩檢和不完全金標(biāo)準(zhǔn)檢驗的敏感度。本文考慮不存在假陽性誤判的情形，因而它們的特異度皆為1。

2.1 模型1及其參數(shù)估計

假定兩種檢驗方法滿足條件獨立性，即Pr(J,S|D)=Pr(J|D)Pr(S|D)時，Nedelman[10]提出了一種模型(記為模型1)，其概率結(jié)構(gòu)如表2所示。

表2 模型1的概率結(jié)構(gòu)

令m=(n11,n10,n01,n00,x,y),則在此概率模型下m=(n11,n10,n01,n00,x,y)的似然函數(shù)為：

L1(m;π,η,θ)=A1(πηθ)n11(π(1-η)θ)n10(πη(1-θ))n01(1-π(η+θ-ηθ))n00(πη)x(1-πη)y=

A1πn-n00+xηn+1+x(1-η)n10θn1+(1-θ)n01(1-π(η+θ-ηθ))n00(1-πη)y

其中A1是與參數(shù)無關(guān)的常數(shù)，其對數(shù)似然函數(shù)為

l1(m;π,η,θ)=C1+(n-n00+x)logπ+(n+1+x)logη+n10log(1-η)+n1+logθ+

n01log(1-θ)+n00log(1-π(η+θ-ηθ))+ylog(1-πη)

(1)

其中C1=logA1為常數(shù)。

令λ1=Pr(S=1|J=1)，λ2=Pr(S=1|J=0)和p=Pr(J=1)，則λ1=θ，λ2=π(1-η)θ/(1-πη)，p=πη，因此對數(shù)似然函數(shù)(1)變?yōu)椋?/p>

l1(m;λ1,λ2,p)=C1+n11logλ1+n01log(1-λ1)+

n10logλ2+n00log(1-λ2)+(n+1+x)logp+(n+0+y)log(1-p)

則λ1、λ2和p的極大似然估計為：

由此易得到，當(dāng)n11n00≥n10n01時，π、η和θ的極大似然估計分別為：

(2)

(3)

2.2 模型2及其參數(shù)估計

假定兩種分類器不具有條件獨立性時， Lie等[12]提出了一種模型(記為模型2)，其概率模型結(jié)構(gòu)如表3所示。

表3 模型2的概率結(jié)構(gòu)

在此概率模型下m=(n11,n10,n01,n00,x,y)的似然函數(shù)為

L2(m;π,η,θ)=A2(π(η-1+θ))n11(π-πη)n10(π-πθ)n01(1-π)n00(πη)x(1-πη)y=

A2πn-n00+x(1-π)n00ηx(1-η)n10(1-θ)n01(η+θ-1)n11(1-πη)y

其對數(shù)似然函數(shù)為：

l2(m;π,η,θ)=C2+(n-n00+x)logπ+n00log(1-π)+xlogη+n10log(1-η)+

n01log(1-θ)+n11log(η+θ-1)+ylog(1-πη)

(4)

其中A2和C2為與參數(shù)無關(guān)的常數(shù)。

在此模型下，由λ1、λ2和p的定義可得λ1=(η+θ-1)/η，λ2=π(1-η)/(1-πη)和p=πη。采用和模型1類似的方法可得π、η和θ的極大似然估計分別為：

(5)

(6)

3 樣本量的確定

本文假設(shè)核實驗證比例為κ，即在κ=n/N情況下考慮實驗樣本量的確定問題。

3.1 基于Wald置信區(qū)間的樣本量

其中zα/2為標(biāo)準(zhǔn)正態(tài)分布的上α/2分位數(shù)。為了將置信區(qū)間寬度控制在2ω以內(nèi)，當(dāng)ω≤π時設(shè)zα/2σ=ω，否則π+zα/2σ=2ω。由此得到基于Wald檢驗的100(1-α)%置信區(qū)間寬度控制在指定寬度2ω所需的樣本量(即N)，其計算公式為：

(7)

其中:對于模型1有

對于模型2有

由此容易發(fā)現(xiàn)：在模型2下的樣本量與不完全無誤判金標(biāo)準(zhǔn)的敏感度無關(guān)。

3.2 基于限制性極大似然估計下方差的Wald置信區(qū)間、似然比置信區(qū)間和Score置信區(qū)間的樣本量

3.2.1 基于限制性極大似然估計下方差的Wald置信區(qū)間

(8)

此方程組沒有顯示解，因此可通過一種迭代方法(如牛頓迭代法)獲得。

其中：a=π0(N-n00)；b=N-y-n00+π0(N-n+0)；c=b2-4a(x+n+1)。

在模型2下，有

3.2.2 基于似然比檢驗的置信區(qū)間

兩種模型下的似然比檢驗統(tǒng)計量為

3.2.3 基于Score檢驗的置信區(qū)間

對于模型1，Score函數(shù)為：

對于模型2， Score函數(shù)為：

則對于假設(shè)檢驗H0:π=π0的Score檢驗統(tǒng)計量為

或

其中I11的表達式見附錄。在原假設(shè)下Tsc漸近服從標(biāo)準(zhǔn)正態(tài)分布，因而基于Score統(tǒng)計量的(1-α)100%的置信區(qū)間[πl(wèi),πu]的上下限可利用迭代方法通過以下方程得到：

Tsc=±z1-α/2

其中πl(wèi),πu分別對應(yīng)方程中取“+”和“-”得到。

3.2.4 樣本量的求解算法

由于基于以上檢驗統(tǒng)計量Tw2、Tl和Tsc的置信區(qū)間都沒有顯表達式，因此考慮以下搜索算法獲得近似樣本量的估計:

步驟1 給定π,η，θ，N和κ的值，產(chǎn)生K組隨機樣本m=(n11,n10,n01,n00,x,y)，其中：在模型1中,(n11,n10,n01,n00)～M(Nκ;πηθ,π(1-η)θ,πη(1-θ),π(1-η)(1-θ)+1-π)；在模型2中,(n11,n10,n01,n00)～M(Nκ;π(η+θ-1),π(1-η),π(1-θ),1-π)。在兩種模型下(x,y)都服從二項分布B(N(1-κ),πη)，其中M(·)表示多項分布。

步驟2 基于步驟1產(chǎn)生的每一個樣本，計算π的(1-α)100%置信區(qū)間，并通過K個置信區(qū)間寬度的平均值計算近似的區(qū)間寬度，記為2ω*(N)。

步驟3 若2ω*(N)小于(大于)2ω，那么減少(增加)N的值，重復(fù)步驟1和步驟2。

步驟4 重復(fù)步驟3直到近似的半?yún)^(qū)間寬度ω*(N)非常接近于給定的區(qū)間寬度ω，則N=minN:|ω*(N)-ω|≤0.001}，即為滿足條件的近似樣本量。

通過以上搜索方法獲得的基于Tw2、Tl和Tsc的置信區(qū)間的樣本量公式分別記為Nw2、Nl和Nsc。

4 模擬研究

為了評價本文所提出方法的效果，考慮將置信水平1-α=0.95下的置信區(qū)間寬度的一半控制在ω=0.05下的樣本量的確定，并在估計的樣本量(即Nw1，Nw2，Nl和Nsc)下，通過模擬分別研究各置信區(qū)間的經(jīng)驗覆蓋概率和經(jīng)驗覆蓋寬度。假設(shè)核實比例即接受不完全金標(biāo)準(zhǔn)檢驗的個體比例κ=0.2，參數(shù)設(shè)置為π=0.2、0.5和0.8，由于不完全無誤判金標(biāo)準(zhǔn)檢驗不劣于篩檢方法，因而對于敏感度設(shè)置為η=0.5和0.7，θ=0.75和0.95，總共3(π的取值)×2(η的取值)×2(θ的取值)=12種參數(shù)組合進行蒙特卡洛模擬。

在給定參數(shù)組合下，估計出樣本量Nw1、Nw2、Nl和Nsc，在估計的樣本量下，產(chǎn)生2 000個隨機樣本m=(n11,n10,n01,n00,x,y)，基于每個統(tǒng)計量計算相應(yīng)的置信區(qū)間，從而得到置信區(qū)間的經(jīng)驗覆蓋概率(ECP)和經(jīng)驗覆蓋寬度(ECW)。模型1的模擬結(jié)果見表4，模型2的模擬結(jié)果見表5。

πηθNl(ECP，ECW)Nsc(ECP，ECW)0．20．50．751435(95．1，0．101)1495(94．2，0．101)0．95944(95．1，0．100)954(95．6，0．099)0．70．75889(94．7，0．100)949(94．6，0．100)0．95654(95．7，0．101)674(94．4，0．100)0．50．50．752710(94．3，0．101)2710(94．1，0．101)0．951611(95．7，0．100)1611(94．5，0．100)0．70．751653(94．9，0．100)1653(95．2，0．101)0．951180(94．7，0．099)1180(95．5，0．099)0．80．50．753115(94．0，0．100)3115(93．7，0．100)0．951389(95．0，0．100)1389(95．5，0．101)0．70．751795(94．9，0．100)1795(94．8，0．100)0．951056(95．2，0．099)1056(95．2，0．098)

表5 基于模型2的95%置信區(qū)間寬度控制在2ω=0.1下的近似樣本量、經(jīng)驗覆蓋概率(%)和經(jīng)驗區(qū)間寬度

πηθNl(ECP，ECW)Nsc(ECP，ECW)0．20．50．75799(95．4，0．099)800(95．3，0．099)0．95799(94．0，0．099)800(95．1，0．098)0．70．75600(95．6，0．099)626(94．9，0．098)0．95609(96．1，0．099)620(95．3，0．098)0．50．50．751370(94．7，0．100)1360(95．2，0．100)0．951370(94．9，0．100)1360(94．4，0．101)0．70．751065(94．7，0．100)1057(94．9，0．101)0．951065(94．8，0．100)1060(95．2，0．101)0．80．50．751026(94．9，0．101)1020(94．8，0．101)0．951025(95．2，0．101)1024(95．7，0．100)0．70．75875(95．1，0．101)868(95．1，0．101)0．95879(94．6，0．101)870(94．8，0．101)

模擬研究表明：在基于限制性極大似然估計下方差的Wald置信區(qū)間、似然比置信區(qū)間和Score置信區(qū)間確定的樣本量下，它們的經(jīng)驗覆蓋概率接近事先給定的置信水平，且半?yún)^(qū)間寬度很好地控制在ω以內(nèi)，因而基于這3種方法確定的樣本量準(zhǔn)確有效。但是，在基于非限制性極大似然估計下方差的Wald置信區(qū)間確定的樣本量下，其經(jīng)驗覆蓋概率不能保證其接近給定的覆蓋概率。其次，隨著患病概率的增大，要使95%的置信區(qū)間寬度控制在2ω以內(nèi)，所需要的樣本量就越大，隨著兩種檢驗的敏感度的增大，所需的樣本量相應(yīng)減少。通過對比4種方法可以看到：模型1基于非限制性極大似然估計下方差的Wald置信區(qū)間、似然比置信區(qū)間和Score置信區(qū)間確定的樣本量很接近，沒有顯著差別,只有當(dāng)π比較小的時候相差較大,限制性極大似然估計下方差的Wald置信區(qū)間與其他3種方法相差較大；而模型2基于限制性極大似然估計下方差的Wald置信區(qū)間、似然比置信區(qū)間和Score置信區(qū)間確定的樣本量相差不大，沒有顯著差異，基于非限制性極大似然估計下方差的Wald置信區(qū)間與這3種方法所計算的樣本量大概相差40。最后，在相同的參數(shù)設(shè)置下，基于模型2所需要的樣本量通常比模型1要小。因此，當(dāng)兩種檢驗不滿足條件獨立性時，應(yīng)考慮模型2的統(tǒng)計分析。

5 實例分析

首先，對于引言中的瘧疾數(shù)據(jù)，采用本文提出的方法進行研究。由于篇幅的限制，僅列出小于1歲的年齡組的分析結(jié)果，其他年齡組的結(jié)果不再一一列出。此年齡組的觀測數(shù)據(jù)為n11=165、n10=8、n01=7、n00=80、x=1736、y=799，即n=260、N=2795。通過數(shù)據(jù)易知：核實數(shù)據(jù)的比例約等于0.09，在模型1中參數(shù)π、η和θ的極大似然估計分別為0.712 7、0.957 8和0.959 3，在模型2中參數(shù)的極大似然估計分別為0.711 5、0.959 5和0.961 0。因此，在樣本量的估計中，令κ=0.1，π=0.70,η=0.95,θ=0.95，置信水平為95%，半?yún)^(qū)間寬度ω=0.05，利用本文提出的樣本量估計方法，分別得到兩種模型下的樣本量的估計值，并計算在估計的樣本量下95%置信區(qū)間的經(jīng)驗覆蓋概率和經(jīng)驗覆蓋寬度。在模型1下，Nw1=786，Nw2=846，Nl=816，Nsc=856，其經(jīng)驗覆蓋概率(覆蓋寬度)分別為92.3%(0.097 7),96.1%(0.100 4),95.6%(0.100 3)和95.1%(0.100 7);在模型2下，Nw1=757，Nw2=800，Nl=799，Nsc=813，其經(jīng)驗覆蓋概率(覆蓋寬度)分別為93.5%(0.097 6),95.1%(0.099 8),95.4%(0.098 7)和95.3%(0.101 4)。由此可見，要使95%的置信區(qū)間寬度控制在0.10內(nèi)，現(xiàn)有的樣本量是滿足的。

其次，對Qiu等[7]使用的移植時年齡小于20歲的再生障礙性貧血患者數(shù)據(jù)進行分析，原始數(shù)據(jù)為：n11=6、n10=4、n01=3、n00=18、x=12、y=49，即n=31、N=92。在k=0.3，π=0.4，η=0.55，θ=0.75且置信水平為95%、半?yún)^(qū)間寬度ω=0.05下，基于4種方法在模型1下的樣本量分別為Nw1=1 424，Nw2=1 534，Nl=1 424，Nsc=1 454，模型2下樣本量分別為Nw1=866，Nw2=837，Nl=841和Nsc=840。

在有金標(biāo)準(zhǔn)的情況下(Qiu等[7])結(jié)果為：Nw1=650，Nw2=656，Nl=657和Nsc=665。通過比較可以看到：有金標(biāo)準(zhǔn)所需的樣本量比沒有金標(biāo)準(zhǔn)所需樣本量小，并且無論在哪種情況下現(xiàn)有的樣本量都不能使95%的置信區(qū)間寬度控制在0.10內(nèi)。

6 討論

本文在不完全無誤判金標(biāo)準(zhǔn)下，基于部分核實數(shù)據(jù)從置信區(qū)間的角度研究檢驗所需的樣本量，分別考慮了兩種Wald檢驗、似然比檢驗和Score檢驗4種方法在兩種模型下樣本量的確定情況，然后對這些方法進行模擬研究，通過經(jīng)驗覆蓋率和區(qū)間寬度評價本文所提出方法的準(zhǔn)確性。通過模擬結(jié)果發(fā)現(xiàn)：在相同的參數(shù)設(shè)置下，基于模型2所需樣本量比模型1小，且樣本量的大小與不完全無誤判金標(biāo)準(zhǔn)檢驗的敏感度θ無關(guān)?；谙拗菩詷O大似然估計方差的Wald檢驗、似然比檢驗和Score檢驗的樣本量公式具有良好的統(tǒng)計性質(zhì)，在估計的樣本量下的置信區(qū)間的經(jīng)驗覆蓋概率和經(jīng)驗覆蓋寬度都很接近給定的置信水平和區(qū)間寬度，因而可推薦用于實際應(yīng)用中。

模型1假設(shè)兩種檢驗滿足條件獨立性，這個條件在實際中不一定滿足，所以在實際應(yīng)用時要先判斷用于診斷的檢驗方法是否滿足條件獨立性。模型2的假設(shè)在實際中同樣可能不成立，所以應(yīng)根據(jù)具體實際情況來選擇相應(yīng)的模型。通過對實際數(shù)據(jù)的分析可以看到：本文研究成果具有實際應(yīng)用價值，可為醫(yī)務(wù)工作者提供指導(dǎo)。本文從經(jīng)典的頻率統(tǒng)計方面對樣本量的問題進行研究，在未來的研究中還可考慮采用貝葉斯等方法。

[1] TENENBEIN A A.A double sampling scheme for estimating from binomial data with misclassifications[J].Journal of the American Statistical Association,1970,65:1350-1361.

[2] YIU C F,POON W Y.Estimating the polychoric correlation from misclassified data[J].British Journal of Mathematical and Statistical Psychology,2008,61:133-161.

[3] TANG M L,QIU S F,POON W Y,et al.Test procedures for disease prevalence with partially validated data[J].Journal of Biopharmaceutical Statistics,2012,22:368-386.

[4] BOESE D H,YOUNG D M,STAMEY J D.Confidence intervals for a binomial parameter based on binary data subject to false-positive misclassification[J].Computational Statatistics and Data Analysis,2006,50:3369-3385.

[5] MORVAN J,COSTE J,ROUX C H,et al.Guillemin F.Prevalence in two-phase surveys:accuracy of screening procedure and corrected estimates[J].Annuals of Epidemiology,2008,18:261-269.

[6] TANG M L,QIU S F,POON W Y.Confidence interval construction for disease prevalence based on partial validation series[J].Computational Statistics and Data Analysis,2012,56:1200-1220.

[7] QIU S F,POON W Y,TANG M L.Sample size determination for disease prevalence studies with partially validated data[J].Statistical Methods in Medical Research,2016,25(1):37-63.

[8] TANG M L,QIU S F,POON W Y.Comparison of disease prevalence in two populations in presence of misclassification[J].Biometrical Journal,2012,54(6):786-807.

[9] QIU S F,POON W Y,TANG M L.Confidence intervals for proportion difference from two independent partially validated series[J].Statistical Methods in Medical Research,2016,25(5):2250-2273.

[10] NEDELMAN J.The prevalence of malaria in Garki,Nideria:double sampling with a fallible expert[J].Biometrics，1988,44(3):635-655.

[11] QIU S F,LIAN H,ZOU G Y,et al.Interval estimation for a proportion using a double sampling scheme with two fallible classifiers[J].Statistical Methods in Medical Research,2016.DOI:10.1177/0962280216681599.

[12] LIU R T,HEUCH I，IRGENS L M.Maximum likelihood estimation of the proportion of congenital malformation using double registration systems[J].Biometrics,1994,50:433-444.

附錄

Score檢驗統(tǒng)計量的有關(guān)推導(dǎo)：

在模型1下， Fisher信息矩陣中各元素分別為：

Fisher信息矩陣的逆矩陣中第一個對角元素為：

在模型2下，F(xiàn)isher信息矩陣中各元素分別為：

Fisher信息矩陣的逆矩陣中第1個對角元素為：