高 偉，馮海林

（西安電子科技大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，西安 710126）

0 引言

涉及醫(yī)學(xué)、工程等領(lǐng)域的可靠性研究中，剩余壽命分位數(shù)是做出相關(guān)決策的關(guān)鍵指標(biāo)[1]，它對壽命分布的尾部特征性依賴性較小，且在偏態(tài)分布下具有魯棒性的特點(diǎn)，能較為全面地刻畫分布特征。右刪失數(shù)據(jù)是一種常見的數(shù)據(jù)類型[2]，具有非對稱和偏態(tài)分布性，而剩余壽命分位數(shù)對此類數(shù)據(jù)有較好的表現(xiàn)[3]。競爭風(fēng)險數(shù)據(jù)在經(jīng)典生存分析中極為重要，它由個體面臨失效時間跨度、導(dǎo)致失效的終止事件等多種潛在結(jié)局產(chǎn)生，其中潛在終止事件被彼此之間稱為“競爭風(fēng)險”事件。

在剩余壽命分位數(shù)的預(yù)測中，由右刪失數(shù)據(jù)的分位數(shù)回歸方法得到的結(jié)果較準(zhǔn)確。分位數(shù)回歸方法被Koenker[4]首次引入，被Ying等[5]、Gelfand等[6]、Wang等[7]進(jìn)一步推廣；在競爭風(fēng)險數(shù)據(jù)類型下，Jeong和Fine[8]進(jìn)一步提出了對剩余壽命分位數(shù)推斷的兩樣本檢驗(yàn)統(tǒng)計量。但這些方法都忽略了協(xié)變量對預(yù)測的影響。盡管Jung等[9]在右刪失數(shù)據(jù)下構(gòu)建了協(xié)變量對剩余壽命分位數(shù)影響的半?yún)?shù)估計模型，但其中的協(xié)變量被考慮為固定的。實(shí)際上，動態(tài)協(xié)變量對剩余壽命的預(yù)測也有顯著影響[10]。目前，尚未有動態(tài)協(xié)變量對剩余壽命預(yù)測影響的研究。

本文基于競爭風(fēng)險下右刪失數(shù)據(jù)，建立了剩余壽命分位數(shù)回歸模型，得到了剩余壽命分位數(shù)和動態(tài)協(xié)變量之間的函數(shù)關(guān)系，證明了提出的估計量的漸近性和一致性。數(shù)據(jù)仿真和實(shí)例分析驗(yàn)證了所提出方法的準(zhǔn)確性和有效性。

1 符號和定義

用T表示個體的失效時間，C表示刪失時間，即從觀測開始到個體刪失的時間，令Yi=min(Ti，Ci)，定義示性函數(shù) δi=I(Ti≤Ci)，用 ε∈(1，2，…，K)表示失效原因，令ηi=δiε=I(Ti≤Ci)ε。通常把影響個體失效時間的重要因素稱之為協(xié)變量，記W∈R(p+1)為固定協(xié)變量向量，即不隨時間變化，Z∈R(p+1)為與時間相關(guān)的動態(tài)協(xié)變量向量，令U=(W，Z)。假定失效原因?yàn)閗(k=1，…，K)的累積發(fā)病函數(shù)為Fk(t|U)=Pr(T≤t，ε=k|U)，它表示給定協(xié)變量U 時，觀察到失效原因?yàn)閗的概率。記在t0時失效原因?yàn)閗的剩余壽命的子分布函數(shù)為 Fk，t0(t|U)=Pr(T-t0≤t，ε=k|U，T＞t0)。假定給定協(xié)變量U時，刪失變量C和(T，ε)是條件獨(dú)立的。設(shè)G0(t|U)=Pr(C＞t，εi=k|U)是刪失分布的條件生存函數(shù)。另外，記觀測到的數(shù)據(jù)為{(Yi，δiεi，Ui)，i=1，…，n}。

2 基于動態(tài)協(xié)變量的剩余壽命分位數(shù)建模及其參數(shù)估計

基于真實(shí)數(shù)據(jù)的分析表明，剩余壽命分位數(shù)與協(xié)變量之間的關(guān)系一般為對數(shù)線性關(guān)系[11]。設(shè)失效原因?yàn)棣?k的剩余壽命τ-分位數(shù)和協(xié)變量之間的關(guān)系是對數(shù)線性的，即：

