王小剛,陳姜猛

(北方民族大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院,寧夏 銀川 750021)

0 引言

相較于傳統(tǒng)的回歸模型，分位數(shù)回歸模型能夠根據(jù)不同分位數(shù)得到不同的估計結(jié)果，既能衡量在平均意義下的估計結(jié)果，又能反映在極端情況下的估計效果，具有靈活性和穩(wěn)健性以及單調(diào)同變性等優(yōu)點，從而它受到越來越多學(xué)者的關(guān)注.在管理、金融和醫(yī)學(xué)等領(lǐng)域的研究中，由于研究提前結(jié)束、研究對象無法繼續(xù)參與實驗(如研究對象工作變遷、意外身亡等事件不能繼續(xù)研究)會導(dǎo)致無法獲取個體準(zhǔn)確的生存時間，即數(shù)據(jù)存在刪失，所以直接使用分位數(shù)回歸會損失有用的信息，使得模型估計效果產(chǎn)生偏差.文獻(xiàn)[1-2]提出了基于刪失數(shù)據(jù)的分位數(shù)回歸模型估計方法，得到有效的估計;Peng Liming等[3]使用與刪失數(shù)據(jù)相關(guān)聯(lián)的鞅結(jié)構(gòu),估計了刪失分位數(shù)回歸模型的參數(shù);Wang Huixia Judy等[4]在假設(shè)生存時間和刪失變量條件獨立的條件下，提出用局部加權(quán)分位數(shù)回歸方法估計刪失分位數(shù)回歸模型參數(shù);Tang Yanlin等[5]為了自動地估計和檢測分位數(shù)之間的共性,提出了一種新的帶有2種分位數(shù)懲罰變量的刪失分位數(shù)回歸算法;張倩倩等[6]對存在刪失的醫(yī)療費用建立了分位數(shù)回歸模型,提出了簡單加權(quán)估計和模擬退火算法,并驗證了模型具有穩(wěn)健性;李忠桂等[7]基于EL法和SEL法對右刪失數(shù)據(jù)構(gòu)造了更高精度的檢驗統(tǒng)計量;孫桂萍等[8]對長度偏差右刪失數(shù)據(jù)剩余壽命分位數(shù)模型提出了參數(shù)估計,推導(dǎo)了估計的相合性和漸近正態(tài)性;張立文等[9]研究了刪失分位數(shù)回歸模型的變點檢驗問題;王江峰等[10]研究了在刪失指標(biāo)隨機缺失下的非參數(shù)刪失回歸模型,構(gòu)造了回歸函數(shù)的復(fù)合分位數(shù)回歸估計,得到估計的大樣本性質(zhì).

在實際應(yīng)用中,刪失數(shù)據(jù)會受到諸如新政策實施、新技術(shù)變革、經(jīng)濟金融危機、全球新冠肺炎疫情等重大事件沖擊,可能會呈現(xiàn)出非線性的結(jié)構(gòu)特征和趨勢,若不考慮變點的存在性則會產(chǎn)生有偏的估計.在此種情況下,由于分段線性刪失分位數(shù)回歸模型既能解釋模型的非線性結(jié)構(gòu)特征,又能保持刪失分位數(shù)回歸模型原有優(yōu)點,所以該模型就成為研究的熱點,也是使用較為廣泛的模型.

在分段線性刪失分位數(shù)回歸模型中,困難之處在于需要對未知的變點位置進(jìn)行估計.常用的變點位置估計方法可以分為2類.一類方法是2步法，即事先給定某個變點初值，將模型按照變點的位置分成前后2個部分，利用已知方法分別估計參數(shù)，然后通過最優(yōu)化目標(biāo)函數(shù)得到變點估計.如常用的格點搜索法[11-13],該方法的優(yōu)點在于不需要對模型分布做假設(shè)，具有一定的靈活性，然而其缺點在于變點估計的精度與運算效率成反向變動關(guān)系，且變點估計較難給出符合實際意義的解釋，因此常用來做變點位置的初步估計.另一類方法是同時估計方法，即通過某種方法將變點參數(shù)轉(zhuǎn)化為模型參數(shù)，再利用已有方法得到模型參數(shù)及變點位置的估計.如常用的線性化技術(shù)[14-15],線性化技術(shù)簡便易行,但難以推導(dǎo)大樣本性質(zhì),且變點估計效果依賴于迭代快慢與質(zhì)量,常會遇到變點收斂速率過慢狀況.線性化方法的缺陷可以通過光滑化的核函數(shù)方法加以彌補[16-18],光滑化的核函數(shù)方法的思想是將目標(biāo)函數(shù)在變點處展開，利用光滑的核函數(shù)替代在目標(biāo)函數(shù)中的不可導(dǎo)項，從而得到變點位置和模型系數(shù)的估計.光滑化的核函數(shù)方法估計更加有效，收斂速度更快，且能夠得到估計的大樣本性質(zhì)，因此光滑化的核函數(shù)方法被廣泛應(yīng)用于變點的位置估計研究.

