張文麗，張 彤，易丹輝＊＊，楊宇飛＊＊

（1.中國人民大學(xué)應(yīng)用統(tǒng)計科學(xué)研究中心，北京 100872；2.中國人民大學(xué)統(tǒng)計學(xué)院 北京 100872；3.中國中醫(yī)科學(xué)院西苑醫(yī)院 北京 100091）

運用Cox模型時打結(jié)數(shù)據(jù)的處理方法探討＊

張文麗1,2，張 彤3，易丹輝1,2＊＊，楊宇飛3＊＊

（1.中國人民大學(xué)應(yīng)用統(tǒng)計科學(xué)研究中心，北京 100872；2.中國人民大學(xué)統(tǒng)計學(xué)院 北京 100872；3.中國中醫(yī)科學(xué)院西苑醫(yī)院 北京 100091）

Cox回歸模型是目前生存分析中最為廣泛使用的方法之一，模型的假設(shè)之一是失效時間不存在打結(jié)情況，即個體之間有著不同的失效時間。在實際應(yīng)用當(dāng)中，生存時間數(shù)據(jù)存在打結(jié)是很常見的。目前有四種常見的處理方法：Exact法，discrete model法，Efron法以及Breslow法。本文研究目的是比較這四種處理方法的優(yōu)劣。本文采用模擬進行比較，設(shè)置了不同的樣本量和打結(jié)程度，比較四種方法在擬合統(tǒng)計量，計算時間，參數(shù)估計精確性等方面的表現(xiàn)，發(fā)現(xiàn)Exact法和discrete model擬合統(tǒng)計量結(jié)果最好，但計算耗時最久；Efron法以及Breslow法運算較快但是在參數(shù)估計方面存在偏差。另外，樣本量和打結(jié)程度也影響處理的結(jié)果，總的來說，當(dāng)結(jié)點數(shù)較小時，四種方法之間差別不大。當(dāng)數(shù)據(jù)量較大或打結(jié)比例較高時，除exact以外的三種近似方法的偏差增加，但運算時間無明顯變化，而exact法的運算時間迅速增加。此時如果估計的精確性沒有估計時間那么重要，Efron法以及Breslow法是不錯的選擇，其中，Efron法更為精確而Breslow方法傾向于低估正確的β值。如果時間上沒有限制，可以選擇Exact法和discretemodel，將得到更為精確的結(jié)果。

生存分析 Cox模型 打結(jié)數(shù)據(jù) 部分似然函數(shù)

Cox回歸模型，或叫相對風(fēng)險模型（Relative Risk Model）是目前生存分析中最為廣泛使用的方法之一。該模型最早由Cox D.R在1972年提出。該模型無需對基準(zhǔn)風(fēng)險函數(shù)做任何的限制，是半?yún)?shù)模型，克服了生存分析傳統(tǒng)的參數(shù)法和非參數(shù)法的局限性，目前廣泛運用于不同治療方法的比較以及各種疾病預(yù)后相關(guān)因素的研究[1]。

Cox模型在應(yīng)用中還存在一些問題，其中之一就是數(shù)據(jù)存在結(jié)點時的處理方法。數(shù)據(jù)存在打結(jié)是指有多個個體有相同的失效時間。在實際應(yīng)用中，由于失效時間往往是以一種離散的方式記錄的，得到的數(shù)據(jù)存在打結(jié)是很普遍的。

失效時間不存在打結(jié)是Cox模型的一個重要假設(shè)，與模型的估計緊密相關(guān)，在該假設(shè)不滿足時仍可以使用Cox模型，但需要對估計的方法進行改進。

目前有四種常見的方式用來處理打結(jié)數(shù)據(jù)，分別是 Exact model、discrete model、Efron 法以及 Breslow法。研究表明，Exact法和discrete model的結(jié)果較為精確，但是計算時間較長。Efron法以及Breslow法用時短，但參數(shù)估計偏差較大。Breslow方法計算比較簡便，是目前大多數(shù)軟件默認(rèn)的處理Cox模型打結(jié)數(shù)據(jù)的方法，但是在R軟件的survival包中，默認(rèn)的方法是Efron法。

