王 娟 王 彤

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

原理與方法

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

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

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

模擬分析

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

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

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

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

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

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

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

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

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

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

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

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

3．違背比例風險假定的模擬分析

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

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

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

實例分析

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

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

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

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

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

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

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

討 論

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

截尾分位數(shù)回歸模型中，回歸系數(shù)向量中^β(τ)第j個元素表示的是固定其他協(xié)變量時第j個協(xié)變量的單位變化引起第τ分位數(shù)的平均變化量。如果在線性截尾分位數(shù)回歸模型中只表現(xiàn)為截距項的不等，即位置的漂移，而協(xié)變量對應的回歸系數(shù)不隨著τ變化，說明總體中不存在異質(zhì)性;如果截尾分位數(shù)回歸模型中協(xié)變量的回歸系數(shù)隨著τ變化，說明總體中存在異質(zhì)性。因此，分位數(shù)回歸模型不只可用來做多因素的統(tǒng)計分析，還能作為一種檢驗異質(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比例風險回歸模型診斷及預測有關問題的研究:〔博士學位論文〕西安:第四軍醫(yī)大學衛(wèi)生統(tǒng)計學教研室，2001．

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

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