本文首先給出了分段線性刪失分位數(shù)回歸模型及估計方法，推導(dǎo)估計量的大樣本性質(zhì)，然后通過數(shù)值模擬驗證了在分段線性刪失分位數(shù)回歸模型中的變點位置及模型系數(shù)估計效果，并將該方法對線性化技術(shù)進(jìn)行了比較.利用藥物濫用數(shù)據(jù)進(jìn)行的實證分析發(fā)現(xiàn)了存在的變點，并給出了合理解釋.

1 分段線性分位數(shù)模型及其估計

1.1 模型定義與估計

給定觀察值(xi,yi,δi),考慮刪失分位數(shù)回歸模型如下:

假設(shè)未知的變點位置為ψ,分段線性刪失分位數(shù)回歸模型定義如下:

(1)

(2)

其中ρτ(k)=k(τ-I(k<0))是檢驗函數(shù).由于變點ψ的存在導(dǎo)致目標(biāo)函數(shù)(2)不可導(dǎo),因此,借鑒文獻(xiàn)[17-18]的方法采用光滑的核函數(shù)K((X-ψ)/h)代替示性函數(shù)I(X>ψ),h(>0)表示窗寬,將核函數(shù)K((X-ψ)/h)代入目標(biāo)函數(shù)(2),得到新的目標(biāo)函數(shù):

(3)

最小化目標(biāo)函數(shù)(3)可得到模型參數(shù)ξ=(β0,β1,β2,β3,γT)T和變點ψ的估計.為了尋求目標(biāo)函數(shù)(3)的最優(yōu)化解,將非線性項β2(Xi-ψ)K((Xi-ψ)/h)在給定初值ψ*處進(jìn)行泰勒展開,即

β2(Xi-ψ)K((Xi-ψ)/h)≈β2Ui+β3Vi,

其中Ui=(Xi-ψ*)K((Xi-ψ*)/h),β3=β2(ψ-ψ*),Vi=-{K((Xi-ψ*)/h)+(Xi-ψ*)K′((Xi-ψ*)/h)/h}.

將Ui、Vi代入目標(biāo)函數(shù)(3),有

(4)

通過最小化目標(biāo)函數(shù)(4)即可更新參數(shù)ξ=(β0,β1,β2,β3,γ)T,估計算法如下.

算法1基于光滑化方法的分段線性刪失分位數(shù)模型估計.

h)/h}.

Step4重復(fù)步驟2和步驟3直到收斂,收斂條件為

參數(shù)估計依賴窗寬的選擇,本文通過最小化交叉驗證函數(shù)選擇最優(yōu)窗寬,即

1.2 漸近性質(zhì)

(5)

定義