本文對四種不同的處理打結(jié)數(shù)據(jù)方法的原理進行討論，并利用模擬討論在不同樣本量水平和打結(jié)程度下，在參數(shù)估計，估計量效率（efficiency of estimators），擬合統(tǒng)計量，計算時間方面的表現(xiàn)。最后結(jié)合具體數(shù)據(jù)進行展示。

1 方法介紹

Cox模型構(gòu)建的思路是，所有研究對象的生存情況是多個影響因素共同作用的結(jié)果，可用風(fēng)險函數(shù)表示，記為 λ(t;x)，其中 x=(x1,x2,…)′是基準(zhǔn)協(xié)變量，在個體進入試驗之前或進入試驗之時已經(jīng)測得，T是絕對連續(xù)的失效時間變量[2]。假設(shè)研究對象中影響生存的因素均不存在，其生存情況用λ0(t)表示，稱為基礎(chǔ)風(fēng)險函數(shù)，那么，研究對象的實際生存情況，即λ(t;x)，是影響因素x在基礎(chǔ)風(fēng)險函數(shù)的基礎(chǔ)上，進一步修改的結(jié)果，即風(fēng)險函數(shù)可表達為：

英國統(tǒng)計學(xué)家D.R.Cox于于1972年首次提出把r(t;x)構(gòu)造為指數(shù)形式，即將風(fēng)險函數(shù)寫作（2.2）式：

其中Z(t)=[Z1(t),…,Zp(t)]′是一個可能時間相依的協(xié)變量向量，是時間t和基準(zhǔn)協(xié)變量的函數(shù)。Cox回歸模型可以估計相對風(fēng)險，例如所建立的回歸模型包含組別變量，該變量是一個二分類變量，以1表示治療組，0表示對照組，則治療組與對照組風(fēng)險函數(shù)的比值為

該比值的含義為，當(dāng)其他變量相同時，在每個時間點，治療組的死亡風(fēng)險都是對照組的eβ倍。

對于定量數(shù)據(jù)的協(xié)變量，如年齡，HR的結(jié)果是：

即該協(xié)變量每增長m單位，HR就要在原來的基礎(chǔ)上乘以exp(mβ)。

Cox模型的假設(shè)之一是失效時間不存在打結(jié)，在該假設(shè)成立的前提下，參數(shù)估計可以通過部分似然（partial likelihood）方法得到[3][4]。

用zl表示第l個個體的解釋變量，排序的失效時間為t1＜…＜tk。用Di={i1,…,idi}表示在ti時刻失效的個體組成的集合，Qi是{i1,…,idi}的di!個排列組成的集合，P=(p1,…pdi)是Qi中的一個元素，用di表示在ti時刻失效的個體數(shù)，用Ri表示在ti時刻的風(fēng)險集。R(tj,P,r)=R(tj)-{p1,…,pr-1}。

根據(jù)Cox（1972）提出的方法，數(shù)據(jù)不存在打結(jié)時，估計β的部分似然的公式為（2.5）式：

利用該方法估計參數(shù)時，似然函數(shù)的計算依賴于事件發(fā)生的順序，個體的失效時間必須是有序的，如果有兩個個體（例如A和B）有著相同的失效時間，在這種情況下，無法確定其中一個（如A）失效時，B是否在該時刻的危險集中。

當(dāng)數(shù)據(jù)存在結(jié)點時，Kalbfleisch and Prentice（2002）[5]提出的exact法考慮在每個存在打結(jié)的時間點上事件發(fā)生的所有可能的排序。在tj處將結(jié)點分解成各種可能的情形后的平均部分似然為：

當(dāng)每一個失效時間點的結(jié)點數(shù)目比較大時，上式的計算量會非常大。此時，可以對似然函數(shù)進行近似。Breslow(1974)[6]提出的近似似然函數(shù)為

