擴(kuò)展Cox模型在非線性生存資料分析中的預(yù)測(cè)能力比較

陳雨軒，韋紅霞，潘建紅，安勝利

1南方醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計(jì)學(xué)系，廣東 廣州 510515；2國(guó)家藥品監(jiān)督管理局藥品審評(píng)中心，北京100022

在醫(yī)學(xué)研究中，研究者們經(jīng)常通過(guò)構(gòu)建回歸模型來(lái)分析變量間的關(guān)系。通?；貧w模型的使用前提假設(shè)是自變量和因變量間或因變量的某種函數(shù)間呈線性關(guān)聯(lián)，但在實(shí)際應(yīng)用中這個(gè)條件并不一定滿足。常見(jiàn)的解決方法是將連續(xù)型變量分類，但分類數(shù)目和節(jié)點(diǎn)位置的選擇往往帶有主觀性，并且分類會(huì)損失信息。因此，一個(gè)更好的解決方法是擬合自變量與因變量之間的非線性關(guān)系。

以往生存分析中對(duì)非線性數(shù)據(jù)的處理通常有兩種，一個(gè)是直接擬合機(jī)器學(xué)習(xí)模型；另一個(gè)就是把連續(xù)型的變量轉(zhuǎn)換成分類變量并以啞變量形式納入傳統(tǒng)生存分析模型。但是后者可能導(dǎo)致協(xié)變量不滿足比例風(fēng)險(xiǎn)假定［1］。為了彌補(bǔ)傳統(tǒng)生存分析模型的應(yīng)用條件限制及缺陷，避免對(duì)變量關(guān)系間的形狀進(jìn)行參數(shù)約束，目前常用平滑工具來(lái)解決這類問(wèn)題。平滑工具中限制性立方樣條（RCS）是目前分析非線性關(guān)系的最常見(jiàn)的方法之一［2-4］。生存分析中，限制性立方樣條結(jié)合Cox模型，不僅可以用于探索非線性關(guān)系，越來(lái)越多的學(xué)者也將其運(yùn)用于構(gòu)建預(yù)測(cè)模型［5-8］。

近年來(lái)機(jī)器學(xué)習(xí)以及深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)的發(fā)展，使得更多的方法可以用于解決復(fù)雜生存數(shù)據(jù)構(gòu)建預(yù)測(cè)模型。深度學(xué)習(xí)屬于機(jī)器學(xué)習(xí)的一個(gè)子域，機(jī)器學(xué)習(xí)方法中，隨機(jī)生存森林方法是最常見(jiàn)的分析方法，當(dāng)有新的機(jī)器學(xué)習(xí)方法提出時(shí)一般會(huì)跟隨機(jī)生存森林做比較，所以本文也考慮納入隨機(jī)生存森林作為參照方法。機(jī)器學(xué)習(xí)生存分析方法，包括隨機(jī)生存森林、條件推斷森林等［9-13］；深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)方法如DeepSurv，DeepHit，Neural net-extended time-dependent cox model等［14］，都可以解決非線性、交互等復(fù)雜情況，其在高維數(shù)據(jù)中有明顯的優(yōu)勢(shì)，并且不需要滿足傳統(tǒng)Cox回歸模型比例風(fēng)險(xiǎn)假定等前提條件，但是機(jī)器學(xué)習(xí)以及神經(jīng)網(wǎng)絡(luò)模型實(shí)施起來(lái)較為復(fù)雜，特別是深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)需要用Python軟件來(lái)實(shí)現(xiàn)，其所花費(fèi)的時(shí)間比傳統(tǒng)模型要多。直接使用機(jī)器學(xué)習(xí)或者深度學(xué)習(xí)方法來(lái)分析，花費(fèi)人力物力的同時(shí)其擬合效果也不一定比傳統(tǒng)生存分析方法好。所以在非線性關(guān)系存在的生存數(shù)據(jù)分析中，各方法的優(yōu)劣還有待研究比較。

1 方法

方法部分主要介紹Cox比例風(fēng)險(xiǎn)回歸模型（Cox）、限制性立方樣條Cox回歸（Cox_RCS）、深度生存神經(jīng)網(wǎng)絡(luò)（Cox_DNN）這3種方法。

1.1 Cox比例風(fēng)險(xiǎn)回歸模型

Cox比例風(fēng)險(xiǎn)回歸模型是1972年由英國(guó)統(tǒng)計(jì)學(xué)家D.R.Cox提出的一種半?yún)?shù)模型［15］，其數(shù)學(xué)公式為：

其中λ(t,Xi)為受試者在一組預(yù)測(cè)變量Xi條件下所預(yù)測(cè)的t時(shí)刻的事件發(fā)生風(fēng)險(xiǎn)，λ0(t)是基線風(fēng)險(xiǎn)函數(shù)，β為預(yù)測(cè)變量的回歸系數(shù)。醫(yī)學(xué)研究人員常使用Cox模型評(píng)估預(yù)后協(xié)變量在死亡或疾病復(fù)發(fā)等事件中的重要性，并隨后告知患者其治療選擇；Cox 比例風(fēng)險(xiǎn)回歸模型使用時(shí)需要滿足兩個(gè)基本假定：（1）比例風(fēng)險(xiǎn)假定：λ(t,Xi)/λ0(t)為固定值，即協(xié)變量對(duì)生存率的影響不隨時(shí)間的改變而改變；（2）對(duì)數(shù)線性假定：對(duì)數(shù)風(fēng)險(xiǎn)比應(yīng)與模型中的連續(xù)型協(xié)變量呈線性關(guān)系。模型中的連續(xù)型變量可能會(huì)對(duì)生存時(shí)間產(chǎn)生非線性風(fēng)險(xiǎn)，如果不考慮這個(gè)風(fēng)險(xiǎn)，可能會(huì)導(dǎo)致結(jié)果偏差。

