(山東科技大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，山東 青島 266590)

肺癌是發(fā)病率和死亡率增長最快、對人類健康和生命威脅最大的惡性腫瘤之一，世界衛(wèi)生組織國際癌癥研究機(jī)構(gòu)發(fā)布最新報(bào)告稱肺癌死亡人數(shù)最多，占預(yù)計(jì)癌癥死亡總?cè)藬?shù)的18.4%。另外，肺癌還具有預(yù)后差的特點(diǎn)，影響其預(yù)后的因素主要包括患者個體相關(guān)因素、腫瘤相關(guān)因素和治療相關(guān)因素[1]。目前，臨床醫(yī)學(xué)主要根據(jù)手術(shù)病理分期判斷預(yù)后，由于考慮影響腫瘤發(fā)生的因素減少，其預(yù)測效果較差[2]，因此，建立一個適用于臨床醫(yī)學(xué)且考慮多因素的肺癌預(yù)后模型具有重要意義。

早期，國內(nèi)外在疾病預(yù)測方面多采用統(tǒng)計(jì)學(xué)中的COX回歸方法構(gòu)建模型。隨著數(shù)據(jù)挖掘技術(shù)被應(yīng)用到醫(yī)學(xué)研究領(lǐng)域，眾多學(xué)者采用機(jī)器學(xué)習(xí)的方法進(jìn)行疾病研究。劉雅琴等[3]使用logistic回歸、決策樹和人工神經(jīng)網(wǎng)絡(luò)方法研究預(yù)后模型的預(yù)測效果，是國內(nèi)機(jī)器學(xué)習(xí)領(lǐng)域研究腫瘤預(yù)測的有效嘗試。Kim等[4]利用支持向量機(jī)預(yù)測了乳腺癌患者術(shù)后5年生存情況。Chen等[5]對4個醫(yī)療機(jī)構(gòu)的非小細(xì)胞肺癌患者，使用人工神經(jīng)網(wǎng)絡(luò)建立了患者生存狀況風(fēng)險(xiǎn)模型。牟冬梅等[6]構(gòu)建了妊娠高血壓綜合征危險(xiǎn)因素決策樹預(yù)測模型。宋一鳴[7]基于SEER數(shù)據(jù)庫使用決策樹、神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)、Logistic回歸、深度神經(jīng)網(wǎng)絡(luò)等分別建立了肺癌患者預(yù)后的相關(guān)研究模型。

復(fù)發(fā)、轉(zhuǎn)移、風(fēng)險(xiǎn)評估及生存情況評價(jià)是腫瘤預(yù)后的主要研究內(nèi)容[8]，本研究針對患者術(shù)后5年后生存情況進(jìn)行研究。選取SEER數(shù)據(jù)庫[9]中部分肺癌患者的數(shù)據(jù)，根據(jù)相關(guān)研究提取影響患者生存情況的預(yù)后因素，通過貝葉斯網(wǎng)絡(luò)方法利用訓(xùn)練集構(gòu)建肺癌預(yù)后模型，其中采用最大最小爬山算法建立模型，并采用貝葉斯估計(jì)進(jìn)行模型參數(shù)學(xué)習(xí)與概率推理，最后將本研究模型與Logistic回歸、人工神經(jīng)網(wǎng)絡(luò)、決策樹及支持向量機(jī)方法在測試集上進(jìn)行分類實(shí)驗(yàn)比較，驗(yàn)證所建立模型的有效性。

1 數(shù)據(jù)來源及變量選擇

數(shù)據(jù)選自美國國立癌癥研究所“監(jiān)測、流行病學(xué)和結(jié)果數(shù)據(jù)庫”(SEER數(shù)據(jù)庫)[9]中2008年至2014年期間被確診為肺癌的患者，其中包括5年內(nèi)直接因癌細(xì)胞致死和隨訪期滿5年仍生存的患者。刪除數(shù)據(jù)缺失嚴(yán)重、記錄錯誤及因非肺癌致死的患者記錄，最終共計(jì)879位患者數(shù)據(jù)。

表1 肺癌患者變量信息

