王純杰，任美慧，肖男男，張 鉞

(長(zhǎng)春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，吉林 長(zhǎng)春 130012)

0 引言

在生物醫(yī)學(xué)、社會(huì)科學(xué)、現(xiàn)代工業(yè)等領(lǐng)域，都需要對(duì)相關(guān)的生存資料進(jìn)行研究.研究的重點(diǎn)主要在生存時(shí)間以及相關(guān)的影響因素，生存分析是解決此類問(wèn)題主要方法.目前，生存分析主要涉及兩個(gè)主要研究方向：(1)描述生存時(shí)間分布，通過(guò)觀測(cè)研究對(duì)象的生存曲線，進(jìn)一步明確某一疾病的人群生存規(guī)律；(2)危險(xiǎn)函數(shù)估計(jì)以及相關(guān)影響因素分析，如在某一疾病的隨訪過(guò)程中，分析影響患者生存的影響因素，進(jìn)一步預(yù)測(cè)患者的生存時(shí)間分布[1].

目前，學(xué)者對(duì)于生存數(shù)據(jù)的分析與預(yù)測(cè)已經(jīng)取得了很大的進(jìn)步，因此，能夠通過(guò)繪制疾病隨訪過(guò)程中的疾病危險(xiǎn)率變化曲線，可以更加直觀地了解疾病的變化情況，并針對(duì)危險(xiǎn)率變化峰值采取相應(yīng)的治療措施.D.Faraggi和R.Simon[2]首次將神經(jīng)網(wǎng)絡(luò)模型應(yīng)用于Cox比例風(fēng)險(xiǎn)模型中，進(jìn)一步分析了影響生存時(shí)間的相關(guān)因素；K.Liestol[3]提出應(yīng)用神經(jīng)網(wǎng)絡(luò)來(lái)構(gòu)建離散時(shí)間模型和分段指數(shù)模型兩種估計(jì)危險(xiǎn)率變化情況的非參數(shù)靈活模型；此后，M.Fornili[4-5]再次利用神經(jīng)網(wǎng)絡(luò)分段指數(shù)模型處理右刪失數(shù)據(jù)，分別對(duì)乳腺癌以及腎細(xì)胞癌患者進(jìn)行實(shí)例分析，對(duì)病人的危險(xiǎn)變化情況進(jìn)行預(yù)測(cè)估計(jì)；L.Zhao[6]將深度學(xué)習(xí)方法應(yīng)用于生存數(shù)據(jù)的處理，同時(shí)估計(jì)出相應(yīng)的生存函數(shù)；進(jìn)而，C.Lee[7]將神網(wǎng)絡(luò)方法應(yīng)用到生存分析的競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型中，對(duì)危險(xiǎn)函數(shù)進(jìn)行了估計(jì)；另外，Y.X.Sun[8]采用分段指數(shù)模型來(lái)描述疾病潛伏過(guò)程，并將其用于處理帶有空間距離和時(shí)變協(xié)變量的疾病數(shù)據(jù).

本文將神經(jīng)網(wǎng)絡(luò)分段指數(shù)模型應(yīng)用于區(qū)間刪失數(shù)據(jù)，由于區(qū)間刪失數(shù)據(jù)缺乏準(zhǔn)確的觀測(cè)時(shí)間點(diǎn)，因此采用多重插補(bǔ)法下，對(duì)已有缺失數(shù)據(jù)進(jìn)行插補(bǔ)，如花琳琳等[9]在不同缺失值處理方法下，對(duì)隨機(jī)缺失數(shù)據(jù)處理效果的比較；紀(jì)忠光等[10]針對(duì)缺失數(shù)據(jù)進(jìn)行非參數(shù)插補(bǔ)；張琳琳等[11]基于數(shù)據(jù)增強(qiáng)算法，將多重插補(bǔ)方法應(yīng)用于失效時(shí)間為區(qū)間刪失型數(shù)據(jù)的加速失效時(shí)間模型的參數(shù)估計(jì)問(wèn)題中，并將有限區(qū)間刪失(但不是右刪失)數(shù)據(jù)插補(bǔ)成確切失效時(shí)間數(shù)據(jù).本文旨在應(yīng)用多重插補(bǔ)法結(jié)合神經(jīng)網(wǎng)絡(luò)分段指數(shù)模型對(duì)美國(guó)宇航員減壓病數(shù)據(jù)進(jìn)行分析，關(guān)注后期可能出現(xiàn)的多峰值變化情況.

1 研究模型與方法

1.1 分段指數(shù)模型

考慮分段指數(shù)模型并將隨訪時(shí)間劃分為J個(gè)不相鄰的區(qū)間，即

Aj=(aj-1，aj]，j=1，2，…，J，

