杜 婧，滕 婧，馬 卞，張奧鑫

華北電力大學(xué) 控制與計算機工程學(xué)院，北京102206

1 引言

乳腺癌在我國女性最常見的惡性腫瘤中排名第二，也是排名第六的惡性腫瘤致死原因[1]。中國女性的乳腺癌發(fā)病率和死亡率在近年來呈現(xiàn)出迅速增長的趨勢，且部分地區(qū)上升幅度顯著，增長較快[2]。我國女性的乳腺癌死亡率持續(xù)上升的現(xiàn)象，與歐美國家的顯著下降趨勢呈現(xiàn)出明顯的對比[3]。為降低患者的死亡率，需要分析挑選出更為有效的治療方案，或是針對患者定制個性化治療方案，這都依賴于有效的預(yù)后分析以及對患者生存情況的預(yù)測。

在癌癥的預(yù)后影響因素的研究中，研究者們發(fā)現(xiàn)有多種因素影響患者的預(yù)后水平，可將它們大致分為以下三類：一是患者的人口統(tǒng)計學(xué)和遺傳學(xué)特征，如患者的發(fā)病年齡、是否攜帶乳腺癌易感基因等；二是疾病特征，如腫瘤的位置、大小、組織學(xué)分級情況等；三是治療方案，如化療、免疫治療等[4]。目前已有大量統(tǒng)計學(xué)方法針對上述因素進(jìn)行分析研究，以求量化判斷這些因素對于患者預(yù)后情況的影響[5-6]。

預(yù)后分析的首要步驟是就是要判斷哪些因素的影響最為顯著。在上述列舉的多種因素中，有些因素之間高度相關(guān)，甚至冗余，無法提供更多的信息。由于對患者預(yù)后情況的跟蹤調(diào)查需要耗費大量的時間、經(jīng)濟成本，因此建立生存預(yù)測模型的首要步驟是選取顯著的預(yù)測特征，使預(yù)測模型盡可能的簡潔，即在得到幾乎相同信息量的前提下，選擇包含特征量最少的模型。目前，常用的特征選擇方法包括前向特征選擇、反向特征選擇，亦或是使用Cox回歸模型做單變量分析，選取影響權(quán)重較大的特征。本文在參考了大量文獻(xiàn)，并比較了各種特征選擇方法之間的優(yōu)缺點后[7-8]，選擇了常用的反向特征選擇方法，選取出最為顯著的因素進(jìn)行預(yù)后分析，其中LNR是最為顯著的疾病特征因素之一。

大量針對癌癥預(yù)后影響因素的研究表明，LNR是癌癥患者預(yù)后分析中最重要的特征之一，與患者預(yù)后的復(fù)發(fā)風(fēng)險尤其相關(guān)，且該項特征對于多種癌癥患者的預(yù)后生存率都有很好的預(yù)測效果。Cheng等人證明，轉(zhuǎn)移淋巴結(jié)比率是宮頸鱗狀細(xì)胞癌患者生存率的獨立預(yù)測因子[9]。Smith等人利用LNR與其他影響因素，使用基于貝葉斯方法的經(jīng)典Cox模型，預(yù)測胰腺癌患者的預(yù)后生存率，得到了相較于其他同類研究更好的效果[4]。在針對乳腺癌患者的研究中，Wu等人使用經(jīng)典Cox回歸模型分析了中山大學(xué)癌癥中心1998年至2007年的2 591份病歷記錄后發(fā)現(xiàn)，LNR水平較低的乳腺癌患者比LNR水平較高的患者有更高的總體生存率、無病生存率以及無轉(zhuǎn)移生存率，并證明LNR是上述生存率的獨立預(yù)測因子[10]。

在臨床上，LNR使用切片上檢測到的陽性淋巴結(jié)個數(shù)除以觀察到的淋巴結(jié)總數(shù)。但檢測所得的LNR值可能與患者實際的總體LNR值有較大偏差。其中最主要的原因是檢測樣本中的淋巴結(jié)總數(shù)僅為局部觀測值。這導(dǎo)致實驗所得的LNR值與患者的實際總體LNR值出入較大。相對而言，在切片檢測過程中觀察到的淋巴結(jié)數(shù)量越多，得到的LNR值也相應(yīng)地越準(zhǔn)確。因此，根據(jù)檢測到的淋巴結(jié)總數(shù)與陽性淋巴結(jié)數(shù)量計算得到的LNR檢測值應(yīng)當(dāng)看作是患者真實LNR值的一個估計[4]。為優(yōu)化LNR的估計精度，本研究引入其他相關(guān)的病理特征，基于邏輯回歸（Logistic Regression）算法對LNR值進(jìn)行估計，以期得到更接近患者真實總體LNR值的估計。本文仿真部分分別比較了基于邏輯回歸的LNR總體估計值與LNR局部切檢值對預(yù)后分析的影響。