根據(jù)腫瘤信息，參考文獻(xiàn)[10,11]和其他相關(guān)研究[2-3,7]中提及的與患者生存相關(guān)的預(yù)后因素，從數(shù)據(jù)庫中導(dǎo)出包含這些因素的16個信息變量，具體如表1所示，其中后四項(xiàng)為連續(xù)型變量，其余為離散型變量。

5年后生存情況是預(yù)后效果的重要評價(jià)指標(biāo)，所以選擇患者術(shù)后5年生存情況為結(jié)果變量(生存時間以月為單位)。生存時間60個月及以上患者生存情況為“生存”(記為1)，低于60個月的患者生存情況為“死亡”(記為0)。

2 肺癌預(yù)后模型的建立

2.1 特征選擇

為提高模型的預(yù)測準(zhǔn)確性，對上述16個信息變量進(jìn)行特征選擇。首先，利用SPSS進(jìn)行卡方檢驗(yàn)，在p<0.05下通過檢驗(yàn)的變量有12個，分別為：婚姻狀況、組織學(xué)分級、腫瘤分期、轉(zhuǎn)移程度、擴(kuò)散程度、淋巴結(jié)累積程度、手術(shù)類型、是否放療、確診年齡、腫瘤大小、淋巴結(jié)受檢數(shù)量及淋巴結(jié)陽性數(shù)量。然后，在卡方檢驗(yàn)基礎(chǔ)上利用SPSS進(jìn)行Logistic回歸分析，在p<0.05下最終篩選出的特征變量有6個，分別為：組織學(xué)分級、腫瘤分期、確診年齡、腫瘤大小、淋巴結(jié)受檢數(shù)量及淋巴結(jié)陽性數(shù)量。篩選結(jié)果如表2所示。

2.2 數(shù)據(jù)離散化

(1)

表2 Logistic回歸分析篩選變量結(jié)果

2.3 模型建立方法

在疾病生存預(yù)測方面，傳統(tǒng)的統(tǒng)計(jì)模型難以計(jì)算后驗(yàn)概率，不能直觀地表示變量之間的關(guān)系，本研究利用貝葉斯網(wǎng)絡(luò)方法建立肺癌預(yù)后模型。

貝葉斯網(wǎng)絡(luò)是一個帶參數(shù)的有向無環(huán)圖，用二元組〈G,Θ〉表示，其中G=(V,E)表示節(jié)點(diǎn)關(guān)系的有向無環(huán)圖，稱為貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)，節(jié)點(diǎn)集合V={X1,X2,…,Xn}表示隨機(jī)變量，有向邊集合E={eij|Xi→Xj,i,j=1,2,…,n}表示變量之間的依賴關(guān)系；Θ={Θ1,Θ2,…,Θn}表示節(jié)點(diǎn)Xi的條件概率，稱為貝葉斯網(wǎng)絡(luò)參數(shù)，節(jié)點(diǎn)Xi的參數(shù)Θi表示其自身和父節(jié)點(diǎn)集Pa(Xi)的條件概率分布，即Θi=P(Xi|Pa(Xi))。另外，任意給定的貝葉斯網(wǎng)絡(luò)都滿足馬爾科夫條件，即?Xi∈V，Xi獨(dú)立于除其父節(jié)點(diǎn)集合Pa(Xi)外的所有非子孫節(jié)點(diǎn)，因此，變量集V=(X1,X2,…,Xn}聯(lián)合概率分布可分解為：

(2)

貝葉斯網(wǎng)絡(luò)模型用有向無環(huán)圖表示變量之間的依賴和獨(dú)立關(guān)系，用條件概率分布刻畫變量對其父節(jié)點(diǎn)的依賴關(guān)系，因此，建立貝葉斯網(wǎng)絡(luò)模型包括兩部分：①確定變量間關(guān)系，找到網(wǎng)絡(luò)結(jié)構(gòu)，即結(jié)構(gòu)學(xué)習(xí)；②確定每個節(jié)點(diǎn)的條件概率表，即參數(shù)學(xué)習(xí)。

2.3.1 結(jié)構(gòu)學(xué)習(xí)方法

