高 偉,馮海林
(西安電子科技大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,西安 710126)
涉及醫(yī)學(xué)、工程等領(lǐng)域的可靠性研究中,剩余壽命分位數(shù)是做出相關(guān)決策的關(guān)鍵指標(biāo)[1],它對(duì)壽命分布的尾部特征性依賴性較小,且在偏態(tài)分布下具有魯棒性的特點(diǎn),能較為全面地刻畫分布特征。右刪失數(shù)據(jù)是一種常見的數(shù)據(jù)類型[2],具有非對(duì)稱和偏態(tài)分布性,而剩余壽命分位數(shù)對(duì)此類數(shù)據(jù)有較好的表現(xiàn)[3]。競(jìng)爭(zhēng)風(fēng)險(xiǎn)數(shù)據(jù)在經(jīng)典生存分析中極為重要,它由個(gè)體面臨失效時(shí)間跨度、導(dǎo)致失效的終止事件等多種潛在結(jié)局產(chǎn)生,其中潛在終止事件被彼此之間稱為“競(jìng)爭(zhēng)風(fēng)險(xiǎn)”事件。
在剩余壽命分位數(shù)的預(yù)測(cè)中,由右刪失數(shù)據(jù)的分位數(shù)回歸方法得到的結(jié)果較準(zhǔn)確。分位數(shù)回歸方法被Koenker[4]首次引入,被Ying等[5]、Gelfand等[6]、Wang等[7]進(jìn)一步推廣;在競(jìng)爭(zhēng)風(fēng)險(xiǎn)數(shù)據(jù)類型下,Jeong和Fine[8]進(jìn)一步提出了對(duì)剩余壽命分位數(shù)推斷的兩樣本檢驗(yàn)統(tǒng)計(jì)量。但這些方法都忽略了協(xié)變量對(duì)預(yù)測(cè)的影響。盡管Jung等[9]在右刪失數(shù)據(jù)下構(gòu)建了協(xié)變量對(duì)剩余壽命分位數(shù)影響的半?yún)?shù)估計(jì)模型,但其中的協(xié)變量被考慮為固定的。實(shí)際上,動(dòng)態(tài)協(xié)變量對(duì)剩余壽命的預(yù)測(cè)也有顯著影響[10]。目前,尚未有動(dòng)態(tài)協(xié)變量對(duì)剩余壽命預(yù)測(cè)影響的研究。
本文基于競(jìng)爭(zhēng)風(fēng)險(xiǎn)下右刪失數(shù)據(jù),建立了剩余壽命分位數(shù)回歸模型,得到了剩余壽命分位數(shù)和動(dòng)態(tài)協(xié)變量之間的函數(shù)關(guān)系,證明了提出的估計(jì)量的漸近性和一致性。數(shù)據(jù)仿真和實(shí)例分析驗(yàn)證了所提出方法的準(zhǔn)確性和有效性。
用T表示個(gè)體的失效時(shí)間,C表示刪失時(shí)間,即從觀測(cè)開始到個(gè)體刪失的時(shí)間,令Yi=min(Ti,Ci),定義示性函數(shù) δi=I(Ti≤Ci),用 ε∈(1,2,…,K)表示失效原因,令ηi=δiε=I(Ti≤Ci)ε。通常把影響個(gè)體失效時(shí)間的重要因素稱之為協(xié)變量,記W∈R(p+1)為固定協(xié)變量向量,即不隨時(shí)間變化,Z∈R(p+1)為與時(shí)間相關(guān)的動(dòng)態(tài)協(xié)變量向量,令U=(W,Z)。假定失效原因?yàn)閗(k=1,…,K)的累積發(fā)病函數(shù)為Fk(t|U)=Pr(T≤t,ε=k|U),它表示給定協(xié)變量U 時(shí),觀察到失效原因?yàn)閗的概率。記在t0時(shí)失效原因?yàn)閗的剩余壽命的子分布函數(shù)為 Fk,t0(t|U)=Pr(T-t0≤t,ε=k|U,T>t0)。假定給定協(xié)變量U時(shí),刪失變量C和(T,ε)是條件獨(dú)立的。設(shè)G0(t|U)=Pr(C>t,εi=k|U)是刪失分布的條件生存函數(shù)。另外,記觀測(cè)到的數(shù)據(jù)為{(Yi,δiεi,Ui),i=1,…,n}。
基于真實(shí)數(shù)據(jù)的分析表明,剩余壽命分位數(shù)與協(xié)變量之間的關(guān)系一般為對(duì)數(shù)線性關(guān)系[11]。設(shè)失效原因?yàn)棣?k的剩余壽命τ-分位數(shù)和協(xié)變量之間的關(guān)系是對(duì)數(shù)線性的,即:
其中 Fk(·)是失效原因?yàn)?ε=k 的CIF ,S(·)是所有原因的生存函數(shù)。
由式(3)得,給定協(xié)變量信息下,失效原因?yàn)棣?k的τ-分位數(shù)函數(shù)表示為:
由于在時(shí)間t0時(shí)τ-分位數(shù)函數(shù)是CIF在Fk(t0|W,時(shí)的逆,因此在t0時(shí),F(xiàn)k(·)和 S(·)都會(huì)影響τ-分位數(shù)。令,其中是給定協(xié)變量信息時(shí)失效原因?yàn)棣?k的累積風(fēng)險(xiǎn)函數(shù),ζ即觀測(cè)時(shí)間的最大值。于是,τ-分位數(shù)存在的必要條件是,即ζ,因此t0被約束。
由式(1)和式(3),知:
假設(shè) (Ti,εi)和Ci是條件獨(dú)立的,則式(4)的左邊第一項(xiàng)等價(jià)于:
而式(4)的左邊的第二項(xiàng)等價(jià)于:
因此,式(4)即為:
其中G(t|W,Z(t0))是給定協(xié)變信息時(shí)刪失分布的條件生存函數(shù)的Kaplan-Meier估計(jì)值[9]。于是:
假設(shè)競(jìng)爭(zhēng)失效原因ε∈{1,2},本文僅對(duì)類型1事件感興趣,即ε=1,通過數(shù)值仿真來檢驗(yàn)所提出方法在有限樣本下的表現(xiàn)??紤]一個(gè)簡(jiǎn)單的情形:
t0=0時(shí)假設(shè) ρ2=0.4,κ1=κ2=1.5,。 用產(chǎn)生類型 1 事件的事件時(shí)間真值產(chǎn)生類型 2事件的事件時(shí)間真值,其中,刪失時(shí)間C~UNIF(0,c),其中c為控制刪失比例的常數(shù),觀察到的生存時(shí)間Yi=min(Ti,Ci) ,其 中假設(shè) w~Bernοulli(0.5),是不隨時(shí)間變化的協(xié)變量,zt0是隨著時(shí)間變化的協(xié)變量,這里采用文獻(xiàn)[10]中的協(xié)變量表達(dá)式zt0=
本文估計(jì)了在刪失比率分別為0、10%、20%、30%的情況下,τ-分位數(shù)分別為0.3、0.5且時(shí)間t0分別為0、0.5、1、1.5時(shí)回歸參數(shù)估計(jì)的經(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í),真值回歸參數(shù)1.40,1.29和的經(jīng)驗(yàn)估計(jì)
表1 t0=0,0.5,1,1.5時(shí),真值回歸參數(shù)1.40,1.29和的經(jīng)驗(yàn)估計(jì)
注: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時(shí)的Bias、SD、SE、CP
表3 τ=0.5時(shí)的Bias、SD、SE、CP
本文利用美國(guó)Channing House數(shù)據(jù)對(duì)所提出的方法進(jìn)行評(píng)估,Channing House是位于美國(guó)加利福尼亞州帕洛阿爾托市的一個(gè)退休中心,Channing House數(shù)據(jù)收集并記錄了從1964年至1975年7月1日之間成員的有關(guān)數(shù)據(jù)。在這期間,總共有97名男性和365名女性在該中心生活。此外,所有成員進(jìn)入和離開退休中心時(shí)的年齡也被記錄。根據(jù)記錄,發(fā)現(xiàn)該數(shù)據(jù)集屬于右刪失數(shù)據(jù)類型,因?yàn)榻Y(jié)束記錄時(shí)還有許多成員依舊存活。僅有46名男性和130名女性在研究期間在Channing House退休中心死亡,由此可得刪失率大約為61.9%。
本文感興趣的是成員的性別差異以及進(jìn)入Channing House退休中心時(shí)的年齡和居住持續(xù)的時(shí)間對(duì)生存時(shí)間的影響。用w=1表示男性個(gè)體,w=0表示女性個(gè)體表示人員進(jìn)入Channing House退休中心時(shí)的年齡和居住持續(xù)的時(shí)間t0。考慮人員的死亡原因?yàn)樗信d趣的類型1事件,考慮回歸模型對(duì)于式(8),在 τ=0.5時(shí)利用網(wǎng)格搜索法求解。
表4總結(jié)了在t0=0、10、20時(shí)的估計(jì)值;表5給出參數(shù)估計(jì)的經(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ù)中可以看出,估計(jì)值有較好的表現(xiàn);表6表明了性別差異對(duì)剩余壽命分位數(shù)的影響,在t0取固定值時(shí),女性比男性生存時(shí)間更長(zhǎng);表7列出了考慮動(dòng)態(tài)協(xié)變量和未考慮動(dòng)態(tài)協(xié)變量條件下的剩余壽命分位數(shù)的估計(jì)值,可以看出成員進(jìn)入退休中心的年齡和居住持續(xù)的時(shí)間對(duì)剩余壽命有一定的影響,這與先前的預(yù)期是一致的。
表4 在t0=0,10,20時(shí)的估計(jì)值
表4 在t0=0,10,20時(shí)的估計(jì)值
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時(shí)的Bias、SD、SE、CP
表6 τ=0.5時(shí),性別差異對(duì)剩余壽命分位數(shù)的影響
表7 τ=0.5時(shí),兩種情形下剩余壽命分位數(shù)的估計(jì)值
由于在生物醫(yī)學(xué)研究領(lǐng)域,生物數(shù)據(jù)經(jīng)常是偏態(tài)分布的,分位數(shù)回歸模型越來越受到研究者的廣泛關(guān)注。本文提出在競(jìng)爭(zhēng)風(fēng)險(xiǎn)下右刪失數(shù)據(jù)的剩余壽命分位數(shù)回歸模型,其主要特點(diǎn)是體現(xiàn)動(dòng)態(tài)協(xié)變量與剩余壽命分位數(shù)的關(guān)系。本文從估計(jì)方程中得到估計(jì)值,并對(duì)其漸進(jìn)性和一致性進(jìn)行推導(dǎo)。對(duì)模型進(jìn)行數(shù)值仿真,證明提出的估計(jì)方法有很好的有限樣本性質(zhì)。最后,將實(shí)際數(shù)據(jù)集應(yīng)用于模型,充分體現(xiàn)了動(dòng)態(tài)協(xié)變量對(duì)于剩余壽命分位數(shù)的影響。下一步研究,可以嘗試研究數(shù)據(jù)類型更為復(fù)雜的情況下,動(dòng)態(tài)協(xié)變量對(duì)于剩余壽命分位數(shù)的影響。