刪失指標隨機缺失下條件分位數的加權核估計

周楊程, 王江峰, 袁汶汶, 張惠利

(浙江工商大學統(tǒng)計與數學學院,浙江杭州310018)

§1 引 言

在統(tǒng)計學中,對條件分位數的估計長期以來備受關注.在描述數據時,應用條件分位數會得到一些很好的性質,例如在處理重尾或離異的數據時,可以更加穩(wěn)健,尤其是條件中位數.令T表示因變量,其分布函數為F(·),X=(X1,···,Xd)>∈Rd為協(xié)變量,其聯合密度函數為v(·).設F(·,·)和f(·,·)分別為(X,T)的聯合分布函數和密度函數,則在X=x條件下，T的條件分布函數F(t|x)可以寫成：

通過(1)式得出條件分位數Qτ(x):

在過去的二三十年中,國內外有大量的文獻研究了條件分位數的估計問題.在獨立完整數據下,Chaudhuri[1]采用了多種非參數方法來估計條件分位數函數;Xiang[2]先構造了條件分布函數的核估計,由此得到了條件分位數的核估計.但在醫(yī)學和工程壽命等領域中,因變量T往往會由于許多原因被右刪失而不能完全觀測到.在右刪失數據下,Dabrowska[3]得到了分位數核估計量的Bahadur型表達式;Xiang[4]對隨機右刪失數據在等覆蓋概率準則下,比較了兩個核分位數估計與樣本分位數估計的性能;Ould-Sa¨?d[5]建立了條件分位數核估計的強一致收斂速度.

在以往研究右刪失數據時,刪失指標通常是可以完全觀測到的,但在很多實際應用中,由于存在種種復雜因素,部分刪失指標會隨機缺失.例如,在醫(yī)學領域中,當死亡原因調查表丟失,死亡原因未被記錄或由于病情復雜無法確定等情況出現時[6],研究某個特定死亡原因導致的死亡時間,那么刪失指標就會存在部分的缺失.當遇見刪失指標缺失的樣本時,最直接最簡單的方法將刪失指標缺失的數據全部去除,并用刪失數據下的統(tǒng)計方法對剩下的完全數據進行推斷,這種方法稱為完全數據(CC)方法.然而,這種CC方法推斷出來的結果常常是不可靠的,特別是當缺失嚴重時,刪除后的數據量也會大大減少,無法很好的反映總體.對于此類數據的處理,首先應當將其缺失機制進行劃分,Little和Rubin[7]提出可以分為三類:完全隨機缺失(MCAR),隨機缺失(MAR)以及不可忽略的缺失(non-ignorable missing).過去有大量文獻關于相對較強條件的MCAR缺失機制做了深入的研究,見Zhou和Sun[8]以及McKeague[9].

本文討論的是更一般的MAR情況,這種情況在許多實際問題中常常碰到,具體實例見Little和Rubin[7].最近,有許多作者對MAR情況做了不少工作.比如,Li和Wang[10]在研究線性回歸模型中提出了加權最小二乘法估計;Brunel等人[11]研究了條件風險率函數的估計問題,并對該估計提出了基于模型選擇方法的非參數自適應策略;Shen和Liang[12]應用分位數回歸方法研究了部分線性變系數模型.然而,在MAR情況下,目前還沒人研究條件分位數的估計問題,更別說是協(xié)變量為多維的情形.在這里,先分別構造條件分布函數的校準加權核估計,插值加權核估計以及逆概率加權核估計;然后根據這些估計分別導出條件分位數的核估計,并建立這些估計的漸近正態(tài)性結果;最后在有限樣本下,對這些估計進行模擬研究.

§2 估計的構造

設{(Xi,Ti),i=1,2,...,n}是一列來自總體(X,T)∈Rd+1的隨機向量序列,在右刪失情況下,由于刪失變量Ci的出現,Ti不能全部觀測到,因此實際觀測到的數據為Yi=min(Ti,Ci)以及δi=I(Ti≤Ci).設Ci的分布函數為G(·),與Xi獨立,并且在條件Xi下也與Ti獨立.

在右刪失數據下,Ould-Sa¨?d[5]給出在一維協(xié)變量下條件分布函數F(t|x)的核估計

Y(1)≤Y(2)≤···≤Y(n)為Y1,Y2,···,Yn的次序統(tǒng)計量,δ(i)為Y(i)相對應的量.

事實上,易知上式(3)可以看作下式的最小解:

本文考慮的是MAR情形,定義缺失指標ξi:當δi可觀測時,ξi=1;反之,ξi=0.假設缺失指標ξi和(Yi,Xi)相關,而與δi獨立,即P(ξi=1|Yi,Xi,δi)=P(ξi=1|Yi,Xi).最后實際觀察到的樣本數據為{Yi,Xi,δiξi,ξi,1≤i≤n}.由于刪失指標δi的隨機缺失,而(3)式和n(·) 都涉及到δi,因此上面的核估計不能直接用.過去在處理δi的隨機缺失時,可以根據Li和Wang[10]的方法,引入條件期望函數,即令φ(Z)=E(δ|Z),這里Z=(X,Y).