如前所述，基于邏輯回歸模型得到的LNR總體估計值作為重要的臨床特征用于預(yù)后分析。目前，預(yù)后分析中廣泛采用經(jīng)典Cox回歸模型預(yù)測患者的生存率。該模型由英國統(tǒng)計學(xué)家Cox于1972年提出[11]，其基本思想是將患者的生存率表示為風(fēng)險函數(shù)，即個體在生存過程中某一單位時間內(nèi)的死亡概率。Cox回歸模型屬于半?yún)?shù)的生存分析模型[12]，與參數(shù)模型相比[13]，其條件更為寬松，不需要生存數(shù)據(jù)預(yù)先滿足某種特定的分布；與參數(shù)模型相比，其檢驗效能相對更高，可在未知數(shù)據(jù)的生存分布和基準(zhǔn)風(fēng)險函數(shù)的情況下，同時得到生存函數(shù)和基準(zhǔn)風(fēng)險函數(shù)。正是這些優(yōu)勢，使得經(jīng)典Cox回歸模型在提出后的幾十年時間里頗受青睞，得到了廣泛的應(yīng)用。然而，在經(jīng)典Cox回歸模型中，協(xié)變量系數(shù)始終為常數(shù)，無法反映隨時間的推移，預(yù)測變量對生存率所產(chǎn)生的動態(tài)影響。本研究基于貝葉斯方法構(gòu)建動態(tài)Cox回歸模型[14]，通過結(jié)合各參數(shù)的先驗知識與觀測數(shù)據(jù)共同推斷參數(shù)的后驗分布并持續(xù)更新，從而更好地捕捉不同時間區(qū)間中，預(yù)測變量對生存率的影響。

本研究總體流程框圖如圖1所示。首先從SEER（The Surveillance，Epidemiology，and End Results）中選取合適的患者樣本。選取條件為：女性乳腺癌患者；確診年齡在20歲到80歲之間；確診時間在2010年至2012年間；檢測到的淋巴結(jié)個數(shù)不小于1。另外考慮到不同的乳腺癌子類型患者之間的總體生存率有所差異[15]，在此僅選取“Her2-/HR+”子類型的乳腺癌患者。篩選后共得到4 402個樣本。樣本的簡要描述如表1所示。之后進(jìn)行反向特征選擇，挑選出用于估計LNR的特征，其余特征用于生存分析。根據(jù)赤池信息準(zhǔn)則（Akaike Information Criterion，AIC）指標(biāo)選取得到的與LNR相關(guān)性最大的特征分別為：淋巴結(jié)總數(shù)、陽性淋巴結(jié)個數(shù)、M分期和N分期，然后將它們作為輸入訓(xùn)練邏輯回歸模型并估計LNR值。在預(yù)后分析部分，使用LNR的總體估計值，以及樣本中的其他特征，包括患者的年齡、腫瘤大小和T分期，基于貝葉斯方法構(gòu)建動態(tài)Cox回歸模型以預(yù)測患者生存率。

圖1 整體流程圖

2 數(shù)據(jù)與方法

2.1 基于邏輯回歸模型的LNR估計

臨床應(yīng)用的LNR值是用切片上檢測到的陽性淋巴結(jié)個數(shù)除以淋巴結(jié)總數(shù)計算所得。但實際情況下，這樣得到的LNR局部切檢值難以準(zhǔn)確地反映患者的總體LNR值。因此本研究綜合考慮患者其他相關(guān)的病理信息與LNR的關(guān)聯(lián)，使用邏輯回歸模型估計患者的總體LNR值。首先，將LNR局部切檢值作為響應(yīng)特征，通過反向特征選擇篩選出相關(guān)病理特征作為輸入特征，再通過擬合邏輯回歸模型，計算協(xié)變量系數(shù)，最后根據(jù)相應(yīng)的協(xié)變量系數(shù)估計患者總體LNR值。

表1 數(shù)據(jù)集樣本特征

