基于貝葉斯方法的乳腺癌預后分析

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

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

1 引言

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

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

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

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

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

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

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

圖1 整體流程圖

2 數據與方法

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

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

表1 數據集樣本特征

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

2.2 基于貝葉斯方法的動態Cox回歸模型

經典Cox回歸模型的基本形式為：

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

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

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

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

在動態模型中，采用貝葉斯方法用于估計模型的參數。貝葉斯方法的基本形式為：

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

式中，n為樣本總數。其中任一樣本i的似然貢獻為：

其中：

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

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

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

圖2 參數間的先驗關系

根據前文所述的貝葉斯框架和參數的動態先驗關系，θ(t)和ω的聯合先驗分布，π(θ(t),ω)正比于下式：

3 仿真結果

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

3.1 LNR估計

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

表2 部分預測變量的系數

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

3.2 生存分析

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

表3 模型特點

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

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

4 討論與總結

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

圖3 預測的訓練集和測試集生存率-時間曲線

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