設Λ(·)為某個分布函數,λ(·)為它的密度函數,f(·|x)為在X=x條件下T的條件密度函數,Kd(·)定義在Rd上的核函數,.這樣有

因此通過(4)式,在MAR情形下,可以構造F(t|x)的估計為下式最優(yōu)解:

對(5)式不難求解,得到F(t|x)的估計為

在實際應用時,(6)式中的φ(Z)和G(·)往往是未知的.先處理φ(Z)的估計問題.注意到δ二進制的，對于二進制數據的建?？梢杂靡粋€指定模型形式來建立,例如Logit模型等.因此假設φ(Z)=φ0(Z;θ0),其中φ0(·;θ0)是已知的函數,θ0∈Rl是未知的參數向量.對于θ0的估計,這里可以沿用Li和Wang[10]的極大似然估計(MLE)方法,即下式的最大解：

故在MAR情形下,在(6)式中分別用φ0(Zi;n)和Gn(Yi)替換φ(Zi)和G(Yi)},得到F(t|x)的校準(CA)加權核估計為:

§3 主要結果

在給出結果之前,先給出一些注記.對于任何分布函數L(·),記

接下來,給出一些條件.

(A1)(i)H(t)在t

(ii)F1(r,s)和條件分布函數F(s|r)在任意的(r,s)∈U(x)×U(t)上具有連續(xù)的二階偏導數.

(A2)Kd(·)為定義在Rd上緊支撐的有界核函數,滿足

(A3)Λ(·)為一個連續(xù)且嚴格增的分布函數,其密度函數λ(·)滿足Rtλ(t)dt=0.

(A4)X的密度函數v(·)在點x上連續(xù),且v(x)>0.

(A5)窗寬hn滿足

(A6)條件密度函數f(s|r)在任意的(r,s)∈U(x)×U(t)上具有連續(xù)的二階偏導,且f(t|x)>0.

(A7)Σi,1(x,t)>0,i=1,2,3,且關于點x上連續(xù).

(A8)I(θ0)為正定矩陣,?φ0(Z;θ)在θ0上連續(xù);L(·)和?(·)是R上緊支撐的核函數,窗寬滿足

注3.1條件(A1)中的aH

§4 模擬研究

在這一節(jié),通過兩個例子(d=1,d=2),模擬研究Qτ(x)的各個估計的性質,內容主要包括:(i)比較Qτ(x)的CC(第一節(jié)提到過),CA,IW以及IPW核估計之間的整體均方誤差(GMSE)以及固定點上的偏移(Bias)和標準差(Sd);(ii)應用正態(tài)Q-Q圖研究各個估計的漸近正態(tài)性效果.

例4.1考慮協(xié)變量為一維的模型

這里的Xi～N(0,1),且與εi相互獨立.研究τ=0.5下Q0.5(x)的各估計的模擬情況.顯

然Q0.5(x)=2+sin(2x)+2exp{?16x}.為了計算Q0.5(x)的各個估計,分下面幾個步驟:

Step 1從模型(11)中,產生一組樣本數據{(Xi,Ti),1≤i≤n}.

Step 2獲取不同的刪失比例(CR):假設Ci～exp(μ)+2,Yi=min(Ti,Ci)以及δi=I(Ti≤Ci),這里μ用來調整獲取不同的CR.

Step 3獲取不同的缺失比例(MR):由于考慮的是MAR情況,假設Logit(E(ξi|Xi,Yi))=α1+α2Xi+α3Yi,這里通過調整(α1,α2,α3)的值獲得不同的MR.

Step 4φ0(Zi;n)和n(Yi)的計算:通過(11)以及步驟1-3,可以模擬一組觀察數據{(Yi,Xi,δiξi,ξi),1≤i≤n}.假設φ0(Zi;θ)服從Logit模型, 即Logit(φ0(Zi;θ))=θ1+θ2Xi+θ3Yi.參數θ=(θ1,θ2,θ3) 的估計根據§2提到的MLE方法得到; 對于n(Yi),選取L(x)=,然后由§2的乘積極限估計獲得.

Step 5計算Q0.5(x)的CA,IW以及IPW估計:選取,Λ(x)為標準正態(tài)分布函數,把φ0(Zi;n)和n(Yi)帶入(7)-(9)式,算出,i=1,2,3;然后根據,計算出Q0.5(x)的CA,IW以及IPW估計.

在比較Q0.5(x)的各估計好壞時,要用到GMSE,其定義為

在表格1和表格2中,分別對εi～N(0,0.52)和εi～t(3)兩種情況進行模擬.考慮CR=10%和40%,MR=10%和40%,n=200和500,基于M=200重復得到關于Q0.5(x)的CA,IW以及IPW估計的GMSE.窗寬hn的取值范圍從0.01到1.0,增量為0.02,分別選擇使GMSE達到最小的窗寬.另外,在x0=0.5上,基于200次重復計算各估計下的偏差和標準差.

表1 在例4.1和εi～N(0,0.52)中,CC,CA,IW以及IPW估計下的GMSE,Bias以及Sd的模擬結果