Efron方法(1977)[7]提出的似然函數(shù)為

另外，當(dāng)數(shù)據(jù)打結(jié)的比例較大時，可以考慮將失效時間看作離散變量。Cox（1975）建議利用離散Logistic模型（也叫條件logistic模型）：

其中dλ0(t)是一個未指定的離散風(fēng)險函數(shù)，在觀測失效時間點t1＜…＜tk有值，將沒有結(jié)點情形下的部分似然進行推廣，得到以下部分似然函數(shù)。

其中Rdi(ti)是從風(fēng)險集R(ti)中挑出di個個體的所有子集組成的集合，l=(l1,…,ldi)是 Rdi(ti)的一個元素，

Kalbfleisch和 Prentice（2002）指出，Efron法和Breslow法對參數(shù)的估計存在偏差，且的方差估計值是不一致的。模擬顯示，當(dāng)diRi的值較大時，Breslow法對參數(shù)的估計會有較大的偏差。結(jié)點數(shù)很少時，三種方法得到相似的結(jié)果，不存在結(jié)點時，三種方法會得到完全一樣的結(jié)果[5]。

2 數(shù)據(jù)模擬

本文采用模擬比較Exact方法和三種近似方法在參數(shù)估計，估計量效率（efficiency of estimators），擬合統(tǒng)計量和計算時間方面的差異。對于計算時間，使用每種方法重復(fù)10次估計，以得到平均計算時間。

生成兩組數(shù)據(jù)，一組是失效時間服從指數(shù)分布的生存時間數(shù)據(jù)，即風(fēng)險函數(shù)為：

另一組是與該組基礎(chǔ)風(fēng)險函數(shù)相同，HR是e-1的生存時間數(shù)據(jù)，風(fēng)險函數(shù)為

另外，是否刪失的設(shè)置與生存時間獨立，即生成的數(shù)據(jù)包括的變量為：生存時間，組別（0,1），是否刪失（0,1）。

首先在每組1 000個個體的樣本量水平上進行模擬。利用數(shù)據(jù)分組制造結(jié)點，分別將數(shù)據(jù)分到k=50,200,500個時間區(qū)間中去來制造高、中、低三種打結(jié)水平。具體方法是用生存時間數(shù)據(jù)落入的區(qū)間的右端點來代替原先的數(shù)據(jù)以制造結(jié)點。

利用SAS9.4中的PHREGPRocedule擬合模型。

當(dāng)樣本量為每組1 000個個體時，在結(jié)點數(shù)最多即k=50，平均每個時間點有20人的情形下，各種方法的計算時間均小于0.5秒，因此不再對計算時間進行記錄和比較。

表1 n=1 000 β=-1 Cox模型模擬結(jié)果

表2 n=1 000 β=-1 Cox模型模擬結(jié)果

表2展示了四種處理方法處理三種打結(jié)水平的SV和三種擬合統(tǒng)計量結(jié)果。其中SV（standardized measuresof variability）的定義為：

該統(tǒng)計量可以用來衡量參數(shù)估計值的有效性（efficiency of estimators）。

可以看到，Breslow方法在三種打結(jié)程度均有著最高的SV值，在參數(shù)估計值的有效性方面表現(xiàn)較差，其他三種方法之間不存在明顯差異。

在擬合統(tǒng)計量方面，discrete model和Exact法的表現(xiàn)最好，其次是Efron法，Breslow法的結(jié)果最差。

表3 n=100 000 β=-1 Cox模型模擬結(jié)果

表4 n=100 000 β=-1 Cox模型模擬結(jié)果

表5 數(shù)據(jù)打結(jié)情況

在每組100,000個個體的樣本量水平上進行模擬，樣本量變大本身也會導(dǎo)致結(jié)點數(shù)增加，設(shè)置了k=100,500,1 000三個打結(jié)水平。

表4展示了四種處理方法處理三種打結(jié)水平的SV和三種擬合統(tǒng)計量結(jié)果。