邏輯回歸模型是一種常用的機器學(xué)習(xí)模型，其形式較為簡單直觀，且具有良好的可解釋性。本文中利用邏輯回歸模型的值域在[0，1]范圍內(nèi)進(jìn)行回歸分析，估計具有相同值域的LNR總體值。本研究中，邏輯回歸模型的基本形式如下：對于樣本i，其響應(yīng)值，即LNR值為yi；經(jīng)過反向特征選擇篩選出4個與LNR相關(guān)的病理特征作為作為輸入特征，分別為陽性淋巴結(jié)個數(shù)X1、淋巴結(jié)總數(shù)X2、M分期X3和N分期X4。根據(jù)邏輯回歸模型，樣本i的LNR值yi與其對應(yīng)的預(yù)測特征Xi的關(guān)系為：其中，β0為常數(shù)項，β1,β2,β3,β4是每個預(yù)測特征對應(yīng)的協(xié)變量系數(shù)，β是上述協(xié)變量系數(shù)組成的向量。Xi,1,Xi,2,Xi,3,Xi,4是樣本i的4個預(yù)測特征值，Xi是上述預(yù)測特征值組成的向量。采用訓(xùn)練集數(shù)據(jù)擬合邏輯回歸模型后，得到所有預(yù)測特征對應(yīng)的協(xié)變量系數(shù)β。之后，根據(jù)系數(shù)β和預(yù)測特征Xi，再代入式（1）可得LNR的總體估計值。

2.2 基于貝葉斯方法的動態(tài)Cox回歸模型

經(jīng)典Cox回歸模型的基本形式為：

在經(jīng)典Cox回歸模型中，不同時間節(jié)點的協(xié)變量系數(shù)βs保持不變。但實際應(yīng)用中，各預(yù)測變量對患者生存率的影響往往是隨時間變化的。為此，基于貝葉斯方法的動態(tài)回歸模型將不同時間點的協(xié)變量系數(shù)表示為βs(t)，并基于貝葉斯方法，結(jié)合生存數(shù)據(jù)求取其后驗分布。該方法由Wang等人提出，本文僅作簡要表述。

基于貝葉斯方法的動態(tài)Cox回歸模型的基本形式為：

其中，Z是所有樣本的預(yù)測變量組成的矩陣，λ0(t)表示t時刻的基準(zhǔn)風(fēng)險，βs(t)表示t時刻的協(xié)變量系數(shù)向量。在該動態(tài)模型中，λ0(t)和βs(t)均假設(shè)為左連續(xù)的階躍函數(shù)。模型中，需要估計基準(zhǔn)風(fēng)險函數(shù)λ0(t)和協(xié)變量系數(shù)向量βs(t)具體的階躍次數(shù)和相應(yīng)的階躍值與階躍時間點，以達(dá)到動態(tài)估計參數(shù)的目的。令：Θ=，該集合包含了所有的未知參數(shù)，使用θ(t)來指代lnλ0(t)或βs(t)中的一個分量。

所有的未知參數(shù)都根據(jù)數(shù)據(jù)樣本估計得到。對于樣本i(i=1,2,…,n)，令Ti表示事件“患者死亡”發(fā)生的時間。若Ti已知，則該樣本數(shù)據(jù)為完整生存數(shù)據(jù)。若僅能確定Ti∈[Li,Ri)，且Ri為有限值，則該樣本數(shù)據(jù)為區(qū)間刪失數(shù)據(jù)；若Ri=∞，則該樣本數(shù)據(jù)為右刪失數(shù)據(jù)。同時，根據(jù)所有樣本的Ti構(gòu)造由等距點組成的時間網(wǎng)格，G={0=s0＜s1＜…＜sK＜sK+1=∞}。令：Δk=sksk-1表示第k個網(wǎng)格區(qū)間的寬度，并計：λk=λ0(sk)，βk=β(sk)。最后，令：Dobs={ }Ti∈[Li,Ri),Zi;i=1,2,…,n，該集合表示所有樣本的生存信息，以及與生存分析有關(guān)的預(yù)測變量的信息。

在動態(tài)模型中，采用貝葉斯方法用于估計模型的參數(shù)。貝葉斯方法的基本形式為：

該式表示，參數(shù)的后驗分布正比于參數(shù)的聯(lián)合先驗分布p(Θ)與樣本似然L(Z|Θ)的乘積。其中，樣本似然L(Z|Θ)可表示為：

式中，n為樣本總數(shù)。其中任一樣本i的似然貢獻(xiàn)為：

其中：