1.2 限制性立方樣條Cox回歸

限制性立方樣條是最常用來(lái)檢驗(yàn)假設(shè)關(guān)系不是線性的一種方法，或者用來(lái)總結(jié)非線性關(guān)系的一種方法。限制性立方樣條函數(shù)只是一個(gè)自變量的變換。因此，它們不僅可以用于普通最小二乘回歸，還可以用于邏輯回歸、生存分析等。

臨床評(píng)估科學(xué)研究所進(jìn)行的兩項(xiàng)研究很好地說(shuō)明了限制性立方樣條的兩種用途。一是識(shí)別非線性關(guān)系，二是構(gòu)建預(yù)測(cè)模型。如Κendzerka等［5］使用限制性立方樣條來(lái)檢驗(yàn)連續(xù)型預(yù)測(cè)變量與心血管事件相關(guān)呼吸暫?；颊咦≡猴L(fēng)險(xiǎn)之間的線性關(guān)系，當(dāng)發(fā)現(xiàn)非線性的關(guān)聯(lián)時(shí)，通過(guò)RCS構(gòu)建預(yù)測(cè)模型實(shí)行風(fēng)險(xiǎn)預(yù)測(cè)。正如前文所述，連續(xù)型變量進(jìn)行分類可能會(huì)丟失重要的信息，因此限制性立方樣條為分析連續(xù)型預(yù)測(cè)變量對(duì)預(yù)測(cè)結(jié)果的影響提供了一個(gè)有用的工具，它在預(yù)測(cè)變量和結(jié)局之間的關(guān)系形式上允許很大的靈活性。

限制性立方樣條Cox 模型，1989 年由Durrleman等［16］學(xué)者提出，用于生存分析中擬合靈活的非線性關(guān)系。在Cox比例風(fēng)險(xiǎn)回歸模型中，考慮變量的非線性效應(yīng)后的模型為：

1.3 深度生存神經(jīng)網(wǎng)絡(luò)Cox模型

DeepSurv是一種深度前饋神經(jīng)網(wǎng)絡(luò)［17］，是類似于Faraggi-Simon網(wǎng)絡(luò)的多層感知器。不同的是加入了多個(gè)隱藏層，以及其他新的技術(shù)，例如權(quán)重衰減正則化、整流線性單位（ReLU）、批處理歸一化（Bacth normalization）、dropout、隨機(jī)梯度下降使用Nesterov動(dòng)量、梯度修剪和學(xué)習(xí)率調(diào)度等。它通過(guò)網(wǎng)絡(luò)θ的權(quán)重參數(shù)來(lái)預(yù)測(cè)患者的協(xié)變量對(duì)其風(fēng)險(xiǎn)率的影響。圖1 是DeepSurv 的網(wǎng)絡(luò)結(jié)構(gòu)，X為輸入到網(wǎng)絡(luò)的患者的基線數(shù)據(jù)。DeepSurv模型不涉及篩選特征變量，網(wǎng)絡(luò)的隱藏層使用了一個(gè)全連接的節(jié)點(diǎn)層（Fully-Connected Layer）和一個(gè)dropout 層，交替堆疊。對(duì)最后一個(gè)dropout 層的輸出進(jìn)行線性組合，最終得到風(fēng)險(xiǎn)率,其估計(jì)了Cox 模型中的對(duì)數(shù)風(fēng)險(xiǎn)函數(shù)h(x),通過(guò)設(shè)置目標(biāo)函數(shù)的平均偏負(fù)對(duì)數(shù)似然來(lái)訓(xùn)練網(wǎng)絡(luò)，其重設(shè)計(jì)的損失函數(shù)如下：

圖1 DeepSurv網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.1 Diagram of the structure of DeepSurv.

其中，NE=1是觀察到結(jié)局的病人例數(shù)，E為生存結(jié)局，T為生存時(shí)間，x為基線特征變量，λ是L2 正則化參數(shù)。

