截尾分位數(shù)回歸及其在生存分析中的應(yīng)用*

王 娟 王 彤

在傳統(tǒng)的生存分析中，多因素分析常采用Cox比例風(fēng)險(xiǎn)回歸模型和加速失效時(shí)間模型，Cox比例風(fēng)險(xiǎn)回歸模型要求資料滿足等比例風(fēng)險(xiǎn)假定和對(duì)數(shù)線性假定，加速失效時(shí)間模型要求資料分布已知。但醫(yī)學(xué)研究中，某些資料不滿足上述要求，特別當(dāng)觀察值中有離群點(diǎn)(outliers)，所擬合的回歸方程因“遷就”這個(gè)(些)離群點(diǎn)而使整體的擬合結(jié)果產(chǎn)生不同程度的偏離，以致影響了穩(wěn)定性;再者，當(dāng)總體存在異方差，即y的條件分布中方差不為常數(shù)時(shí);或者研究興趣在于分布在尾端區(qū)域數(shù)據(jù)的建模，亦不宜用Cox比例風(fēng)險(xiǎn)回歸模型和加速失效時(shí)間模型。本文介紹一種適合于當(dāng)生存資料不滿足上述條件時(shí)的回歸方法——截尾分位數(shù)回歸。

原理與方法

所謂截尾分位數(shù)回歸(censored regression quantiles，CRQ)，就是根據(jù)估計(jì)(生存時(shí)間的對(duì)數(shù)形式)的分位數(shù)，其模型為:

QY|X(τ，x)=xβ(τ)

與加速失效時(shí)間模型不同的是，QY|X(τ，x)表示給定x的條件下，y的第τ分位數(shù)的估計(jì)值。其中τ的取值范圍為0＜τ＜1，當(dāng)τ=0.5時(shí)，截尾分位數(shù)回歸就是中位數(shù)回歸(median regression)。截尾分位數(shù)回歸中，參數(shù)估計(jì)一般用加權(quán)最小一乘(weighted least absolute，WLA)準(zhǔn)則，即使目標(biāo)函數(shù):

模擬分析

1．滿足參數(shù)模型假定的模擬分析

如果生存時(shí)間服從指數(shù)分布，既能建立加速失效時(shí)間模型log(t)=+x+ε，又能建立比例風(fēng)險(xiǎn)回歸模型 h(t，x)=h0(t)exp(x'β)，且兩個(gè)模型中回歸參數(shù)的相互關(guān)系為= － β〔2〕。

模擬方法如下〔3〕:模擬的樣本含量為200，截尾比例為20%、40%、60%。首先產(chǎn)生服從(0，1)均勻分布的隨機(jī)數(shù)S，令生存函數(shù)S(t)=S，同時(shí)產(chǎn)生服從N(0，1)的隨機(jī)數(shù)作為自變量 x。令 β=1，利用 t=產(chǎn)生相應(yīng)的生存時(shí)間t。產(chǎn)生服從B(200，0.2)的隨機(jī)數(shù)作為截尾指示變量，這時(shí)截尾比例為20%。每次得到1000個(gè)模擬樣本，分別采用Cox比例風(fēng)險(xiǎn)回歸模型(Cox)、加速失效時(shí)間模型(AFT)和截尾分位數(shù)回歸模型(CRQ)，最后計(jì)算出這1000個(gè)樣本的回歸系數(shù)的均數(shù)和標(biāo)準(zhǔn)差。模擬結(jié)果見(jiàn)表1。

表格中出現(xiàn)缺項(xiàng)的原因是因?yàn)殡S著截尾比例的不斷增加，越來(lái)越多的截尾生存時(shí)間會(huì)隨機(jī)地出現(xiàn)在生存時(shí)間分布的上游，直到觀察結(jié)束時(shí)，生存率未降到0，即生存分布不完全，無(wú)法估計(jì)高分位數(shù)水平對(duì)應(yīng)的回歸方程。