q(ωi,θ)=?max{Ci,g(ωi,θ)}/?θ=diag(I(g(ωi,

h)/h)T,

其中θ0為真值,ρ′τ(k)=τ-I(k<0)表示ρτ的1階導(dǎo)數(shù).極小化目標(biāo)函數(shù)(5)等價于求解方程

為了證明其漸近性，給出如下正則性條件.

(A3)Xi的密度函數(shù)Pi(x)在有界緊集合[M1,M2]上連續(xù);

(A4)E(‖Zi‖3)<∞,其中‖·‖表示Euclidean

范數(shù);

(A5)存在2個正定矩陣C0和D0,分別滿足

容易得到?δ>0,?N,當(dāng)n>N時,有

P(Yi-max{Ci,g(ωi,θ)}<ε)<δ,

用核函數(shù)K((X-ψ)/h)代替示性函數(shù)I(X>ψ),有I(Xi>ψ)=K((Xi-ψ)/h)+o(h),h→0.

可得

由正則性假設(shè)(A1)知有

對于有限的N和任意小的δ,當(dāng)n→∞時,上述不等式右邊4項都收斂于0,即

結(jié)合Li Chenxi等[12]證明的一致性，得到ln,τ(θ)的連續(xù)性.引理1得證.

其中

vi(θ,θ0)=ρ′τ(Yi-max{Ci,g(ωi,θ)})q(ωi,θ)-ρ′τ(Yi-max{Ci,g(ωi,θ0)})q(ωi,θ0).

證引理2的證明與文獻(xiàn)[19]中的引理4.6的證明相似.為證明引理2，需驗證文獻(xiàn)[19]中引理4.6的條件(B1)、(B3)和(B5)滿足即可.條件(B1)顯然成立.

對于條件(B3),有

‖vi(θ,θ0(τ))‖≤‖ρ′τ(Yi-max{Ci,g(ωi,θ)})(q(ωi,θ)-q(ωi,θ0))‖+‖(ρ′τ(Yi-max{Ci,g(ωi,θ)})-ρ′τ(Yi-max{Ci,g(ωi,θ0)}))q(ωi,θ0)‖=‖I1i‖+‖I2i‖,

其中

‖I1i‖=‖ρ′τ(Yi-max{Ci,g(ωi,θ)})(q(ωi,θ)-q(ωi,θ0))‖=‖ρ′τ(Yi-max{Ci,g(ωi,θ)})q′(ωi,θ0)(θ-θ0)‖=op(1),

所以有E(‖I1i‖2|ωi)=op(1).同時,

‖I2i‖=‖(ρ′τ(Y-max{Ci,g(ωi,θ)})-ρ′τ(Yi-max{Ci,g(ωi,θ0)}))q(ωi,θ0)‖=‖(I(Yi≤max{Ci,

g(ωi,θ)})-I(Yi≤max{Ci,g(ωi,θ0)}))q(ωi,θ0)‖≤‖q(ωi,θ0)‖I(max{Ci,g1(ωi,θ,θ0)}≤Yi≤max{Ci,g2(ωi,θ,θ0)}),

這里的g1(ωi,θ,θ0)和g2(ωi,θ,θ0)分別為

g1(ωi,θ,θ0)=min{max{Ci,g(ωi,θ)},max{Ci,g(ωi,θ0)}},

g2(ωi,θ,θ0)=max{max{Ci,g(ωi,θ)},max{Ci,g(ωi,θ0)}}.

又由于g1(ωi,θ,θ0)≤g2(ωi,θ,θ0),有

‖max{Ci,g(ωi,θ)}-max{Ci,g(ωi,θ0)}‖=‖qT(ωi,θ0)(θ-θ0)+Rn‖≤‖q(ωi,θ0)(θ-θ0)‖≤an‖q(ωi,θ0)‖,

即

E(‖I2i‖2|ωi)≤‖q(ωi,θ0)‖2E(I(g1(ωi,θ,θ0)≤Yi≤g2(ωi,θ,θ0)))≤‖q(ωi,θ0)‖2·fτ,ωi(ξi)‖max{Ci,g(ωi,θ)}-max{Ci,g(ωi,θ0)}‖≤anfτ,ωi(ξi)‖q(ωi,θ0)‖3,

即條件(B3)滿足.

下面驗證條件(B5).對于任意的常數(shù)C,令

有

定理1在正則條件(A1)～(A5)下,有

其中協(xié)方差陣中的C0和D0的估計分別為

Δn是平滑參數(shù),借鑒文獻(xiàn)[21]的方法選擇如下:

其中φ(·)和Φ(·)分別表示標(biāo)準(zhǔn)正態(tài)分布的密度函數(shù)和分布函數(shù).

證在引理1和引理2的條件下,可得

因此有

2 蒙特卡羅模擬

Y=max{C,T},δ=I(T≥C),QT(τ,θ|X,Z)=β0+β1X+β2(X-ψ)++ZTγ.

在固定刪失和隨機刪失下、在不同分位數(shù)下的同方差和異方差數(shù)值模擬結(jié)果分別如表1和表2所示.表1和表2的上部分和下部分分別為在樣本量n=200和n=500下估計量的評價指標(biāo),即偏差(Bs)、標(biāo)準(zhǔn)差(SD) 、標(biāo)準(zhǔn)誤(ES) 、95%的覆蓋率(CP) 和置信區(qū)間長度(AW).

從表1可以看出:在固定刪失情況下，當(dāng)分位數(shù)為0.5且同方差時，模型估計的Bs和SD都很小，且ES和SD非常接近，這說明估計效果是有效的，AW較小，且CP都接近于95%.從橫向來看，在0.5分位數(shù)下，異方差的參數(shù)估計Bs和SD均較小，ES和SD很接近，這說明即使在異方差情況下模型估計效果也較好.另外,異方差的置信區(qū)間長度AW和覆蓋率CP的估計結(jié)果要比在同方差情況下的估計結(jié)果更差.當(dāng)分位數(shù)為0.3和0.7時，模型估計的Bs、SD要比在分位數(shù)為0.5時的估計結(jié)果更差，ES和SD的接近程度較好，AW與CP的效果與在分位數(shù)為0.5時的效果差異不大，這說明估計結(jié)果具有穩(wěn)健性.

當(dāng)樣本量從200增長到500時，可以發(fā)現(xiàn)無論是估計的Bs、SD、ES還是AW都呈現(xiàn)減小趨勢，但CP呈現(xiàn)增大趨勢，這也驗證了定理1的結(jié)果.同樣，ES與SD非常接近，CP更接近于95%，這也驗證了估計具有可靠性.

表1 在固定刪失下的模擬結(jié)果

在隨機刪失下的模擬結(jié)果如表2所示.從表2可以看出:在同方差、相同樣本容量、相同分位數(shù)情況下，固定刪失的估計效果要比隨機刪失的估計效果更好.Bs、SD、ES之間的關(guān)系以及AW和CP之間的特征與表1中所得的結(jié)論相同.

表2 在隨機刪失下的模擬結(jié)果

表2(續(xù))

總體上，在其他情況不變的條件下，同方差下的估計效果要優(yōu)于異方差下的估計效果，更大的樣本容量的估計效果更好,這驗證了參數(shù)估計滿足相合性，固定刪失的估計效果優(yōu)于隨機刪失的估計效果，0.5的分位數(shù)的估計效果優(yōu)于其他分位數(shù)的估計效果.蒙特卡羅模擬結(jié)果表明本文的光滑化方法具有有效性和穩(wěn)健性.

為對比光滑化方法的可行性，對本文方法和線性化技術(shù)得到的模擬結(jié)果進(jìn)行比較.當(dāng)樣本容量為1 000且分位數(shù)為0.5時的模擬結(jié)果如表3和表4所示.

表3 在固定刪失下線性化技術(shù)與本文方法的模擬結(jié)果對比

表4 在隨機刪失下線性化技術(shù)與本文方法的模擬結(jié)果對比

從表3和表4中可以看出所有的Bs都非常小，因此2種方法的估計是漸近一致的.2種方法的ES與SD非常接近，這表明估計具有相合性.此外，2種方法的CP均接近95%.但在相同顯著性水平下，本文方法在同方差情況下略好于線性化技術(shù)，在異方差情況下優(yōu)于線性化技術(shù)，本文方法估計的AW比線性化技術(shù)估計的AW更小，因此，本文方法在估計變點位置及模型參數(shù)的有效性上要優(yōu)于線性化技術(shù)，且本文方法的收斂速率較高,而線性化技術(shù)在模擬時可能會存在不收斂情形.

3 實證分析

在醫(yī)學(xué)領(lǐng)域中，藥物的濫用會給個人、家庭及社會帶來巨大的危害性.若不采取有效的預(yù)防措施和控制，則與之有關(guān)的疾病將會很快在全球泛濫成災(zāi)，所有國家都將處于這種危險之中，因此研究藥物濫用復(fù)發(fā)時間的影響因素和變化趨勢有著重要的實際意義，也有利于對藥物濫用復(fù)發(fā)時間的控制和采取有效的預(yù)防措施.但是在對藥物濫用復(fù)發(fā)治療時，由于病人的藥物沒有復(fù)發(fā)而研究結(jié)束了或在研究結(jié)束前有病人提前結(jié)束治療，所以導(dǎo)致數(shù)據(jù)存在刪失.另外,考慮到模型會受到藥物復(fù)發(fā)治療的影響因素以及藥物復(fù)發(fā)的反應(yīng)快慢的沖擊而產(chǎn)生變點，因此可通過用含變點的分段線性刪失分位數(shù)回歸模型研究藥物濫用復(fù)發(fā)時間的影響因素，并利用本文方法進(jìn)行估計.

文獻(xiàn)[22]給出了關(guān)于UIS藥物治療的研究數(shù)據(jù)，原始數(shù)據(jù)的樣本容量為628，包含575個無缺失的數(shù)據(jù).響應(yīng)變量為復(fù)發(fā)時間間隔，協(xié)變量為隨機化治療ZT,ZT=1表示接受了180 d的治療,ZT=0表示接受了90 d的治療，ZA為登記年齡，ZB為抑郁評分(取值為0～54)，ZV為最近靜脈用藥史，ZV=1表示用藥，ZV=0表示沒有用藥，ZN表示既往藥物治療次數(shù)(取值0～40),由于ZN在一定情況下會受到治療效應(yīng)解釋的影響，所以將ZN分成ZN1和ZN22個指標(biāo)，ZR表示種族，ZR=0為白人，ZR=1表示非白人，ZS表示治療地點，ZS=0表示治療點A，ZS=1表示治療點B，XF表示治療時間，用治療天數(shù)TL定義，XF=TL/90為短治療，XF=TL/180為長治療.在研究結(jié)束前有病人提前結(jié)束治療，因此數(shù)據(jù)存在刪失.選取協(xié)變量ZN1、ZN2、ZV、ZT、XF[23]構(gòu)建分段線性刪失分位數(shù)回歸模型.

從藥物復(fù)發(fā)時間間隔(取對數(shù))和XF的散點圖(見圖1)可以看出:XF值在0.5前后復(fù)發(fā)時間間隔呈現(xiàn)不同特點.在0.5之前復(fù)發(fā)時間間隔較大，但在0.5之后逐漸變小，呈現(xiàn)出非線性特征,這表明模型可能存在變點.因此建立分段線性刪失分位數(shù)回歸來分析藥物濫用復(fù)發(fā)時間的治療數(shù)據(jù)，模型如下：

Yi=max(Ci,Ti),δi=I(Ti≥Ci),i=1,2,…,575,

QTi(τ|Xi,Zi)=β0+β1XFi+β2(XFi-ψ)++γ1ZN1i+γ2ZN2i+γ3ZVi+γ4ZTi,

其中Yi為藥物復(fù)發(fā)時間間隔，δi為刪失指標(biāo)，θ=(β0,β1,β2,γ1,γ2,γ3,γ4,ψ)T為未知參數(shù).在分位數(shù)τ=0.3、0.5、0.7時利用本文方法估計了變點位置及模型參數(shù)，模型估計結(jié)果如表5所示.

由表5可知：在0.5的分位數(shù)下，所有參數(shù)估計的P值都小于0.05，通過了顯著性檢驗.變點系數(shù)β不為0表示存在變點，變點位置為0.498.其他協(xié)變量的系數(shù)都為正，這表示其他協(xié)變量與復(fù)發(fā)時間間隔正相關(guān).在0.3分位數(shù)下的估計與在0.5分位數(shù)下的估計相同，變點位置為0.485，與在0.5分位數(shù)下的情況相似.在0.7的分位數(shù)下，隨機化治療ZT沒通過顯著性檢驗，其他協(xié)變量都通過了檢驗，且也存在變點，變點位置為0.814.在不同分位數(shù)下的分段線性刪失回歸模型的估計都表明:復(fù)發(fā)時間間隔與治療時間之間存在變點，前一半治療時間(在0.5的分位數(shù)時為0.498)與后一半治療時間的治療效果有顯著差異，即前一半治療時間的復(fù)發(fā)時間間隔更長.

表5 估計結(jié)果

為更清晰地得到復(fù)發(fā)時間間隔與治療時間之間的變化特征，在復(fù)發(fā)時間和XF散點圖上畫出不同分位數(shù)下的折線圖(見圖1).

圖1 藥物復(fù)發(fā)時間間隔與治療時間散點圖

從圖1可以看出：當(dāng)分位數(shù)為0.3和0.5時，治療時間在0.5附近(估計值分別為0.485和0.498)有顯著差異.當(dāng)分位數(shù)為0.7時，治療時間在0.814處存在變點.變點前的復(fù)發(fā)時間間隔隨著治療時間的增加而增加，變點后的復(fù)發(fā)時間間隔明顯比變點前的復(fù)發(fā)時間更小，也就是說，在前半段治療時間里復(fù)發(fā)時間間隔更長，即變點之前的治療比變點之后的治療更有效.

4 結(jié)論

本文基于光滑化的核函數(shù)方法研究了分段線性刪失分位數(shù)回歸模型變點位置及模型系數(shù)的估計問題，利用核函數(shù)替代在目標(biāo)函數(shù)中的非光滑項，通過迭代算法得到變點位置和模型系數(shù)估計，推導(dǎo)了估計量的漸近性質(zhì).蒙特卡羅模擬結(jié)果驗證了估計效果的有效性和穩(wěn)健性,同時利用該方法與線性化技術(shù)進(jìn)行了比較，發(fā)現(xiàn)在置信區(qū)間較小的情況下，該方法的估計效果更好.針對藥物濫用的數(shù)據(jù)實證分析表明：復(fù)發(fā)時間間隔與治療時間存在正向影響，復(fù)發(fā)時間在0.5的分位數(shù)下存在變點，變點位置為0.498，即前一半時間的治療比后一半時間的治療更加有效.在0.3和0.7的分位數(shù)下也發(fā)現(xiàn)變點，且變點前后的治療效果與在0.5分位數(shù)時的治療效果相近.