Breslow方法在三種打結(jié)程度均有著最高的SV值，在參數(shù)估計值的有效性方面表現(xiàn)較差，其他三種方法之間不存在明顯差異。在擬合統(tǒng)計量方面，the discrete model和Exact法的表現(xiàn)最好，Efron法和Breslow法的結(jié)果幾乎是另外兩種方法的3倍。

4 實際應(yīng)用

某醫(yī)院采用隨機對照臨床研究方法，納入晚期結(jié)直腸癌患者60例，經(jīng)過數(shù)據(jù)處理有效數(shù)據(jù)共53例。其中23人的治療結(jié)局為死亡，30人的治療結(jié)局為未死亡，即有56.6%的數(shù)據(jù)右刪失。兩組均采用常規(guī)治療(營養(yǎng)支持、化療、對癥、中醫(yī))，治療組在此基礎(chǔ)上加用某中藥,對照組加用安慰劑膠囊，治療一段時間后進行隨訪，觀察兩組患者的生存期情況。

數(shù)據(jù)中主要考慮的指標(biāo)包括：組別（治療組=1，對照組=0），性別（女=1，男=0），年齡，患病階段（階段1、階段2、階段3、階段4），是否死亡，OS生存期（單位：月）。其中OS生存期數(shù)據(jù)存在結(jié)點，具體情況如表5所示。這種情況下，運用Cox模型比較兩組的治療，必須考慮打結(jié)數(shù)據(jù)的處理。表6至表10是四種方法處理該數(shù)據(jù)的結(jié)果。

表7是p值結(jié)果。由表可以看出，Breslow得到的p值較大。

表8是參數(shù)估計的標(biāo)準(zhǔn)誤。由表可以看出，Exact法，Breslow法和Efron法的結(jié)果無明顯差異，Discrete法得到的標(biāo)準(zhǔn)誤與其他三種方法相比略大一些。

表9是SV結(jié)果，由表可以看出，Breslow法在四個協(xié)變量上的結(jié)果均是最差的，另外三種方法不存在明顯差異。

表9是擬合統(tǒng)計量結(jié)果，由表可以看出，Exact法結(jié)果最好，discrete法與Exact法差別不大，Breslow法和Efron法結(jié)果較差。

該數(shù)據(jù)有較多結(jié)點，但由于樣本量小，四種方法的計算時間都很短，推薦使用Exact法以得到更精確的結(jié)果，如果軟件（如R）中不包含該方法，則推薦Efron法，該方法估計的結(jié)果與Exact法最接近，而其余兩種方法偏差較大。

5 結(jié)論

樣本量較小的情形下（n≤1000），不同程度的打結(jié)，四種方法的估計時間都很短，此時Breslow法會低估參數(shù)絕對值且偏差較大，Efron法表現(xiàn)好于Breslow法，discrete法會高估參數(shù)絕對值。但還是建議使用Exact法來獲得最為精確的估計，R軟件的survival包沒有包括該方法，可以使用SAS軟件的PHREG PRocedule，選擇“exact”即可。

樣本量較大的情況下（n≥100,000），Exact法的計算時間迅速增加，打結(jié)程度高時個人的電腦可能出現(xiàn)內(nèi)存不足無法利用SAS運算該方法的情形，此時可以考慮使用discrete法，計算時間不到Exact法的一半?？紤]到樣本量較大，各種方法的偏差都小，更推薦使用Breslow法和Efron法，在樣本量為100,000且打結(jié)程度最高時（平均每個時間點有1 000個個體），運算時間仍不超過1秒鐘，其中Efron法更加精確。如果時間上沒有限制，可以使用Exact法或discrete法，可以得到更好的擬合統(tǒng)計量結(jié)果和更準(zhǔn)確的參數(shù)估計結(jié)果。

表6 實際數(shù)據(jù)擬合Cox回歸模型

表6 實際數(shù)據(jù)擬合Cox回歸模型

?