模擬結(jié)果顯示:對(duì)滿足比例風(fēng)險(xiǎn)假定和對(duì)數(shù)線性假定的數(shù)據(jù)，Cox比例風(fēng)險(xiǎn)回歸模型和截尾分位數(shù)回歸模型擬合的回歸系數(shù)與真值1非常接近，估計(jì)效率也相似。由于資料中不存在異質(zhì)問(wèn)題，所以不同分位數(shù)水平對(duì)應(yīng)的回歸系數(shù)非常接近。隨著截尾比例的不斷增加，回歸系數(shù)可能會(huì)偏離真值1，估計(jì)的效率逐漸降低。該數(shù)據(jù)模擬實(shí)際上符合指數(shù)回歸模型，采用加速失效時(shí)間模型估計(jì)的回歸系數(shù)與真值1最接近，標(biāo)準(zhǔn)差最小，效率最高。

表1 不同截尾比例下回歸系數(shù)的模擬分析

2．參數(shù)條件下存在離群點(diǎn)的模擬分析

模擬方法如下:模擬的樣本含量為200，模擬的截尾比例分別為20%，40%，60%，模擬的離群點(diǎn)比例為5%，10%，15%，20%。首先產(chǎn)生服從(0.5，1.5)均勻分布的隨機(jī)數(shù)作為自變量，再產(chǎn)生服從W(1，1)的威布爾分布隨機(jī)數(shù)作為誤差項(xiàng)，產(chǎn)生服從B(200，0.2)的隨機(jī)數(shù)作為截尾指示變量，這時(shí)截尾比例為20%。如擬定的模型是無(wú)離群點(diǎn)模型，利用log t=1+βx+ε(令β=1)便可以得到相應(yīng)的對(duì)數(shù)生存時(shí)間log t。如擬定的模型是有離群點(diǎn)模型，離群點(diǎn)比例為5%，先從200例中隨機(jī)抽取5%的對(duì)數(shù)生存時(shí)間數(shù)據(jù)，被抽中的對(duì)數(shù)生存時(shí)間數(shù)值在原值的基礎(chǔ)上加5，作為離群點(diǎn)的對(duì)數(shù)生存時(shí)間，而原樣本中其他個(gè)體的對(duì)數(shù)生存時(shí)間保持不變。由此構(gòu)造離群點(diǎn)比例為5%的有離群點(diǎn)樣本。每次得到1000個(gè)模擬樣本，分別采用Cox比例風(fēng)險(xiǎn)回歸模型(Cox)、加速失效時(shí)間模型(AFT)以及截尾分位數(shù)回歸模型(CRQ)，最后計(jì)算出這1000個(gè)樣本的回歸系數(shù)的均數(shù)和標(biāo)準(zhǔn)差。模擬結(jié)果見(jiàn)表2～4。

表2 不同比例離群點(diǎn)下回歸系數(shù)的模擬比較(截尾比例=20%)

表3 不同比例離群點(diǎn)下回歸系數(shù)的模擬比較(截尾比例=40%)

模擬結(jié)果表明，當(dāng)離群點(diǎn)存在時(shí)，Cox比例風(fēng)險(xiǎn)回歸模型擬合的回歸系數(shù)明顯偏離真值1，出現(xiàn)了偏性，并且隨著離群點(diǎn)比例不斷增加，偏性越明顯。當(dāng)截尾比例較低(≤40%)時(shí)，截尾分位數(shù)回歸模型(τ=0.5)擬合的回歸系數(shù)較加速失效時(shí)間模型接近真值1，效率也比較高。但是截尾比例繼續(xù)增加時(shí)，截尾分位數(shù)回歸(τ=0.5)擬合的回歸系數(shù)偏性比加速失效時(shí)間模型大，標(biāo)準(zhǔn)差也偏大，這是因?yàn)殡S著截尾比例的增加，就容易在生存分布的上游出現(xiàn)一定比例的截尾，從而造成了在較高水平的分位數(shù)函數(shù)估計(jì)有偏，效率也降低，但是在低水平的分位數(shù)函數(shù)擬合良好。

表4 不同比例離群點(diǎn)下回歸系數(shù)的模擬比較(截尾比例=60%)

3．違背比例風(fēng)險(xiǎn)假定的模擬分析