利用最大最小爬山(Max-Min hill-climbing, MMHC)算法對貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行學(xué)習(xí)。該算法是Tsamardinos等[12]于2006年提出的一種經(jīng)典的貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)算法，結(jié)合了依賴分析和評分搜索等方法，分為兩個階段進(jìn)行學(xué)習(xí)：第一階段利用MMPC(max-min parents and children)算法確定出每個節(jié)點(diǎn)的候選父子節(jié)點(diǎn)集，構(gòu)建出貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)的無向框架；第二階段利用貪婪爬山算法對已經(jīng)得到的網(wǎng)絡(luò)結(jié)構(gòu)的框架進(jìn)行搜索評分，找出使評分函數(shù)最大的網(wǎng)絡(luò)結(jié)構(gòu)。

MMPC算法是從給定數(shù)據(jù)集中利用最大-最小啟發(fā)式策略確定目標(biāo)變量T的候選父子節(jié)點(diǎn)(candidate parents and children，CPC)集，分為兩個階段。第一階段通過定義一個關(guān)聯(lián)度函數(shù)來確定其他變量與目標(biāo)變量T在給定CPC下的條件依賴程度，函數(shù)值越大表示變量間的條件依賴關(guān)系越強(qiáng)；當(dāng)函數(shù)值為零時，表示變量間沒有依賴關(guān)系，也就是條件獨(dú)立。最大最小啟發(fā)式策略每次選擇與目標(biāo)變量T在給定CPC條件下最小關(guān)聯(lián)度最大的那個變量進(jìn)入CPC，當(dāng)除了CPC中變量所有其他變量與目標(biāo)變量T在給定CPC條件下都條件獨(dú)立時，第一階段停止。第二階段檢驗(yàn)候選父子節(jié)點(diǎn)集CPC中的變量，移去不該有變量，即對于CPC中的變量X，如果存在CPC的子集S使得Assoc(X,T|S)，則將變量X從CPC中移去。

變量X與T在給定變量集Z下的關(guān)聯(lián)度函數(shù)定義為：

(3)

(4)

其中，S表示變量集Z的子集。MMPC算法如下所示：

算法1：MMPC算法

輸入：目標(biāo)變量T，數(shù)據(jù)集D

輸出：目標(biāo)變量T的候選父子節(jié)點(diǎn)集CPC

第一階段：

1：令CPC≠φ；

2：WhileCPC不再變化 do

3： 〈F,assocF〉=MaxMinHeuristic(T,CPC)

4： ifassocF≠0 then

5：CPC=CPC∪F

6： end if

7：end

第二階段：

8：for 任意X∈CPC

9： if 存在S?CPC，使Assoc(X,T|S)=0，即Ind(X,T|S)then

10：CPC=CPC{X}

11： end if

12：end for

13：返回CPC

子程序MaxMinHeuristic(T,CPC)

輸入：目標(biāo)變量T，CPC子集

輸出：以CPC為條件集，與T的最小關(guān)聯(lián)度最大的變量

14：assocF=maxX∈VMinAssoc(X,T|CPC)

15：F=arg maxX∈VMinAssoc(X,T|CPC)

16：返回 〈F,assocF〉

MMHC算法第二階段利用貪婪爬山搜索在結(jié)構(gòu)空間中搜索評分最高的網(wǎng)絡(luò)結(jié)構(gòu)，評分函數(shù)采用BDeu評分。該階段的貪婪爬山搜索從空圖開始，每一步搜索的過程是：首先在不產(chǎn)生有向環(huán)的情況下，對當(dāng)前所得模型分別執(zhí)行一次加邊、減邊、轉(zhuǎn)邊操作得到一系列候選模型，并計(jì)算出每個候選模型的評分；然后將最大評分的候選模型與當(dāng)前模型比較，若最大評分的候選模型評分大，則將其作為下一個當(dāng)前模型繼續(xù)搜索，否則停止搜索并返回當(dāng)前模型[13]。

在MMHC算法中，貪婪爬山搜索將每個節(jié)點(diǎn)的搜索空間限制在其候選父子節(jié)點(diǎn)集上，即僅考慮當(dāng)Y∈CPCX時添加邊Y→X，此約束顯著降低了搜索空間的復(fù)雜性，提高了算法的效率。MMHC算法如下：