在使用深度神經(jīng)網(wǎng)絡(luò)時(shí)，無(wú)論網(wǎng)絡(luò)有多少層，每一層節(jié)點(diǎn)的輸入都是上層網(wǎng)絡(luò)輸出的線性組合，需要引入非線性激活函數(shù)來(lái)擬合非線性關(guān)系來(lái)提高模型的表達(dá)能力，目前為止激活函數(shù)種類已超過(guò)20種，本文所用的激活函數(shù)有Sigmoid、Tanh、ReLU和LeakyReLU，它們的數(shù)學(xué)形式分別為,本文中a=0.01 。激活函數(shù)是向神經(jīng)網(wǎng)絡(luò)中引入非線性因素，通過(guò)激活函數(shù)神經(jīng)網(wǎng)絡(luò)就可以擬合各種曲線。激活函數(shù)主要分為飽和激活函數(shù)（Saturated Neurons）和非飽和函數(shù)（Onesided Saturations）。Sigmoid和Tanh是飽和激活函數(shù)，而ReLU以及其變種LeakyReLU為非飽和激活函數(shù)。非飽和激活函數(shù)主要有如下優(yōu)勢(shì)：（1）可以解決梯度消失問(wèn)題；（2）可以加速收斂。Sigmoid極容易導(dǎo)致梯度消失問(wèn)題。假設(shè)神經(jīng)元輸入Sigmoid的值特別大或特別小，對(duì)應(yīng)的梯度約等于0，即使從上一步傳導(dǎo)來(lái)的梯度較大，該神經(jīng)元權(quán)重(w)和偏置(bias)的梯度也會(huì)趨近于0，計(jì)算費(fèi)時(shí)并且導(dǎo)致參數(shù)無(wú)法得到有效更新。ReLU激活函數(shù)的提出就是為了解決梯度消失問(wèn)題；但ReLU會(huì)存在神經(jīng)元“死亡”問(wèn)題。而LeakyReLU的提出解決了ReLU的神經(jīng)元“死亡”問(wèn)題。

2 模擬與實(shí)例

2.1 模擬

本研究中的Cox 回歸使用R 軟件“survival”包實(shí)現(xiàn)。限制性立方樣條Cox 回歸的擬合使用“rms”包完成。隨機(jī)生存森林（RSF）實(shí)現(xiàn)使用“randomForestSRC”包完成。以上均基于R 軟件的4.0.3版本。DeepSurv（Cox_DNN）深度神經(jīng)網(wǎng)絡(luò)的擬合應(yīng)用Python 軟件3.7.0 版本完成，具體使用的是Pysurvival 庫(kù)中的NonLinearCoxPHModel 函數(shù)。

2.1.1 生存數(shù)據(jù)模擬設(shè)置 在模擬研究中，生成生存時(shí)間常使用參數(shù)模型，常見(jiàn)的有參數(shù)PH模型、AFT模型。參數(shù)PH 模型，其參數(shù)假定的分布包括指數(shù)分布、Weibull分布和Gompertz分布等［18,19］。產(chǎn)生生存數(shù)據(jù)時(shí)，指數(shù)分布是最常見(jiàn)的分布假設(shè)，但實(shí)際數(shù)據(jù)往往并不滿足指數(shù)分布，如醫(yī)藥領(lǐng)域常用Gompertz分布來(lái)描述數(shù)據(jù)［20］。因此本研究用Gompertz比例風(fēng)險(xiǎn)模型來(lái)產(chǎn)生滿足比例風(fēng)險(xiǎn)假設(shè)的生存數(shù)據(jù)。同時(shí)本研究也將利用Weibull加速失效時(shí)間模型（AFT）產(chǎn)生不滿足比例風(fēng)險(xiǎn)假設(shè)的生存數(shù)據(jù)［21,22］。使用Gaussian非線性關(guān)聯(lián)，其公式如下：

其中a表示Gaussian曲線尖峰高度，本文模擬時(shí)取值為1；b表示Gaussian曲線中心位置，本文模擬時(shí)取值為0；γ表示Gaussian曲線寬度，本文模擬時(shí)取值為0.5。

模擬一：

基于Gompertz比例風(fēng)險(xiǎn)模型生成生存時(shí)間，假設(shè)生存時(shí)間與特征變量之間的關(guān)系如下：

其中α,λ是分布參數(shù)。U為服從[0,1]的均勻分布的隨機(jī)數(shù)。時(shí)間T的生存函數(shù)可表示為：

本研究使用Gaussian非線性Gompertz比例風(fēng)險(xiǎn)模型來(lái)產(chǎn)生滿足比例風(fēng)險(xiǎn)假設(shè)的生存數(shù)據(jù)，公式如下：

綜上，通過(guò)以下過(guò)程產(chǎn)生Gaussian 非線性Gompertz PH生存數(shù)據(jù)：（1）先確定樣本量n，并根據(jù)生存分布的假設(shè)，設(shè)定好生存分布參數(shù)取值，Gompertz比例風(fēng)險(xiǎn)模型的分布參數(shù)α,λ，模擬時(shí)分別取0.5，1；（2）通過(guò)軟件產(chǎn)生服從[0,1]均勻分布的隨機(jī)數(shù)U和服從[-1,1]均勻分布的協(xié)變量X；（3）設(shè)定好回歸系數(shù)β的取值，線性變量的系數(shù)β取值都為1，非線性變量的系數(shù)β#在模擬時(shí)分別取值為1，3，6；根據(jù)公式5去模擬產(chǎn)生生存時(shí)間Ti；（4）根據(jù)事先指定的刪失時(shí)間變量L的分布，以及設(shè)置的刪失比例F%，本文設(shè)置的刪失比例分別有10%，20%，40%，60%，80%；通過(guò)迭代確定不同刪失分布的參數(shù)，生成刪失時(shí)間的n個(gè)時(shí)間點(diǎn)Li；（5）通過(guò)比較Ti、Li的大小，得到每個(gè)個(gè)體的截尾指示變量δi=Ι[Ti≤Li]，其中Ι[·] 為示性函數(shù)。若Ti≤Li，則δi=1，表示終點(diǎn)事件發(fā)生，否則δi=0，表示刪失。生存時(shí)間ti=min(Ti,Li)。最終得到含隨機(jī)刪失的生存數(shù)據(jù)（ti,δi）。