表7 實際數(shù)據(jù)擬合Cox回歸模型p

表8 實際數(shù)據(jù)擬合Cox回歸模型SE

表9 實際數(shù)據(jù)擬合Cox回歸模型SV

1 陳兵,駱福添.生存分析中的回歸模型.中國衛(wèi)生統(tǒng)計,2006,23(5):462-465.

2 金丕煥,陳峰.醫(yī)用統(tǒng)計方法.復(fù)旦大學(xué)出版社,2009:378-385.

3 Cox DR.Regression Models and Life-Tables.JRoy Stat Soc,1972,34(2):187-220.

4 Cox D.R.Partial likelihood.Biometrika,1975,62(2):269-276.

5 Kalbfleisch J D,Prentice R L.Marginal likelihoods based on Cox's regression and lifemodel.Biometrika,1973,60:267-279.

6 Breslow N.Covariance analysis of censored survival data,Biometrics,1974,30:89-99.

7 Efron B.The efficiency of Cox's likelihood function for censored data.J Am Stat Assoc,1977,72,557-565.

8 Hertz P I,Rockhill B.Validity and Efficiency of Approximation Methodsfor Tied Survival Timesin Cox Regression.Biometrics,1997,53(3):1151-1156.

9 Borucka J.Methods of Handling Tied Events in the Cox Proportional Hazard Model.Ieee,2014,2(2):92-106.

Discussion on Methodsfor Tied Survival Timesin Cox Model

Zhang Wenli1,2,Zhang Tong3,Yi Danhui1,2,Yang Yufei3
(1.Center for Applied Statisticsof Renmin University of China,Beijing 100872,China;2.School of Statistics,Renmin University of China,Beijing 100872,China;3.Xiyuan Hospital of China Academy of Chinese Medical Sciences,Beijing 100091,China)

Cox regression model is one of the most widely used methods in the survival analysis.One assumption of this model is that there is no tie in the failure times,that is,individual has different failure times.In practical applications,the existence of ties in time data is very common.In this paper,four common methods of dealing with ties in Cox model,including Exact method,discrete model method,Efron method and Breslow method,were compared with simulation.The results showed that Exact method and discrete model were the best,but they took the longest time.Efron method and Breslow method were faster but there was a greater deviation in parameter estimation.Moreover,the sample amount and ties degree also affect the results.In general,when there are a few ties,the difference between four methods was small;and in the case of large datasets or a large number of ties,the biasof three approximation methodsincreased except Exact method.However,there was no significant change on computational time.While the computational time of the Exact method increased rapidly.Therefore,if the estimation precision is not as important as the estimation time,Efron method and Breslow method will be good choices.Efron method is more preferably as it is more precise.And Breslow method tends to underestimate the trueβ.If there is no limit in time,Exact method and discrete model can be chosen to achieve more accurate results.

Survival analysis,Cox model,tied data,partial likelihood function

10.11842/wst.2017.09.007

R33

A

2017-05-18

修回日期：2017-08-23

＊ 中國人民大學(xué)2017年度‘中央高校建設(shè)世界一流大學(xué)（學(xué)科）和特色發(fā)展引導(dǎo)專項資金’，負(fù)責(zé)人：易丹輝；和教育部人文社會科學(xué)重點研究基地重大項目(16JJD910002)：基于大數(shù)據(jù)的精準(zhǔn)醫(yī)學(xué)生物統(tǒng)計分析方法及其應(yīng)用研究，負(fù)責(zé)人：？？。

＊＊ 通訊作者：易丹輝，中國人民大學(xué)教授，博士生導(dǎo)師，主要研究方向：風(fēng)險管理與保險、預(yù)測與決策。楊宇飛，博士生導(dǎo)師，中國中醫(yī)科學(xué)院西苑醫(yī)院腫瘤診治部主任、主任醫(yī)師，主要研究方向：中西醫(yī)結(jié)合癌癥治療。

（責(zé)任編輯：張娜娜，責(zé)任譯審：王 晶）