其中a0=0，aJ=+∞.由此，可以假定風(fēng)險(xiǎn)函數(shù)在每個(gè)區(qū)間段內(nèi)均為常數(shù)，即

h(ti，xi)=λj(xi)，aj-1

在分段指數(shù)模型下，相應(yīng)的生存函數(shù)可以改寫(xiě)為

為了簡(jiǎn)單起見(jiàn)，當(dāng)滿足t≤a1時(shí)，上述生存函數(shù)的第一項(xiàng)可以看作是1.

1.2 數(shù)據(jù)插補(bǔ)法

針對(duì)于不完全數(shù)據(jù)的分析，通常采用插補(bǔ)法進(jìn)行分析(包括單次插補(bǔ)和多重插補(bǔ)兩種方法).單次插補(bǔ)是對(duì)缺失數(shù)據(jù)依據(jù)某種規(guī)則進(jìn)行一次填充，插補(bǔ)之后可以得到一個(gè)完整的數(shù)據(jù)集[12].多重插補(bǔ)則是對(duì)缺失值進(jìn)行多次補(bǔ)充，插補(bǔ)次數(shù)一般為5～20次，因此可以得到相應(yīng)個(gè)數(shù)的完整數(shù)據(jù)集.在此基礎(chǔ)上，針對(duì)每個(gè)數(shù)據(jù)集進(jìn)行分析，最終得出更為準(zhǔn)確的分析結(jié)果.采用imputeCens包[13]對(duì)區(qū)間刪失數(shù)據(jù)進(jìn)行多重插補(bǔ)，具體操作步驟如下：

(1)通過(guò)ic_par()函數(shù)擬合區(qū)間刪失數(shù)據(jù)的參數(shù)回歸模型，參數(shù)dist采用指數(shù)分布，滿足分段指數(shù)模型；

(2)輸入Ⅱ型觀測(cè)區(qū)間和左刪失數(shù)據(jù)代入模型，進(jìn)行數(shù)據(jù)處理；

(3)應(yīng)用imputeCens()函數(shù)，對(duì)區(qū)間刪失數(shù)據(jù)進(jìn)行插補(bǔ).

通過(guò)上述操作，可以將觀測(cè)區(qū)間以及左刪失數(shù)據(jù)轉(zhuǎn)換成具體的觀測(cè)時(shí)間點(diǎn)，與原數(shù)據(jù)構(gòu)成右刪失數(shù)據(jù)，對(duì)此進(jìn)行相應(yīng)的模型分析.

1.3 神經(jīng)網(wǎng)絡(luò)模型及估計(jì)方法

針對(duì)右刪失數(shù)據(jù)，為了能夠準(zhǔn)確地估計(jì)腫瘤患者的預(yù)后風(fēng)險(xiǎn)，采用了一種基于似然函數(shù)優(yōu)化的方法.假設(shè)有n名受試者(i=1，2，…，n)，其中第i位受試者的生存時(shí)間設(shè)為ti，且xi=(xi1，xi2，…，xip)表示個(gè)體i的p維協(xié)變量，則n個(gè)獨(dú)立個(gè)體的似然函數(shù)[14]為

對(duì)于上述似然函數(shù)，引入分段指數(shù)模型，則相應(yīng)的似然函數(shù)改寫(xiě)為

其中：Ji表示個(gè)體i所落入的最后一個(gè)區(qū)間；

下面對(duì)危險(xiǎn)函數(shù)h(t；x)進(jìn)行建模，利用神經(jīng)網(wǎng)絡(luò)構(gòu)建協(xié)變量與時(shí)間的函數(shù)關(guān)系，采用一種更靈活的模型解決了協(xié)變量之間的相依關(guān)系.采用以下前饋人工神經(jīng)網(wǎng)絡(luò)對(duì)危險(xiǎn)函數(shù)進(jìn)行建模，關(guān)系式為

其中

該網(wǎng)絡(luò)分為三層：輸入層是p維協(xié)變量和感興趣時(shí)間t；隱藏層包括H個(gè)神經(jīng)元；最后輸出層為J個(gè)神經(jīng)元(與分段指數(shù)模型的劃分區(qū)間數(shù)相同).另外，輸入層與隱藏層分別設(shè)定偏置單位“1”，相應(yīng)地選取Sigmoid函數(shù)作為激活函數(shù).

為優(yōu)化神經(jīng)網(wǎng)絡(luò)模型參數(shù)，選取極大似然函數(shù)的負(fù)對(duì)數(shù)作為網(wǎng)絡(luò)損失函數(shù)，通過(guò)最小化損失函數(shù)進(jìn)一步得出最優(yōu)化的網(wǎng)絡(luò)參數(shù)為