其中 Fk(·)是失效原因?yàn)?ε=k 的CIF ，S(·)是所有原因的生存函數(shù)。

由式（3）得，給定協(xié)變量信息下，失效原因?yàn)棣?k的τ-分位數(shù)函數(shù)表示為：

由于在時間t0時τ-分位數(shù)函數(shù)是CIF在Fk(t0|W，時的逆，因此在t0時，F(xiàn)k(·)和 S(·)都會影響τ-分位數(shù)。令，其中是給定協(xié)變量信息時失效原因?yàn)棣?k的累積風(fēng)險函數(shù)，ζ即觀測時間的最大值。于是，τ-分位數(shù)存在的必要條件是，即ζ，因此t0被約束。

由式（1）和式（3），知：

假設(shè) (Ti，εi)和Ci是條件獨(dú)立的，則式（4）的左邊第一項等價于：

而式（4）的左邊的第二項等價于：

因此，式（4）即為：

其中G(t|W，Z(t0))是給定協(xié)變信息時刪失分布的條件生存函數(shù)的Kaplan-Meier估計值[9]。于是：

3 仿真模擬

假設(shè)競爭失效原因ε∈{1，2}，本文僅對類型1事件感興趣，即ε=1，通過數(shù)值仿真來檢驗(yàn)所提出方法在有限樣本下的表現(xiàn)?？紤]一個簡單的情形：

t0=0時假設(shè) ρ2=0.4，κ1=κ2=1.5，。 用產(chǎn)生類型 1 事件的事件時間真值產(chǎn)生類型 2事件的事件時間真值，其中，刪失時間C～UNIF(0，c)，其中c為控制刪失比例的常數(shù)，觀察到的生存時間Yi=min(Ti，Ci) ，其 中假設(shè) w～Bernοulli(0.5)，是不隨時間變化的協(xié)變量，zt0是隨著時間變化的協(xié)變量，這里采用文獻(xiàn)[10]中的協(xié)變量表達(dá)式zt0=

本文估計了在刪失比率分別為0、10%、20%、30%的情況下，τ-分位數(shù)分別為0.3、0.5且時間t0分別為0、0.5、1、1.5時回歸參數(shù)估計的經(jīng)驗(yàn)偏差（Bias）、經(jīng)驗(yàn)標(biāo)準(zhǔn)差（SD）、平均標(biāo)準(zhǔn)差（SE）和95%Wald型置信區(qū)間的經(jīng)驗(yàn)覆蓋率（CP）。

表1 t0=0，0.5，1，1.5時，真值回歸參數(shù)1.40，1.29和的經(jīng)驗(yàn)估計

表1 t0=0，0.5，1，1.5時，真值回歸參數(shù)1.40，1.29和的經(jīng)驗(yàn)估計

注：c%表示刪失比例。

t0 0 c%β1，t0 1，t0 β1，t0 True 1，t0 1.5 1，t0 1，t0 α(1)0.5 0.0937 0.0963 0.0789 0.1175 0.1381 0.1187 0.1665 0.1278 0.0956 0.1479 0.1268 0.1329 0.1768 0.1120 0.1407 0.1586 α(1)10 20 30 0 10 20 30 0 10 20 30 0 10 20 30 τ=0.3 α(0)1.61 1.61 1.61 1.61 1.51 1.51 1.51 1.51 1.40 1.40 1.40 1.40 1.29 1.29 1.29 1.29 τ=0.5 α(0)1 0 1.5813 1.5746 1.5696 1.5499 1.4670 1.4422 1.4106 1.3942 1.3768 1.3756 1.3510 1.3119 1.2170 1.1494 1.1092 0.9422 α(1)1，t0 0.0944 0.0881 0.0839 0.1172 0.1837 0.1334 0.1718 0.1457 0.1236 0.1438 0.1181 0.1723 0.1994 0.1891 0.1952 0.1178 1.6445 1.6471 1.6543 1.6885 1.4919 1.5063 1.5238 1.5425 1.3814 1.4276 1.4453 1.4733 1.2304 1.2816 1.4012 1.4717 0.0615 0.0933 0.0751 0.0751 0.1352 0.1127 0.0694 0.1100 0.1151 0.0718 0.1014 0.1458 0.1073 0.1447 0.1051 0.1274 0.0615 0.0933 0.0751 0.0751 0.1352 0.1127 0.0694 0.1100 0.1151 0.0718 0.1014 0.1458 0.1073 0.1447 0.1051 0.1274 0.0387 0.1019 0.1053 0.0627 0.1880 0.1454 0.1135 0.1255 0.1097 0.1491 0.1353 0.1342 0.1902 0.1367 0.1470 0.1737