上式中，I(·)為指示函數(shù)，若·為真，則I(·)=1，否則I(·)=0。

前文提到，在基于貝葉斯方法的動態(tài)Cox回歸模型中，參數(shù)Θ是根據(jù)生存數(shù)據(jù)動態(tài)階躍變化的。且變化的次數(shù)與變化程度需經(jīng)貝葉斯方法推理計算。假設(shè)參數(shù)的階躍次數(shù)為J，且J服從從1到K的離散的均勻分布。對于固定的J，階躍的時間點0＜τ1＜τ2＜…＜τJ=sK是從時間網(wǎng)格集合G中隨機選取的。對于參數(shù)集合Θ中的所有分量，其先驗假設(shè)為：

其中，ω服從參數(shù)為α和η的逆伽馬分布，是模型中的超參數(shù)。參數(shù)c、α和η用于估計不同時間區(qū)間的參數(shù)，且都是預(yù)先定義的常數(shù)。每個時間區(qū)間的協(xié)變量系數(shù)的先驗分布都與前一個時間區(qū)間該參數(shù)的先驗分布假設(shè)相關(guān)。如圖2所示是參數(shù)的先驗分布假設(shè)間的關(guān)系。圓框中的參數(shù)是預(yù)先定義的常數(shù)，方框中的參數(shù)僅知其服從的特定分布，具體數(shù)值需要估計。

圖2 參數(shù)間的先驗關(guān)系

根據(jù)前文所述的貝葉斯框架和參數(shù)的動態(tài)先驗關(guān)系，θ(t)和ω的聯(lián)合先驗分布，π(θ(t),ω)正比于下式：

3 仿真結(jié)果

本研究運用R Studio 1.0.143進(jìn)行數(shù)據(jù)處理和分析，使用的R語言版本為3.4.4。在SEER中篩選得到了4 402個樣本。在實際仿真中，將這些樣本隨機劃分為訓(xùn)練集與測試集。其中訓(xùn)練集包含3 000個樣本，測試集包含剩余的1 402個樣本。

3.1 LNR估計

經(jīng)過反向特征選擇后，共有4個特征用于估計LNR值，分別是陽性淋巴結(jié)個數(shù)、淋巴結(jié)總數(shù)、M分期和N分期。使用這4個特征時，邏輯回歸模型有最低的AIC值1968。AIC值低說明該模型用較少的參數(shù)獲得了足夠的擬合度。邏輯回歸模訓(xùn)練后的部分預(yù)測特征的協(xié)變量系數(shù)列于表2中。從表中可以看出，淋巴結(jié)總數(shù)和陽性淋巴結(jié)個數(shù)這兩個特征的標(biāo)準(zhǔn)差最低，說明這兩個特征與LNR的關(guān)系最為顯著。

表2 部分預(yù)測變量的系數(shù)

為了判斷是否出現(xiàn)過擬合問題，分別計算了訓(xùn)練集和測試集數(shù)據(jù)擬合邏輯幾率回歸模型后的MSE。經(jīng)過計算，訓(xùn)練集數(shù)據(jù)擬合模型后的MSE值為0.018，測試集數(shù)據(jù)的MSE值為0.020。這兩者在同一數(shù)量級上且差距較小。據(jù)此判斷，訓(xùn)練后的邏輯幾率回歸模型沒有出現(xiàn)過擬合現(xiàn)象。在后續(xù)的計算過程中，使用該模型估計的LNR值是可行的。

3.2 生存分析

為了檢驗LNR的總體估計值對患者生存情況的預(yù)測效果，在生存分析部分，本研究分別使用了兩個數(shù)據(jù)集。兩者均包含患者的生存信息、T和N分期信息、確診年齡和腫瘤大小；唯一不同之處在于，數(shù)據(jù)集1使用LNR局部切檢值，數(shù)據(jù)集2使用LNR總體估計值。再者，為驗證基于貝葉斯方法的動態(tài)Cox回歸模型有能夠更好地反映預(yù)測變量在不同時間階段對患者生存率的影響，本文分別使用了經(jīng)典Cox回歸模型和動態(tài)Cox回歸模型。因此，不同的數(shù)據(jù)集與不同的生存分析模型兩兩組合，共產(chǎn)生了4個模型模型以比較它們之間的結(jié)果差異。4個模型的主要特點如表3所示。對所有模型均設(shè)置500次吉布斯采樣。

表3 模型特點