模擬二：

基于加速失效時(shí)間模型(AFT)生成生存時(shí)間，假設(shè)對(duì)數(shù)生存時(shí)間與特征變量之間的關(guān)系如下：

其中σ為尺度參數(shù)，εi為隨機(jī)誤差。對(duì)于對(duì)數(shù)線性生存模型來(lái)說(shuō)，時(shí)間T的生存函數(shù)可表示為εi的生存函數(shù)：

若指定隨機(jī)誤差項(xiàng)εi的分布，則稱為參數(shù)加速失效時(shí)間模型。

設(shè)置誤差項(xiàng)服從極值分布，即f(ε)=exp(ε-exp(ε))，尺度參數(shù)σ≠1 時(shí)，對(duì)應(yīng)為Weibull回歸模型。本研究使用Weibull加速失效模型來(lái)產(chǎn)生不滿足比例風(fēng)險(xiǎn)的Gaussian非線性生存數(shù)據(jù)，公式如下：

與前述一樣，通過(guò)以下過(guò)程產(chǎn)生Gaussian非線性Weibull AFT 生存數(shù)據(jù)：（1）先確定樣本量n，以及Weibull AFT模型的分布參數(shù)σ，取值為0.1；（2）通過(guò)軟件產(chǎn)生服從[0,1]均勻分布的隨機(jī)數(shù)ε和服從[-1,1]均勻分布的協(xié)變量X；（3）設(shè)定好回歸系數(shù)β的取值，β#在模擬時(shí)分別取值為1，3，6，其余β值都為1；根據(jù)公式10去模擬產(chǎn)生生存時(shí)間Ti；第4和5步與前述Gaussian非線性Gompertz PH生存數(shù)據(jù)產(chǎn)生步驟一致。

綜上，模擬數(shù)據(jù)主要分為兩類，一類是滿足PH假定的數(shù)據(jù)集（模擬一），另一類是不滿足PH假定的數(shù)據(jù)集（模擬二）。模擬總共設(shè)置10個(gè)變量，其中與生存結(jié)局有關(guān)的非線性變量1個(gè)（x1），線性變量1個(gè)（x2），與生存結(jié)局無(wú)關(guān)的線性變量8個(gè)，樣本量分別設(shè)置為200，500，1000。

2.1.2 模型參數(shù)設(shè)置 模擬時(shí)，限制性立方樣條節(jié)點(diǎn)數(shù)統(tǒng)一采用3節(jié)點(diǎn)。

對(duì)于神經(jīng)網(wǎng)絡(luò)參數(shù)的設(shè)置，由于基線特征變量的數(shù)量不多，因此在尋找神經(jīng)網(wǎng)絡(luò)最優(yōu)參數(shù)時(shí)，各個(gè)參數(shù)范圍分別為隱藏層神經(jīng)元個(gè)數(shù)：3、5、7、10、15個(gè)；隱藏層層數(shù)：1 到5 層；激活函數(shù)類型：Sigmoid、Tanh、ReLU 和LeakyReLU；學(xué)習(xí)率：1e-5,1e-4,1e-3,1e-2,1e-1。使用dropout和正則化技巧避免模型的過(guò)擬合，設(shè)定dropout比率為0.5，L2 正則化參數(shù)為1e-4；默認(rèn)使用Adaptive Momentum（adam）優(yōu)化器。

在隨機(jī)生存森林中，有兩個(gè)重要的參數(shù)，分別是節(jié)點(diǎn)預(yù)選的變量個(gè)數(shù)（通常默認(rèn)最佳的變量個(gè)數(shù)為全部變量個(gè)數(shù)的平方根）和隨機(jī)生存森林中樹(shù)的數(shù)量。通常來(lái)說(shuō)，隨機(jī)森林中樹(shù)的數(shù)目決定了整個(gè)隨機(jī)生存森林運(yùn)行效果和時(shí)間。在擬合隨機(jī)生存森林前，必須先確定生存樹(shù)的個(gè)數(shù)，本文中當(dāng)生存樹(shù)的數(shù)量大于500時(shí)，錯(cuò)誤率已趨于穩(wěn)定，所以在本文的所有模擬研究中，隨機(jī)生存森林模型的生存樹(shù)棵數(shù)定為500，節(jié)點(diǎn)預(yù)選變量數(shù)為全部變量的平方根，分裂準(zhǔn)則為log-rank。