表2 τ=0.3時的Bias、SD、SE、CP

表3 τ=0.5時的Bias、SD、SE、CP

4 實(shí)例應(yīng)用

本文利用美國Channing House數(shù)據(jù)對所提出的方法進(jìn)行評估，Channing House是位于美國加利福尼亞州帕洛阿爾托市的一個退休中心，Channing House數(shù)據(jù)收集并記錄了從1964年至1975年7月1日之間成員的有關(guān)數(shù)據(jù)。在這期間，總共有97名男性和365名女性在該中心生活。此外，所有成員進(jìn)入和離開退休中心時的年齡也被記錄。根據(jù)記錄，發(fā)現(xiàn)該數(shù)據(jù)集屬于右刪失數(shù)據(jù)類型，因?yàn)榻Y(jié)束記錄時還有許多成員依舊存活。僅有46名男性和130名女性在研究期間在Channing House退休中心死亡，由此可得刪失率大約為61.9%。

本文感興趣的是成員的性別差異以及進(jìn)入Channing House退休中心時的年齡和居住持續(xù)的時間對生存時間的影響。用w=1表示男性個體，w=0表示女性個體表示人員進(jìn)入Channing House退休中心時的年齡和居住持續(xù)的時間t0。考慮人員的死亡原因?yàn)樗信d趣的類型1事件，考慮回歸模型對于式（8），在 τ=0.5時利用網(wǎng)格搜索法求解。

表4總結(jié)了在t0=0、10、20時的估計值；表5給出參數(shù)估計的經(jīng)驗(yàn)偏差（Bias）、經(jīng)驗(yàn)標(biāo)準(zhǔn)差（SD）、平均標(biāo)準(zhǔn)差（SE）和95%Wald型置信區(qū)間的經(jīng)驗(yàn)覆蓋率（CP），從這些數(shù)據(jù)中可以看出，估計值有較好的表現(xiàn)；表6表明了性別差異對剩余壽命分位數(shù)的影響，在t0取固定值時，女性比男性生存時間更長；表7列出了考慮動態(tài)協(xié)變量和未考慮動態(tài)協(xié)變量條件下的剩余壽命分位數(shù)的估計值，可以看出成員進(jìn)入退休中心的年齡和居住持續(xù)的時間對剩余壽命有一定的影響，這與先前的預(yù)期是一致的。

表4 在t0=0，10，20時的估計值

表4 在t0=0，10，20時的估計值

t0 01 0 20 τ=0.5 α(0)1，t0 4.1963 4.0020 3.7430 α(1)1，t0-0.3948 0.0655 0.0920 β1，t0 0.2421 0.0492 0.0439

表5 τ=0.5時的Bias、SD、SE、CP

表6 τ=0.5時，性別差異對剩余壽命分位數(shù)的影響

表7 τ=0.5時，兩種情形下剩余壽命分位數(shù)的估計值

5 結(jié)束語

由于在生物醫(yī)學(xué)研究領(lǐng)域，生物數(shù)據(jù)經(jīng)常是偏態(tài)分布的，分位數(shù)回歸模型越來越受到研究者的廣泛關(guān)注。本文提出在競爭風(fēng)險下右刪失數(shù)據(jù)的剩余壽命分位數(shù)回歸模型，其主要特點(diǎn)是體現(xiàn)動態(tài)協(xié)變量與剩余壽命分位數(shù)的關(guān)系。本文從估計方程中得到估計值，并對其漸進(jìn)性和一致性進(jìn)行推導(dǎo)。對模型進(jìn)行數(shù)值仿真，證明提出的估計方法有很好的有限樣本性質(zhì)。最后，將實(shí)際數(shù)據(jù)集應(yīng)用于模型，充分體現(xiàn)了動態(tài)協(xié)變量對于剩余壽命分位數(shù)的影響。下一步研究，可以嘗試研究數(shù)據(jù)類型更為復(fù)雜的情況下，動態(tài)協(xié)變量對于剩余壽命分位數(shù)的影響。