本文中使用對數(shù)偽邊際似然（Log Pseudo Marginal Likelihood，LPML）作為生存分析模型的評價指標(biāo)。對模型來說，該指標(biāo)數(shù)值越大，說明樣本越支持該模型。4個模型的LPML值分別是：-710.25，-708.79，-701.14以及-698.49。正如表3所示，相比于Model_1，Model_2的LPML值較大。同樣，Model_4的LPML值大于Model_3的LPML值。根據(jù)這一數(shù)值比較，可以看出，使用LNR總體估計值有相對較好的模型擬合效果。另外，Model_3和Model_4的LPML值均高于Model_1和Model_2的LPML值。這一結(jié)果說明基于貝葉斯方法的動態(tài)Cox回歸模型相較于經(jīng)典Cox回歸模型，能夠更好地擬合生存數(shù)據(jù)。

為了更直觀判斷模型的效果，繪制了4個模型對訓(xùn)練集和測試集數(shù)據(jù)分析后得到的總體生存率-時間曲線。與之相比較的是利用Kaplan-Meier方法（后文簡稱KM方法）繪制的患者總體生存率-時間曲線。與KM曲線更加接近，則說明該模型的預(yù)測效果相對更好。圖3中所示的是4個模型所預(yù)測的訓(xùn)練集和測試集數(shù)據(jù)的生存率-時間曲線。從圖中可看出，使用經(jīng)典Cox回歸模型的Model_1和Model_2在訓(xùn)練集和測試集上的表現(xiàn)均遜色于使用動態(tài)Cox模型的Model_3和Model_4，使用動態(tài)Cox模型預(yù)測得到的總體生存率-時間曲線更接近實際情況。另外，Model_3與Model_4在訓(xùn)練集和測試集上得到的曲線較為接近，但仍然可以看出，Model_4曲線有相對較好的預(yù)測效果，并且Model_4有更低的LPML值。這說明，使用LNR的總體估計值能夠更為準(zhǔn)確地預(yù)測生存率-時間特性曲線。

4 討論與總結(jié)

本文指出了乳腺癌預(yù)后分析中存在的兩個重要問題提出了相應(yīng)的解決方案：一是實驗檢測所得的LNR值受觀測誤差影響較大，對后續(xù)的生存分析過程產(chǎn)生偏差的；二是經(jīng)典的Cox回歸模型無法捕捉到不同時間區(qū)間內(nèi)，癌癥相關(guān)因素對于患者生存率的動態(tài)影響。針對第一個問題，本研究首先使用邏輯回歸估計LNR總體值，之后使用LNR總體估計值與其他預(yù)測變量信息以及生存數(shù)據(jù)擬合基于貝葉斯方法的動態(tài)Cox回歸模型。與使用LNR局部切檢值相比，使用估計值能夠降低較小的淋巴結(jié)檢測總數(shù)對于LNR值的影響，以及患者間個體差異對于LNR的影響。針對第二個問題，使用基于貝葉斯方法的動態(tài)Cox回歸模型能夠更好地捕捉不同時間階段，預(yù)測特征對于患者生存率的影響；同時，在估計系數(shù)的過程中，使用了能夠結(jié)合參數(shù)先驗分布的貝葉斯方法，這使得對于參數(shù)的預(yù)測更為準(zhǔn)確。

圖3 預(yù)測的訓(xùn)練集和測試集生存率-時間曲線

本文中所使用的數(shù)據(jù)集是SEER中的部分女性乳腺癌患者數(shù)據(jù)，算法實現(xiàn)使用R語言。為了驗證文中所述方法的性能，使用LPML值作為衡量模型性能的指標(biāo)。仿真結(jié)果表明，使用LNR估計值和基于貝葉斯方法的動態(tài)Cox回歸方法的模型有最高的LPML值，這說明數(shù)據(jù)最為支持該模型。另外，為了驗證模型的預(yù)測效果，使用KM方法計算測試集數(shù)據(jù)的生存率-時間曲線，并以此作為基準(zhǔn)，與使用邏輯回歸模型估計LNR和基于貝葉斯方法的動態(tài)Cox生存分析模型所預(yù)測的生存率-時間曲線進(jìn)行比較。結(jié)果顯示，兩個曲線有較多重合之處，且趨勢一致。在未來的研究中，可以繼續(xù)探索LNR這一指標(biāo)對于癌癥患者生存率的預(yù)測效果，并嘗試使用其他的機器學(xué)習(xí)方法，如決策樹、隨機森林等，來估計LNR值。