有理論表明〔4〕，在參數(shù)加速失效時(shí)間模型家族中，只有威布爾回歸模型(包括指數(shù)分布)既屬于加速失效時(shí)間模型，又屬于比例風(fēng)險(xiǎn)模型。對(duì)數(shù)正態(tài)回歸模型和對(duì)數(shù)logistic回歸模型不滿足比例風(fēng)險(xiǎn)假定。此次模擬采用對(duì)數(shù)正態(tài)回歸模型來(lái)產(chǎn)生數(shù)據(jù)。

模擬方法如下:模擬的樣本含量為200，模擬的截尾比例分別為20%，40%，60%。首先產(chǎn)生服從(0，1)均勻分布的隨機(jī)數(shù)作為自變量x，再產(chǎn)生服從N(0，1)正態(tài)分布的隨機(jī)數(shù)作為誤差項(xiàng)ε，產(chǎn)生服從B(200，0.2)的隨機(jī)數(shù)作為截尾指示變量，這時(shí)截尾比例為20%。利用t=exp(1+βx+ε)(令β=1)便可以得到相應(yīng)的生存時(shí)間t。每次得到1000個(gè)模擬樣本，分別采用Cox比例風(fēng)險(xiǎn)回歸模型(Cox)、截尾分位數(shù)回歸模型(CRQ)，最后計(jì)算出這1000個(gè)樣本的回歸系數(shù)的均數(shù)和標(biāo)準(zhǔn)差。模擬結(jié)果如下(在模擬結(jié)果中，考慮到兩模型中的回歸系數(shù)不具有可比性，便將Cox模型中的回歸系數(shù)按式QCox(τ|x)=H－10(－log(1－τ)exp(－x'iβ))進(jìn)行了轉(zhuǎn)換，表5中列出的就是轉(zhuǎn)換后的結(jié)果)。模擬結(jié)果見(jiàn)表5。

模擬結(jié)果表明，Cox估計(jì)的回歸系數(shù)有偏;CRQ擬合的結(jié)果接近真值1，且效率相比Cox模型偏高。隨著截尾比例的不斷增加，CRQ估計(jì)的效率有下降的趨勢(shì)。

實(shí)例分析

導(dǎo)尿及留置導(dǎo)尿管是臨床上診斷、治療各種危重病人的常用護(hù)理措施之一，但長(zhǎng)期留置導(dǎo)尿管的多數(shù)病人會(huì)不同程度地出現(xiàn)導(dǎo)尿管引流不暢及尿液從導(dǎo)尿管滲漏等問(wèn)題，且并發(fā)癥較多，其中最嚴(yán)重的是尿路感染。因此尋找尿管誘發(fā)尿路感染的影響因素是我們迫切需要解決的問(wèn)題，從而為預(yù)防感染提供一定的科學(xué)依據(jù)。

表5 不同截尾比例下回歸系數(shù)的模擬分析

某醫(yī)院泌尿外科的臨床醫(yī)師搜集了76例配備有便攜式透析設(shè)備的腎衰病人。記錄了這些病人從開(kāi)始插入導(dǎo)管到感染的時(shí)間(以天計(jì)算)，如果直到研究結(jié)束仍未出現(xiàn)感染，或因?yàn)槟承┰?感染除外)中途移除導(dǎo)管的病人，視為截尾。隨訪結(jié)束時(shí)，共有58例患者出現(xiàn)了感染。此外，還搜集了病人的年齡、性別、疾病類型以及衰弱評(píng)分四項(xiàng)指標(biāo)。這四項(xiàng)指標(biāo)的賦值情況和基本統(tǒng)計(jì)表見(jiàn)表6。利用76例腎衰病人擬合截尾分位數(shù)回歸模型，在

表6 76例腎衰病人生存資料預(yù)后因素及其基本統(tǒng)計(jì)量

α=0.10水平上，利用手動(dòng)向后篩選變量法，入選的變量為性別(sex)、疾病類型GN(disease GN)、疾病類型

AN(disease AN)、疾病類型PKD(disease PKD)以及衰弱評(píng)分(frail)，交互項(xiàng)均無(wú)意義。模型擬合結(jié)果顯示在不同的分位數(shù)函數(shù)上不僅表現(xiàn)為影響變量個(gè)數(shù)的不同，還表現(xiàn)為相同自變量對(duì)應(yīng)的回歸系數(shù)值大小不等。這里，只列舉τ=0.1，0.5，0.9三個(gè)分位數(shù)水平對(duì)應(yīng)的回歸方程。

log t0．1=8.91 －2.73GN －2.273AN+32.183PKD

log t0．5=54.06 －98.35sex －5GN － 69.87AN +15.90PKD－81.61frail

log t0．9=373.24 －96.81GN －141.49AN+156.70PKD－176.82frail

從圖1可以看出性別對(duì)log t影響的回歸系數(shù)隨著τ的變化而變化(先增加后下降)，疾病類型AN(以other為參考)對(duì)log t影響的回歸系數(shù)的絕對(duì)值隨著τ的增加而增加，衰弱評(píng)分對(duì)log t影響回歸系數(shù)的絕對(duì)值隨著τ的增加而增加。即隨著生存時(shí)間的增加，說(shuō)明性別對(duì)生存時(shí)間的保護(hù)作用先逐漸增強(qiáng)后又逐漸減弱，疾病類型AN(以other為參考)和衰弱評(píng)分對(duì)生存時(shí)間的威脅性越來(lái)越強(qiáng)。而從Cox模型擬合的效果來(lái)看，各回歸系數(shù)幾乎不隨τ變化，低估了性別、疾病類型AN和衰弱評(píng)分對(duì)生存時(shí)間的影響作用。從專業(yè)角度看，截尾分位數(shù)回歸擬合的結(jié)果更接近實(shí)際情況。

討 論

在生存分析中，截尾分位數(shù)回歸模型一般用于如下情況:(1)當(dāng)數(shù)據(jù)有離群值，為削弱其對(duì)回歸模型的影響;(2)當(dāng)y的方差不是常數(shù)方差，即存在異方差，或者存在其他類型的異質(zhì)性問(wèn)題。此時(shí)，用中位數(shù)回歸模型估計(jì)給定x時(shí)y的平均水平(中位數(shù))，用其他分位數(shù)回歸模型估計(jì)相應(yīng)的容許區(qū)間或參考值范圍。截尾分位數(shù)回歸模型的特性亦類似于百分位數(shù)。如在截尾分位數(shù)回歸中，中位數(shù)回歸模型較其他百分位數(shù)回歸模型穩(wěn)定，越是接近0%和100%的百分位數(shù)回歸模型越易受離群值和截尾值的影響，越是不穩(wěn)定。因此，在用截尾分位數(shù)回歸模型確定y的容許區(qū)間時(shí)，宜用70%，80%或90%的區(qū)間，而不用95%，98%或99%的區(qū)間〔5〕。

截尾分位數(shù)回歸模型中，回歸系數(shù)向量中^β(τ)第j個(gè)元素表示的是固定其他協(xié)變量時(shí)第j個(gè)協(xié)變量的單位變化引起第τ分位數(shù)的平均變化量。如果在線性截尾分位數(shù)回歸模型中只表現(xiàn)為截距項(xiàng)的不等，即位置的漂移，而協(xié)變量對(duì)應(yīng)的回歸系數(shù)不隨著τ變化，說(shuō)明總體中不存在異質(zhì)性;如果截尾分位數(shù)回歸模型中協(xié)變量的回歸系數(shù)隨著τ變化，說(shuō)明總體中存在異質(zhì)性。因此，分位數(shù)回歸模型不只可用來(lái)做多因素的統(tǒng)計(jì)分析，還能作為一種檢驗(yàn)異質(zhì)性的診斷方法。

1．Portnoy S．Censored regression quantiles．J．Amer．Statist．Assoc，2003，98:1001-1012．

2．Mara Tableman，Jong Sung Kim．Survival Analysis Using S:Analysis of Time-to-event Data．New York:Chapman ＆ Hall/CRC，2004．

3．余紅梅．Cox比例風(fēng)險(xiǎn)回歸模型診斷及預(yù)測(cè)有關(guān)問(wèn)題的研究:〔博士學(xué)位論文〕西安:第四軍醫(yī)大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室，2001．

4．Kalblfeisch JD，Prentice RL．The Statistical Analysis of Failure Time Data．New York:Wiley，1980．

5．季莘，陳峰．百分位數(shù)回歸及其應(yīng)用．中國(guó)衛(wèi)生統(tǒng)計(jì)，1998，15(6):9-11．