算法2：MMHC算法

輸入：數(shù)據(jù)集D

輸出：有向無環(huán)圖

1：for 所有變量X∈Vdo

2：CPCX=MMPC(X,D)

3：end for

4：從空圖出發(fā)執(zhí)行貪婪爬山搜索的3個搜索算子加邊、減邊和轉(zhuǎn)邊。

當(dāng)且僅當(dāng)Y∈CPCX時，添加有向邊Y→X。

5：返回最高得分的有向無環(huán)圖

2.3.2 參數(shù)學(xué)習(xí)方法

參數(shù)學(xué)習(xí)在統(tǒng)計(jì)學(xué)中主要有最大似然估計(jì)和貝葉斯估計(jì)兩種基本方法，本研究采用貝葉斯估計(jì)[13]對貝葉斯網(wǎng)絡(luò)參數(shù)進(jìn)行學(xué)習(xí)。

設(shè)一個貝葉斯網(wǎng)絡(luò)有n個節(jié)點(diǎn)V={X1,X2,…,Xn}，其中節(jié)點(diǎn)Xi有ri種取值，其父節(jié)點(diǎn)π(Xi)的取法有qi種組合。若Xi無父節(jié)點(diǎn)，則qi=1。該貝葉斯網(wǎng)絡(luò)的參數(shù)為：

θijk=P(Xi=k|π(Xi)=j)(i=1,2,…,n;j=1,2,…,qi;k=1,2,…,ri)。

(5)

用θ表示所有θijk組成的參數(shù)向量。設(shè)D={D1,D2,…,Dm}是一組關(guān)于貝葉斯網(wǎng)絡(luò)的獨(dú)立同分布的完整數(shù)據(jù)，則θ的似然函數(shù)為：

(6)

其中Nijk表示數(shù)據(jù)集D中滿足Xi=k和π(Xi)=j的樣本數(shù)量。假設(shè)參數(shù)θ的先驗(yàn)概率分布服從狄利克雷分布Dir(αij1,αij2,…,αijri)(i=1,2,…n;j=1,2,…,qi)，則：

(7)

(8)

從而，p(θ|D)～Dir(Nij1+αij1,Nij2+αij2,…,Nijri+αijri)(i=1,2,…,n;j=1,2,…,qi)，因此，參數(shù)θ的貝葉斯估計(jì)為[13]：

(9)

2.4 模型建立與結(jié)果分析

將最終保留的879條完整觀測記錄的數(shù)據(jù)集按照7∶3的比例分為訓(xùn)練集和測試集，其中訓(xùn)練集樣本為615個，測試集樣本為264個。訓(xùn)練集用來構(gòu)建預(yù)后模型，測試集用來預(yù)測性能，對預(yù)后模型進(jìn)行評價(jià)。

實(shí)驗(yàn)環(huán)境基本配置為CPU 2.53 GHz、RAM 2.00 GB，操作系統(tǒng)為Windows 7，在MATLAB 7.0上利用貝葉斯網(wǎng)絡(luò)工具箱Full BNT1.0.4。對貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)的學(xué)習(xí)，利用MATLAB編程，最終得到的肺癌預(yù)后模型如圖1所示，其中，7個節(jié)點(diǎn)為表2所示的6個特征變量及1個結(jié)果變量，節(jié)點(diǎn)之間的連線表明變量間的相互影響關(guān)系。實(shí)驗(yàn)結(jié)果顯示，腫瘤大小和組織學(xué)分級通過影響腫瘤分期間接地影響患者的生存情況；而確診時的年齡、腫瘤分期、淋巴結(jié)受檢數(shù)量以及淋巴結(jié)陽性數(shù)量直接影響患者的生存情況，這一結(jié)論符合醫(yī)學(xué)實(shí)際。

1-確診時年齡；2-腫瘤大??；3-組織學(xué)分級；4-腫瘤分期；5-淋巴結(jié)受檢數(shù)量；6-淋巴結(jié)陽性數(shù)量；7-生存情況。