在預(yù)測(cè)能力比較的模擬研究中，蒙特卡洛模擬的次數(shù)設(shè)置為1000次，產(chǎn)生數(shù)據(jù)總樣本量設(shè)置為200、500、1000。刪失率設(shè)置為10%、20%、40%、60%。為了使得四種方法的預(yù)測(cè)準(zhǔn)確度具有可比性并避免過(guò)擬合現(xiàn)象，對(duì)每次模擬產(chǎn)生的數(shù)據(jù)集，我們隨機(jī)選取其中的70%作為訓(xùn)練集擬合模型，剩余的30%作為測(cè)試集，得到模型的預(yù)測(cè)能力評(píng)價(jià)指標(biāo)一致性指數(shù)（C-index）和積分布萊爾評(píng)分（Integrated Brier Score,IBS），最后用所有模擬結(jié)果的C-index平均值和IBS平均值分別作為評(píng)估指標(biāo)。

2.1.3 模型評(píng)估 在模型預(yù)測(cè)能力評(píng)估中，使用C-index和IBS這兩個(gè)指標(biāo)來(lái)評(píng)價(jià)：其中C-index為主要評(píng)價(jià)指標(biāo)用于評(píng)價(jià)模型的預(yù)測(cè)區(qū)分度，IBS為次要評(píng)價(jià)指標(biāo)用于評(píng)價(jià)模型的預(yù)測(cè)校準(zhǔn)度。對(duì)于一個(gè)疾病預(yù)測(cè)模型，應(yīng)先考慮區(qū)分度，如果模型區(qū)分度較差，不能區(qū)分不同風(fēng)險(xiǎn)人群，那么此模型就失去臨床應(yīng)用價(jià)值，再繼續(xù)評(píng)價(jià)校準(zhǔn)度也無(wú)意義了。

一致性指數(shù)(C-index)最早是由范德堡大學(xué)生物統(tǒng)計(jì)教授Frank E Harrell Jr 1996年提出，主要用于計(jì)算生存分析中的Cox模型預(yù)測(cè)值與真實(shí)之間的區(qū)分度，常用在評(píng)價(jià)患者預(yù)后模型的預(yù)測(cè)精度。C-index的取值范圍為[0.5,1]，值越接近于1，表明模型預(yù)測(cè)準(zhǔn)確度越高。Cindex=0.5，表明模型預(yù)測(cè)為隨機(jī)預(yù)測(cè)，預(yù)測(cè)能力低。可以用以下公式計(jì)算［6,23］：

布萊爾評(píng)分用于評(píng)估在給定時(shí)間t的預(yù)測(cè)生存函數(shù)的準(zhǔn)確性,它表示觀察到的生存狀態(tài)和預(yù)測(cè)的生存概率之間的平均平方距離［24］，假設(shè)樣本量為N，?i∈[1,N],(xi,ti)分別為第i個(gè)個(gè)體的預(yù)測(cè)變量和生存時(shí)間，為預(yù)測(cè)的生存函數(shù)，此時(shí)布萊爾評(píng)分計(jì)算公式為：

當(dāng)存在刪失時(shí)，需要使用逆概率刪失加權(quán)法對(duì)平方距離進(jìn)行加權(quán)來(lái)調(diào)整得分，調(diào)整后計(jì)算公式為：

IBS 提供了在所有可用時(shí)間的模型性能的總體計(jì)算。

2.2 模擬結(jié)果

2.2.1 滿足PH假定的情況 圖2展示了變量滿足比例風(fēng)險(xiǎn)假定時(shí)，不同刪失率（設(shè)置為0.1，0.2，0.4和0.6）和樣本量（設(shè)置為200，500，1000）情況下4種模型的預(yù)測(cè)區(qū)分度指標(biāo)C-index和預(yù)測(cè)校準(zhǔn)度指標(biāo)IBS的變化情況，C-index值越高，模型的預(yù)測(cè)區(qū)分度就越高，模型的預(yù)測(cè)能力越好。圖2的右縱坐標(biāo)為生存數(shù)據(jù)產(chǎn)生時(shí)非線性變量的系數(shù)值（β#）及其余線性變量的系數(shù)取值均為1。其中Cox為Cox比例風(fēng)險(xiǎn)回歸模型，Cox_DNN為深度生存神經(jīng)網(wǎng)絡(luò)Cox模型，Cox_RCS為限制性立方樣條Cox模型，RSF為隨機(jī)生存森林模型。在所有模擬條件下（樣本量200-1000，刪失率10%～60%），刪失率和樣本量對(duì)Cox_RCS模型預(yù)測(cè)區(qū)分度的穩(wěn)定性影響不大，此時(shí)的Cox_RCS 模型預(yù)測(cè)區(qū)分度始終最好的；Cox_DNN 受刪失率的影響較大，當(dāng)刪失率＞40%時(shí)，Cox_DNN的預(yù)測(cè)區(qū)分度大幅降低；但當(dāng)樣本量較大（本文≥1000）刪失率不高（＜40%）時(shí)，Cox_DNN 與Cox_RCS的C-index接近；Cox_DNN在樣本量足夠大、刪失率越低時(shí)預(yù)測(cè)能力會(huì)越高。RSF的預(yù)測(cè)區(qū)分度略遜于其他模型。

圖2 4種方法在非線性變量系數(shù)β=1且滿足PH模擬數(shù)據(jù)集中的C-index(上)以及IBS（下）Fig.2 Concordance index (top) and integrated Brier score (bottom) of Cox,Cox_DNN,Cox_RCS and RSF for PH datasets when the nonlinear variable coefficient β=1.