采用梯度下降法來(lái)最小化損失函數(shù)，并選定合適的學(xué)習(xí)率，通過(guò)計(jì)算，M.Akaike等[15]提出的評(píng)價(jià)指標(biāo)AIC值，來(lái)選擇最優(yōu)模型結(jié)構(gòu)；M.Stone[16]已經(jīng)證明AIC指標(biāo)與交叉驗(yàn)證的結(jié)果是等價(jià)的(具有計(jì)算優(yōu)勢(shì))，即

AIC=D+2β，

其中D為似然函數(shù)的負(fù)對(duì)數(shù)值，β為待估參數(shù)的個(gè)數(shù).

2 實(shí)例分析

本文主要研究高空飛行時(shí)減壓病的危險(xiǎn)率變化，數(shù)據(jù)NASAs Hypobaric Decompression Sickness(HDSD)來(lái)源于thregI包，該數(shù)據(jù)集記錄了美國(guó)空軍在低壓脫氮情況下肺動(dòng)脈內(nèi)氣栓(VGE IV)發(fā)病時(shí)間、年齡、性別、TR值(減壓前組織氮分壓與減壓終壓值的比值)、下肢運(yùn)動(dòng)等，以此來(lái)分析減壓病(DCS)在觀測(cè)時(shí)間內(nèi)的發(fā)病危險(xiǎn)率.數(shù)據(jù)集內(nèi)包含238個(gè)樣本，探測(cè)發(fā)生DCS患者自減壓至終壓的失效時(shí)間內(nèi)DCS危險(xiǎn)變化情況.

2.1 數(shù)據(jù)預(yù)處理

在HDSD中，每位受試者VGE發(fā)生時(shí)間在某一固定的檢測(cè)時(shí)間段內(nèi)，并未記錄準(zhǔn)確的觀測(cè)時(shí)間，針對(duì)區(qū)間Ⅱ型刪失數(shù)據(jù)，本文采取多重插補(bǔ)法，應(yīng)用imputeCens包同時(shí)依據(jù)分段指數(shù)分布規(guī)律預(yù)測(cè)出可能的生存時(shí)間分布，進(jìn)而將觀測(cè)區(qū)間以及左刪失數(shù)據(jù)轉(zhuǎn)換成具體的觀測(cè)時(shí)間點(diǎn)，與原數(shù)據(jù)構(gòu)成右刪失數(shù)據(jù)，將插補(bǔ)的準(zhǔn)確觀測(cè)時(shí)間以及相應(yīng)協(xié)變量應(yīng)用到神經(jīng)網(wǎng)絡(luò)分段指數(shù)模型中，進(jìn)一步分析減壓病的危險(xiǎn)變化趨勢(shì).

2.2 分析結(jié)果

減壓病是人體暴露在一定程度的低氣壓下而出現(xiàn)的一種特殊病癥，由于處于低氣壓下，使得體內(nèi)的氮?dú)庠跍p壓時(shí)出現(xiàn)過(guò)飽和現(xiàn)象，進(jìn)而生成氣泡而導(dǎo)致的病癥.根據(jù)受試者DCS發(fā)病的時(shí)間依從性信息(相關(guān)協(xié)變量)以及終檢時(shí)間(暴露期間)是否發(fā)病等信息，研究受試者在低壓暴露時(shí)間內(nèi)的危險(xiǎn)率變化情況.

根據(jù)所研究的數(shù)據(jù)集，為了選擇更加合適的網(wǎng)絡(luò)結(jié)構(gòu)，以M.Stone所提出的AIC準(zhǔn)則作為評(píng)價(jià)指標(biāo)，分別考慮了6，8，12個(gè)隱層單元數(shù)和學(xué)習(xí)率0.1，0.01以及0.001，根據(jù)AIC準(zhǔn)則的值，選擇具有較好性能的模型，相應(yīng)的結(jié)果如表1所示，因此選用6個(gè)隱層單元和學(xué)習(xí)率0.1來(lái)構(gòu)建神經(jīng)網(wǎng)絡(luò).

表1 AIC值分布Tab.1 AIC value distribution

在分析中，將隨訪時(shí)間0～6 h劃分為12個(gè)區(qū)間段，主要考慮了以下協(xié)變量：年齡、性別、TR值(即氮?dú)膺^(guò)飽和系數(shù))、是否運(yùn)動(dòng)以及VGE IV級(jí)對(duì)減壓病DCS危險(xiǎn)率變化的影響，具有不同特征的受試者在低壓情況下，一段隨訪時(shí)間內(nèi)，各類人群的危險(xiǎn)變化情況.