進(jìn)一步，對圖1得到的預(yù)后模型進(jìn)行貝葉斯網(wǎng)絡(luò)參數(shù)學(xué)習(xí)與推理，利用測試集實(shí)現(xiàn)對患者生存情況的預(yù)測，從而評價(jià)該模型的性能。貝葉斯網(wǎng)絡(luò)的參數(shù)學(xué)習(xí)與推理過程均利用貝葉斯網(wǎng)絡(luò)工具箱FullBNT-1.0.4實(shí)現(xiàn)。最終實(shí)驗(yàn)結(jié)果顯示在264個測試集樣本中，預(yù)測正確的達(dá)202例，預(yù)測準(zhǔn)確率為76.52%，表明由MMHC算法構(gòu)建的肺癌預(yù)后模型對肺癌患者5年后生存情況的預(yù)測準(zhǔn)確性良好，可以用于對肺癌患者生存情況的預(yù)測。

3 對比試驗(yàn)

在疾病預(yù)測方面，目前常用的有Logistic回歸、人工神經(jīng)網(wǎng)絡(luò)、決策樹及支持向量機(jī)等機(jī)器學(xué)習(xí)方法[7]。為了進(jìn)一步研究MMHC算法構(gòu)建的貝葉斯網(wǎng)絡(luò)預(yù)后模型的優(yōu)良性，以預(yù)測準(zhǔn)確率為標(biāo)準(zhǔn)，將本模型與Logistic回歸、人工神經(jīng)網(wǎng)絡(luò)、決策樹及支持向量機(jī)等方法在測試集上進(jìn)行分類實(shí)驗(yàn)比較。具體地在WEKA[14]上選擇上述四種方法對應(yīng)的Logistic、J48、Multilayer Perceptron及SMO四個算法，采用十折交叉驗(yàn)證的方法對測試集數(shù)據(jù)進(jìn)行分類，與本算法在預(yù)測準(zhǔn)確率及其他性能指標(biāo)方面作比較，結(jié)果如表3所示。

由表3可知，提出的預(yù)后模型在預(yù)測準(zhǔn)確率、精確度和ROC曲線下面積的結(jié)果均好于其他方法，說明在本研究的肺癌數(shù)據(jù)上貝葉斯網(wǎng)絡(luò)模型是最優(yōu)的。傳統(tǒng)的疾病預(yù)后模型以統(tǒng)計(jì)學(xué)中COX回歸、Logistic回歸為主，但統(tǒng)計(jì)學(xué)方法通常要求變量之間滿足獨(dú)立性等條件，無法處理變量間共線性的問題，因此存在局限性。貝葉斯網(wǎng)絡(luò)模型是一種概率圖模型，通過有向邊和條件概率形象地刻畫出變量間的依賴關(guān)系，能夠進(jìn)行有效地概率推理且預(yù)測準(zhǔn)確率高，可以應(yīng)用于疾病預(yù)測。

表3 不同算法的預(yù)測準(zhǔn)確率及性能指標(biāo)

4 結(jié)論

利用貝葉斯網(wǎng)絡(luò)方法建立肺癌預(yù)后模型，對患者術(shù)后5年生存情況進(jìn)行研究。首先對變量進(jìn)行特征選擇，最終選擇影響患者生存情況的6項(xiàng)預(yù)后因素；然后利用MMHC算法在訓(xùn)練集上建立肺癌預(yù)后模型，在測試集上對患者進(jìn)行5年后生存情況預(yù)測。實(shí)驗(yàn)結(jié)果顯示，利用MMHC算法建立的肺癌預(yù)后模型的預(yù)測準(zhǔn)確率達(dá)76.52%，高于目前常用的Logistic回歸、人工神經(jīng)網(wǎng)絡(luò)、決策樹及支持向量機(jī)方法。但是，本研究未對數(shù)據(jù)集中所有的變量進(jìn)行研究，只是根據(jù)腫瘤信息文獻(xiàn)提取了與生存預(yù)測相關(guān)的16個變量，故研究的模型變量具有一定的主觀性與局限性。在未來的研究中，可以對更多的變量進(jìn)行系統(tǒng)的研究，提高模型的準(zhǔn)確性。