IBS值越低，模型的預(yù)測(cè)校準(zhǔn)度就越好。IBS的結(jié)果與C-index一致；刪失率和樣本量對(duì)Cox_RCS模型預(yù)測(cè)校準(zhǔn)度的穩(wěn)定性影響不大，此時(shí)的Cox_RCS模型預(yù)測(cè)校準(zhǔn)度始終是最好的；RSF的預(yù)測(cè)校準(zhǔn)度表現(xiàn)不佳，Cox_DNN 受刪失率的影響大，刪失率大于40%時(shí)Cox_DNN的預(yù)測(cè)校準(zhǔn)度大幅降低。

總之，在滿足PH假定的數(shù)據(jù)中，在預(yù)測(cè)區(qū)分度和測(cè)校準(zhǔn)度上，Cox_RCS模型預(yù)測(cè)能力表現(xiàn)最好；在樣本量較高（≥1000）、刪失率較低（＜40%）時(shí)，Cox_DNN模型預(yù)測(cè)能力與Cox_RCS相當(dāng)，RSF預(yù)測(cè)表現(xiàn)不佳。

2.2.2 不滿足PH假定的情況 圖3展示了變量不滿足比例風(fēng)險(xiǎn)假定時(shí)，不同刪失率（設(shè)置為0.1，0.2，0.4 和0.6）和樣本量（設(shè)置為200，500，1000）情況下4種模型的預(yù)測(cè)區(qū)分度指標(biāo)C-index和預(yù)測(cè)校準(zhǔn)度指標(biāo)IBS的變化情況，C-index值越高，IBS值越低，則模型的預(yù)測(cè)能力就越好。圖3的右縱坐標(biāo)為生存數(shù)據(jù)產(chǎn)生時(shí)非線性變量的系數(shù)值（β#）及其余線性變量的系數(shù)取值均為1。

圖3 4種方法在非線性變量系數(shù)β=1且不滿足PH模擬數(shù)據(jù)集中的C-index（上）和IBS（下）Fig.3 Concordance index(top)and integrated Brier score (bottom) of Cox,Cox_DNN,Cox_RCS and RSF for non-PH datasets when the nonlinear variable coefficient β=1.

在預(yù)測(cè)區(qū)分度上，Cox模型的預(yù)測(cè)區(qū)分度低于其他3種方法，表明Cox模型在不滿足比例風(fēng)險(xiǎn)的數(shù)據(jù)中擬合效果不佳；樣本量為200，刪失率在0.1～0.6 時(shí)，Cox_RCS的C-index值在0.9左右波動(dòng)，其余3個(gè)模型的值都在0.85左右波動(dòng)，Cox_RCS的預(yù)測(cè)區(qū)分度高于其他3個(gè)模型；但是隨著樣本量增大，Cox_DNN的Cindex 值也在隨之增大；當(dāng)樣本量大于等于500 時(shí)，Cox_DNN的預(yù)測(cè)區(qū)分度最高；當(dāng)樣本量達(dá)到1000時(shí)，無(wú)論刪失率變化如何，Cox_DNN的C-index值能保持在0.9以上；當(dāng)樣本量低于1000時(shí)，隨著刪失率增大，Cox_DNN的C-index值在逐漸降低。這也說(shuō)明了在不滿足比例風(fēng)險(xiǎn)假定的數(shù)據(jù)集中，模型擬合的不穩(wěn)定性。

在預(yù)測(cè)校準(zhǔn)度上，Cox_RCS模型的預(yù)測(cè)校準(zhǔn)度高于其他3種方法，Cox模型的預(yù)測(cè)校準(zhǔn)度最低；隨著樣本量的增加，Cox_DNN的預(yù)測(cè)校準(zhǔn)度也在逐漸升高，刪失率的增大會(huì)使Cox_DNN的預(yù)測(cè)校準(zhǔn)度降低，但刪失率對(duì)另外3個(gè)模型的預(yù)測(cè)校準(zhǔn)度影響不大；RSF的預(yù)測(cè)校準(zhǔn)度高于Cox模型，低于其他兩種模型。

總之，在不滿足比例風(fēng)險(xiǎn)的數(shù)據(jù)中，Cox_RCS在模型預(yù)測(cè)校準(zhǔn)度上占優(yōu)但在預(yù)測(cè)區(qū)分度上表現(xiàn)不穩(wěn)定。Cox_DNN的預(yù)測(cè)校準(zhǔn)度受刪失率和樣本量的影響大；樣本量越大（本文中≥500）、刪失率越低（本文中＜40%），則Cox_DNN的預(yù)測(cè)校準(zhǔn)度越高。Cox和RSF的預(yù)測(cè)能力都表現(xiàn)不佳。

2.3 實(shí)例

通過(guò)實(shí)例數(shù)據(jù)（WHAS500）比較四種方法（Cox,Cox_RCS,Cox_DNN和RSF）的優(yōu)劣。為防止過(guò)擬合現(xiàn)象發(fā)生，將隨機(jī)選取數(shù)據(jù)集的70%作為訓(xùn)練集訓(xùn)練模型，剩余的30%作為測(cè)試集，用C-index和IBS評(píng)價(jià)模型在測(cè)試集中的預(yù)測(cè)區(qū)分度和校準(zhǔn)度，對(duì)該過(guò)程重復(fù)1000次，采用箱圖展示預(yù)測(cè)指標(biāo)C-index和IBS的結(jié)果（表1）。