通過表1和表2,可以得到一些結論:(1)當樣本量n增大時,無論從GMSE,Bias還是sd的結果,CC,CA,IW以及IPW估計模擬效果都更好;而當CR和MR增大時,這四種估計模擬結果會變差;(2)在其他條件相同下,從GMSE以及Sd的結果來看,CA估計都最優(yōu),與定理3.1的結論相一致;但是從Bias的模擬結果,發(fā)現IW估計最優(yōu);(3)從CR=10%,MR=40%和CR=40%,MR=10%這兩組模擬的結果發(fā)現,CR的影響比MR的影響似乎更大;(4)當缺失比例很大時,即MR=40%,CC方法比CA,IW以及IPW估計方法模擬出來的結果都要差.

表2 在例4.1和εi～t(3)中,CC,CA,IW以及IPW 估計下的GMSE,Bias以及Sd的模擬結果

接下來,在εi～N(0,0.52)下,通過Q-Q正態(tài)圖來研究CA,IW以及IPW估計的漸近正態(tài)性效果.在圖1-圖3中,對n=500和M=500,分別在(CR=10%,MR=40%),(CR=40%,MR=10%)以及(CR=40%,MR=40%)下對Q0.5(x0)(x0=0.5)的CA,IW以及IPW估計作了QQ正態(tài)圖.從圖1-圖3看,發(fā)現CA估計的漸近正態(tài)性效果似乎更好,這和定理3.1中的結論(CA估計的漸近方差最小)一致;另外,可以發(fā)現隨著CR和MR的增大,各估計的漸近正態(tài)效果越差;而通過圖1和圖2的比較,可以看出CR的變化對各估計的漸近正態(tài)性的影響要比MR的變化所受的影響要大,這個結論和表1的結論一致.

圖1 在例4.1中,在CR=10%和MR=40%,n=500下,從左到右分別表示CA,IW以及IPW估計的正態(tài)Q-Q圖

例4.2考慮一個協(xié)變量為二維的模型

這里的X1i～U(?1,1),X2i～U(?1,1),和εi～N(0,0.52)是相互獨立的.此處的估計目標為條件分位數.選擇二元核函數

圖2 在例4.1中,在CR=40%和MR=10%,n=500下,從左到右分別表示CA,IW以及IPW估計的正態(tài)Q-Q圖

圖3 在例4.1中,在CR=40%和MR=40%,n=500下,從左到右分別表示CA,IW以及IPW估計的正態(tài)Q-Q圖

仍然沿用例4.1的步驟,可以計算Q0.5(x)的各估計.表3選擇CR=10%和40%,MR=10%和40%,n=200和500,基于M=200重復得到關于Q0.5(x)的CA,IW以及IPW估計的GMSE.窗寬hn的取值范圍仍然從0.01到1.0,增量為0.02,分別選擇使GMSE達到最小的窗寬.另外,在x0=(0.4,0.4)上,基于M=200次重復計算各估計下的偏差和標準差.

根據表3中的數據發(fā)現:(1)在其他條件相同下,隨著CA和MR增大,每個估計在GMSE,Bias以及Sd方面的模擬效果會越差;而當樣本容量增大時,各估計的模擬效果越好;(2)在GMSE和Sd方面,CA估計的模擬效果最好;在Bias方面,IW估計模擬的效果比CA和IPW估計的效果要好;(3)CR的變化對各估計模擬效果的影響比MR的變化所帶來的影響要大.這些結論和表1和表2的結論是一致的.

例4.3為了能解決實際問題,選取了肺癌患者臨床數據.在這組數據中,一共有228個肺癌患者,其中刪失的數據為63個.由于這組數據刪失指標是完整的,為了得到本文的觀察數據,對這組數據做一些處理，即對刪失指標進行部分隨機缺失.這樣228個患者中有151個患者的存活時間未被刪失(δ=1,ξ=1),有56個患者的存活時間被刪失了(δ=0,ξ=1),最后還有21個患者的刪失指標被隨機缺失了(ξ=0).另外,作者關心的2項指標是存活時間(time)和患者年齡(age).設T=log(time),X=age/100,X與T的散點圖見圖4中的左圖.

用CC,CA和IW方法對實際數據中的Y=log(time)關于X=age/100作條件中位數曲線(τ=0.5)的模擬(見圖4中的右圖).從圖4種的右圖發(fā)現,用CC和IW兩種方法模擬的效果差不多,在40-50歲之間,患者年齡越大,存活時間反而越長,50-70歲之間,患者存活時間和年齡之間變化不大,到了70歲以后,患者年齡越大,存活時間越短.對于CA方法模擬的曲線似乎更合理,尤其是在40-50歲之間,患者存活時間和年齡之間的變化沒有CC和IW兩種方法那么明顯.

表3 在例4.2中,CC,CA,IW以及IPW估計下的GMSE,Bias以及Sd的模擬結果

圖 4 左圖為例3中X=age/100與T=log(time)之間散點圖;右圖為用CC,CA和IW方法估計下條件中位數曲線的模擬圖

§5 定理的證明

在證明定理前先給出一個引理.