據(jù)美國(guó)潛水與高氣壓醫(yī)學(xué)會(huì)組織研究的相關(guān)資料顯示，發(fā)生減壓病的時(shí)間一般在剛進(jìn)入低氣壓環(huán)境時(shí)，尤其是進(jìn)入低壓環(huán)境1 h以內(nèi)，危險(xiǎn)率呈現(xiàn)遞增趨勢(shì)，隨著時(shí)間的增長(zhǎng)危險(xiǎn)率會(huì)出現(xiàn)下降趨勢(shì)，而后又會(huì)出現(xiàn)上升趨勢(shì)，直至出現(xiàn)平穩(wěn)狀態(tài).圖1是通過(guò)神經(jīng)網(wǎng)絡(luò)分段指數(shù)模型擬合的238位受試者的生存曲線，可以發(fā)現(xiàn)本文所提出的模型擬合的曲線與K-M估計(jì)得出的生存曲線基本重合，由此可以說(shuō)明該模型具有較好的評(píng)估性能.

圖1 生存曲線擬合圖Fig.1 Fitting diagram of survival curve

為進(jìn)一步分析危險(xiǎn)率變化情況，分別將受試者按照不同協(xié)變量劃分為不同的組別，如：年齡(劃分為三個(gè)區(qū)間段：20～30歲、30～40歲以及40～50歲)、性別、TR值(將比值劃分為1.6以下、1.6～1.7以及1.7以上)、是否運(yùn)動(dòng)等連續(xù)變量和分類變量對(duì)危險(xiǎn)率的影響.由圖2結(jié)果發(fā)現(xiàn)，處于不同年齡段的受試者在低壓情況下危險(xiǎn)率變化情況基本一致，說(shuō)明年齡對(duì)低壓環(huán)境下減壓病的危險(xiǎn)率并不造成影響；另外，相較于男性而言，女性在低壓環(huán)境中危險(xiǎn)率系數(shù)更低，更具有耐受性.

圖2 不同年齡段以及不同性別人群的危險(xiǎn)率變化曲線Fig.2 Risk curves of different age groups and different genders

另一方面，還分析了TR值以及在低壓環(huán)境中是否運(yùn)動(dòng)對(duì)危險(xiǎn)率變化的影響，TR值是組織內(nèi)氮分壓與減壓終壓值的比值，有關(guān)資料表明，TR比值一般在1.53～1.63范圍內(nèi)波動(dòng)，如果TR值過(guò)高，相應(yīng)的發(fā)病危險(xiǎn)率也會(huì)提高.由圖3結(jié)果顯示，當(dāng)TR值處于1.6水平以下，受試者的危險(xiǎn)率較低，但中間仍會(huì)出現(xiàn)危險(xiǎn)率的波動(dòng)變化；而當(dāng)危險(xiǎn)率處于1.6～1.7水平之間，受試者在最開(kāi)始時(shí)刻就會(huì)面臨著較大的危險(xiǎn)率，直至5 h以后趨于平緩；如果TR值超過(guò)1.7以上，危險(xiǎn)率會(huì)出現(xiàn)小高峰，相較于前者，具有更高的危險(xiǎn)率.

圖3 不同TR值含量下的危險(xiǎn)率變化曲線Fig.3 Variation curve of risk rate under different TR values

在低壓環(huán)境中，運(yùn)動(dòng)會(huì)增加減壓病發(fā)生的概率，因此在低壓情況下，活動(dòng)的人會(huì)比不活動(dòng)的人具有更高的危險(xiǎn)率，因?yàn)檫\(yùn)動(dòng)是與氮含量相關(guān)的，運(yùn)動(dòng)越多會(huì)導(dǎo)致組織內(nèi)的氮含量也越多，導(dǎo)致氮?dú)膺^(guò)飽和系數(shù)并不斷增加，進(jìn)而具有更高的危險(xiǎn)率，圖4結(jié)果也說(shuō)明了這一結(jié)論.

圖4 低壓情況下是、否為運(yùn)動(dòng)群體的危險(xiǎn)率變化曲線Fig.4 Curves of risk changes in the active and inactive groups under low pressure

3 結(jié)語(yǔ)

本文通過(guò)多重插補(bǔ)法將區(qū)間Ⅱ型刪失數(shù)據(jù)轉(zhuǎn)換為右刪失數(shù)據(jù)，并依托于神經(jīng)網(wǎng)絡(luò)分段指數(shù)模型對(duì)減壓病的危險(xiǎn)率變化進(jìn)行估計(jì)，與傳統(tǒng)的比例風(fēng)險(xiǎn)非線性模型相比，神經(jīng)網(wǎng)絡(luò)模型更具靈活性，解決了協(xié)變量之間的交互作用，進(jìn)而可以直接得到隨訪時(shí)間內(nèi)的個(gè)體危險(xiǎn)率的變化情況.