表1 實(shí)例數(shù)據(jù)集中神經(jīng)網(wǎng)絡(luò)參數(shù)調(diào)節(jié)范圍Tab.1 Adjustment range of neural network parameters in the example dataset

伍斯特心臟病研究數(shù)據(jù)集（WHAS）包含了500名經(jīng)歷過(guò)急性心肌梗死患者，是右刪失生存數(shù)據(jù)［25］，WHAS500 數(shù)據(jù)集下載網(wǎng)址：https://github.com/rfcooper/whas/blob/master/whas500.csv。500名患者在觀察時(shí)間內(nèi)死亡215人，其余患者刪失，刪失率為57%。生存時(shí)間分布如圖4所示，中位生存時(shí)間為1627 d（生存時(shí)間范圍：1～2358 d）。表2列出了該數(shù)據(jù)集生存時(shí)間lenfol和刪失變量fstat和14個(gè)基線變量的基本信息。

圖4 數(shù)據(jù)WHAS500的生存曲線圖Fig.4 Survival curves of WHAS500.

表2 WHAS數(shù)據(jù)集的基線指標(biāo)信息分布Tab.2 Characteristics of the covariates in Dataset WHAS[Mean±SD(min～max)or n(%)]

變量選擇通過(guò)RSF最小深度法變量重要性篩選得出，刪除分布嚴(yán)重不均衡的變量sho、av3和重要性排名最末的變量afb和cvd，最終選擇納入模型的變量有10個(gè)age,hr,sysbp,diasbp,bmi,los,gender,chf,miord,mitype。用Schoenfeld 殘差檢查比例風(fēng)險(xiǎn)，結(jié)果顯示sysbp，mitype這2個(gè)變量不滿足比例風(fēng)險(xiǎn)假定。

非線性檢測(cè)顯示bmi,los這兩個(gè)變量為非線性預(yù)測(cè)變量，經(jīng)過(guò)對(duì)比確定bmi的節(jié)點(diǎn)個(gè)數(shù)為3，los的節(jié)點(diǎn)個(gè)數(shù)為4。RSF樹(shù)的個(gè)數(shù)在1500之后模型趨于穩(wěn)定，因此RSF樹(shù)的個(gè)數(shù)設(shè)為1500，節(jié)點(diǎn)為預(yù)測(cè)變量個(gè)數(shù)的對(duì)數(shù)取整數(shù)3。表1展示神經(jīng)網(wǎng)絡(luò)參數(shù)調(diào)節(jié)范圍。通過(guò)1000 次對(duì)比確定Cox_DNN（DeepSurv）神經(jīng)網(wǎng)絡(luò)參數(shù)分別設(shè)置隱藏層神經(jīng)元個(gè)數(shù)為20，隱藏層層數(shù)為3，激活函數(shù)為tanh，學(xué)習(xí)率為0.001，其余參數(shù)選擇默認(rèn)結(jié)果。

從預(yù)測(cè)區(qū)分度上看，C-index值越高模型預(yù)測(cè)區(qū)分度越好，Cox_RCS的預(yù)測(cè)區(qū)分度整體高于其他3個(gè)模型，其次是Cox_DNN，但Cox_DNN有較多數(shù)值較小的離群值，表現(xiàn)出模型預(yù)測(cè)的不穩(wěn)定，由前面模擬可知在刪失率較高（＞40%）時(shí)，Cox_DNN模型預(yù)測(cè)效果不佳，而WHAS500的刪失率為57%，結(jié)果與模擬相符合；Cox和RSF的預(yù)測(cè)區(qū)分度相差不大（圖5）。從預(yù)測(cè)校準(zhǔn)度看，IBS 值越低模型預(yù)測(cè)校準(zhǔn)度越好，從圖中可看出Cox_RCS 的預(yù)測(cè)校準(zhǔn)度整體高于其他3 個(gè)模型，而Cox_DNN的預(yù)測(cè)校準(zhǔn)度最低。綜上所述，在含有非線性預(yù)測(cè)變量的數(shù)據(jù)中，使用Cox_RCS構(gòu)建預(yù)測(cè)模型，其整體的預(yù)測(cè)準(zhǔn)確度都高于其他3類模型。

圖5 WHAS數(shù)據(jù)集C-index(左)和IBS(右)的箱圖Fig.5 Concordance index(left)and integrated Brier scores(right)of WHAS dataset.

3 討論

盡管近年來(lái)陸續(xù)有學(xué)者提出深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)方法以解決生存分析模型中非線性擬合問(wèn)題，并且證明了在高維樣本數(shù)據(jù)中深度神經(jīng)網(wǎng)絡(luò)方法優(yōu)于傳統(tǒng)的線性Cox回歸，但對(duì)于擴(kuò)展的Cox模型如限制性立方樣條Cox模型與其他方法之間的比較研究還較少。并且目前還沒(méi)有研究綜合考慮樣本量結(jié)合刪失率時(shí)，各方法的優(yōu)劣情況。

從模擬和實(shí)例研究來(lái)看，4個(gè)模型的預(yù)測(cè)區(qū)分度指標(biāo)（C-index）和預(yù)測(cè)校準(zhǔn)度指標(biāo)（IBS），表現(xiàn)最好的是Cox_RCS。綜合模擬與實(shí)例的結(jié)果，Cox_RCS在滿足特定條件時(shí)的擬合效果優(yōu)于其它3個(gè)模型。因此，對(duì)于低維且存在非線性變量的數(shù)據(jù)，協(xié)變量若滿足比例風(fēng)險(xiǎn)假定，則推薦使用Cox_RCS。若不滿足比例風(fēng)險(xiǎn)假定，高刪失率（本文中≥40%）情況下可使用Cox_RCS，若為低刪失率（本文中＜40%）且大樣本量（本文中≥500），推薦使用Cox_DNN，否則使用Cox_RCS。

當(dāng)非線性變量存在時(shí)，傳統(tǒng)的擴(kuò)展Cox 模型如Cox_RCS便能較好的進(jìn)行預(yù)測(cè)，關(guān)于各個(gè)模型的運(yùn)行時(shí)間和難易程度，Cox_RCS的實(shí)現(xiàn)難點(diǎn)是需要篩選出具體的非線性變量，并確認(rèn)非線性變量是哪幾個(gè)，而機(jī)器學(xué)習(xí)和神經(jīng)網(wǎng)絡(luò)方法都無(wú)需識(shí)別出具體的非線性變量，而是直接把變量納入模型進(jìn)行擬合。Cox_RCS的優(yōu)點(diǎn)是可以直觀了解協(xié)變量的表現(xiàn)形式。隨機(jī)生存森林是機(jī)器學(xué)習(xí)樹(shù)模型，由樹(shù)節(jié)點(diǎn)不斷分裂構(gòu)成，機(jī)器學(xué)習(xí)方法構(gòu)建模型前需要調(diào)參，如隨機(jī)生存森林需要確定樹(shù)分裂的個(gè)數(shù)以及節(jié)點(diǎn)數(shù)以達(dá)到樹(shù)模型穩(wěn)定。深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)模型需要對(duì)更多的模型參數(shù)進(jìn)行調(diào)優(yōu)，通?？梢詫⒄{(diào)優(yōu)超參數(shù)分為3類：網(wǎng)絡(luò)參數(shù)、優(yōu)化參數(shù)、正則化參數(shù)。網(wǎng)絡(luò)參數(shù)包括神經(jīng)網(wǎng)絡(luò)的隱藏層層數(shù)，每一層的神經(jīng)元個(gè)數(shù)，激活函數(shù)等；優(yōu)化參數(shù)一般指學(xué)習(xí)率（learning rate）、批樣本數(shù)量（batch size）、不同優(yōu)化器的參數(shù)等；正則化參數(shù)的設(shè)置可以防止模型過(guò)擬合，一般包括權(quán)重衰減系數(shù)，丟棄法比率（dropout）等。因此深度學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)更為復(fù)雜，所花費(fèi)的時(shí)間也是最多的。

從結(jié)果解釋性上看，Cox和Cox_RCS的結(jié)果解讀起來(lái)較為直觀易懂，因?yàn)槠錁?gòu)建的是回歸模型，模型擬合的形式較機(jī)器學(xué)習(xí)方法簡(jiǎn)單，我們能根據(jù)每個(gè)變量的回歸系數(shù)直接解釋其對(duì)生存結(jié)局的影響。而Cox_DNN和RSF等機(jī)器學(xué)習(xí)方法使用前需要不斷探究最優(yōu)解的參數(shù)，通過(guò)建立模型來(lái)預(yù)測(cè)最終結(jié)果。但我們無(wú)法從機(jī)器學(xué)習(xí)的方法中直接解讀每個(gè)變量與生存結(jié)局的關(guān)系，其結(jié)果解釋性較差。

傳統(tǒng)擴(kuò)展Cox模型的優(yōu)點(diǎn)是模型較為簡(jiǎn)單，變量擬合的形式清楚明白，能輸出具體每個(gè)變量的回歸系數(shù)，并從中可以得到每個(gè)變量對(duì)生存時(shí)間的影響大??；缺點(diǎn)是無(wú)法在高維數(shù)據(jù)中擬合。隨機(jī)生存森林的優(yōu)點(diǎn)是可以高效識(shí)別并篩選變量，可以輸出變量的重要性排名，適用于高維數(shù)據(jù)；其缺點(diǎn)是無(wú)法直觀的得到變量與結(jié)局指標(biāo)的具體表現(xiàn)形式，并且數(shù)據(jù)的擬合模型在預(yù)測(cè)上并沒(méi)有優(yōu)于其他模型。深度生存神經(jīng)網(wǎng)絡(luò)方法優(yōu)點(diǎn)是可以很好處理大樣本數(shù)據(jù)，并且可以識(shí)別圖形進(jìn)行建模；其缺點(diǎn)是數(shù)據(jù)刪失程度對(duì)模型影響大，模型較為復(fù)雜，并且調(diào)參數(shù)需要花費(fèi)大量時(shí)間。因此實(shí)際應(yīng)用時(shí)需要結(jié)合數(shù)據(jù)，選擇適合的方法。