從SHAP到概率*
——可解釋性機(jī)器學(xué)習(xí)在糖尿病視網(wǎng)膜病變靶向脂質(zhì)組學(xué)研究中的應(yīng)用

金東鎮(zhèn) 郭城楠 彭 芳 趙淑珍 李慧慧 夏喆錚 車明珠 王亞楠 張澤杰 毛廣運(yùn),2△

【提 要】 目的 基于可解釋性機(jī)器學(xué)習(xí)算法構(gòu)建糖尿病視網(wǎng)膜病病變(diabetic retinopathy,DR)的早期識(shí)別模型,并探討SHAP(SHapley Additive exPlanations)在脂質(zhì)組學(xué)數(shù)據(jù)中的應(yīng)用。方法 基于本項(xiàng)目組的DR靶向脂質(zhì)組學(xué)數(shù)據(jù),通過可解釋性機(jī)器學(xué)習(xí)的方法進(jìn)行特征篩選;在建立糖尿病視網(wǎng)膜病變的早期識(shí)別模型后,通過全局、特征和個(gè)體三個(gè)層面對(duì)模型進(jìn)行解釋,并將SHAP值轉(zhuǎn)換成概率以增強(qiáng)可解釋的能力。結(jié)果 本研究篩選出了5種內(nèi)源性脂質(zhì)代謝物,構(gòu)建了一個(gè)性能較為優(yōu)秀的糖尿病視網(wǎng)膜病變的早期識(shí)別模型,并成功使用SHAP及概率解鎖了模型。結(jié)論 脂質(zhì)代謝物質(zhì)可以應(yīng)用于糖尿病視網(wǎng)膜病變的早期識(shí)別;SHAP在進(jìn)行黑盒模型的解鎖時(shí)表現(xiàn)出色,且有較高的實(shí)踐應(yīng)用價(jià)值。

糖尿病視網(wǎng)膜病變(diabetic retinopathy,DR)是以視網(wǎng)膜新生血管性增殖為特征的一種致盲性眼病,是糖尿病(diabetes mellitus,DM)最常見、最主要的并發(fā)癥之一,也是勞動(dòng)年齡人口失明的主要原因之一[1-2]。隨著糖尿病的患病率增加、治療手段進(jìn)步、糖尿病患者壽命延長(zhǎng),伴有DR的患者急劇增加,預(yù)計(jì)在2030年全球糖尿病視網(wǎng)膜病變的患者將達(dá)到1.91億人[3]。DR的經(jīng)濟(jì)影響巨大,相比于普通糖尿病患者,DR患者的醫(yī)療花費(fèi)增加了兩倍[4-5],視力障礙和失明也會(huì)對(duì)個(gè)人生活質(zhì)量以及他們所生活的社會(huì)經(jīng)濟(jì)狀況造成毀滅性的影響。早期識(shí)別對(duì)疾病的有效預(yù)防和控制至關(guān)重要且意義重大[6],同時(shí)有研究表明血脂異常與DR的發(fā)生發(fā)展密切相關(guān)[7],可能是DR發(fā)生發(fā)展的早期信號(hào)。而脂質(zhì)組學(xué)(lipidomics)作為代謝組學(xué)的一個(gè)重要分支,已經(jīng)被認(rèn)為是有效識(shí)別潛在脂質(zhì)生物標(biāo)志物的優(yōu)勢(shì)技術(shù)。

隨著機(jī)器學(xué)習(xí)(machine learning,ML)技術(shù)的蓬勃發(fā)展,ML目前已成為高維度大數(shù)據(jù)分析的主流技術(shù)。相較于傳統(tǒng)的組學(xué)分析方法,機(jī)器學(xué)習(xí)方法如隨機(jī)森林(random forest,RF)和支持向量機(jī)(support vector machines,SVM)已經(jīng)在復(fù)雜的組學(xué)數(shù)據(jù)分析中顯現(xiàn)出明顯的優(yōu)勢(shì)[8]。然而正是對(duì)預(yù)測(cè)結(jié)果精確程度的追求,機(jī)器學(xué)習(xí)算法所構(gòu)建的模型往往較為復(fù)雜,犧牲了一定的可解釋性[9],即由于復(fù)雜模型只能給出輸出而常被稱為“黑匣子模型”。為了解鎖復(fù)雜的機(jī)器學(xué)習(xí)模型,SHAP(SHapley Additive exPlanations)作為一種解釋個(gè)體預(yù)測(cè)特征貢獻(xiàn)的方法在2017年被提出[10],并且發(fā)展迅速。簡(jiǎn)單來(lái)說,SHAP解釋方法根據(jù)聯(lián)盟博弈理論計(jì)算某個(gè)個(gè)體不同特征的SHAP值(可以表示貢獻(xiàn),值越大貢獻(xiàn)越大),來(lái)解釋個(gè)體的預(yù)測(cè),其中個(gè)體的特征值充當(dāng)聯(lián)盟中的參與者。在此基礎(chǔ)上,SHAP創(chuàng)新性地將SHAP值的解釋表現(xiàn)為一種可加的特征歸因方法,可表現(xiàn)為:

yi=ybase+f(xi1)+f(xi2)+f(xi3)+…+f(xik)

yi:表示模型的預(yù)測(cè);f(xik):表示第i個(gè)樣本第k個(gè)特征對(duì)最終預(yù)測(cè)的貢獻(xiàn),即對(duì)應(yīng)的SHAP值。

相較于其他模型解釋方法,SHAP源于博弈論,有堅(jiān)實(shí)理論基礎(chǔ);由于其所具有局部準(zhǔn)確性(local accuracy)、缺失(missingness)以及一致性(consistency)等優(yōu)良性質(zhì),可以通過SHAP值聚合獲取模型全局的解釋而廣受歡迎[11]。

本文在課題組前期研究的基礎(chǔ)上,將機(jī)器學(xué)習(xí)方法應(yīng)用于糖尿病視網(wǎng)膜病變的靶向脂質(zhì)組學(xué)數(shù)據(jù)分析中,在構(gòu)建早期識(shí)別模型的同時(shí),使用SHAP解鎖模型并進(jìn)一步將SHAP轉(zhuǎn)換成概率,以獲得其生物學(xué)和醫(yī)學(xué)的解釋,為DR早期診斷提供科學(xué)依據(jù)。

資料與方法

1.數(shù)據(jù)來(lái)源

該數(shù)據(jù)來(lái)自本項(xiàng)目組“基于代謝組學(xué)技術(shù)的糖尿病視網(wǎng)膜病變?cè)缙谧R(shí)別研究”,以糖尿病視網(wǎng)膜病變(DR)者為病例,以2型糖尿病(T2D)患者為對(duì)照,基于年齡、性別、體質(zhì)指數(shù)(body mass index,BMI)、血壓分級(jí)和糖化血紅蛋白(glycosylated hemoglobin,HbA1c),采用傾向性評(píng)分匹配(propensity score matching,PSM)法,按1∶1的比例匹配出69對(duì)研究對(duì)象。所有受試者均接受詳細(xì)的全身檢查,包括身高、體重、坐位血壓(blood pressure,BP)、空腹血糖等常規(guī)檢查以及相關(guān)眼科檢查等;而人口學(xué)指標(biāo)通過詳細(xì)的問卷調(diào)查獲取,問卷均由統(tǒng)一受過標(biāo)準(zhǔn)化操作流程(standardized operation procedures,SOP)培訓(xùn)的研究者對(duì)受試者或其家屬進(jìn)行一對(duì)一的詢問。脂質(zhì)組學(xué)數(shù)據(jù)主要包含了22種本課題組的前期研究中被證實(shí)是差異性脂類代謝物,本項(xiàng)目通過進(jìn)一步的靶向脂質(zhì)組學(xué)分析對(duì)其進(jìn)行分析。

2.統(tǒng)計(jì)分析

統(tǒng)計(jì)分析主要包括數(shù)據(jù)預(yù)處理、特征選擇、模型的構(gòu)建及評(píng)估、模型的解釋四個(gè)方面。

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

缺失值的存在會(huì)對(duì)數(shù)據(jù)挖掘結(jié)果產(chǎn)生較大偏倚,降低統(tǒng)計(jì)分析的效率,因此需要對(duì)缺失數(shù)據(jù)進(jìn)行合理的刪減和填補(bǔ)?？紤]到組學(xué)數(shù)據(jù)缺失的三種原因[12],本研究在刪除缺失比例超過20%的特征以后,一般資料用多重填補(bǔ)的方式進(jìn)行填補(bǔ),而脂類物質(zhì)用最小值的一半進(jìn)行填補(bǔ)[13]。

(2)特征選擇

對(duì)于一般資料來(lái)說,定量資料使用獨(dú)立樣本t檢驗(yàn)或Wilcoxon秩和檢驗(yàn)比較兩組間的差異;定性資料則用卡方檢驗(yàn)或Fisher精確概率法進(jìn)行組間的比較。兩組之間存在統(tǒng)計(jì)學(xué)差異的特征被認(rèn)為可以被納入模型。而在進(jìn)行脂質(zhì)物質(zhì)的篩選時(shí),我們使用了lightGBM模型并通過SHAP獲取各種脂質(zhì)物質(zhì)的重要性排序,隨后按照特征重要性排序?qū)⒅|(zhì)依次進(jìn)入模型并計(jì)算其曲線下面積(area under the curve,AUC),并考慮奧卡姆剃刀原理[14],選取數(shù)量合適的脂質(zhì)物質(zhì)作為特征納入模型。

(3)模型的構(gòu)建

依據(jù)上述方式選擇的一般特征以及脂質(zhì)特征,基于訓(xùn)練集構(gòu)建了DR的早期識(shí)別模型——lightGBM,并進(jìn)行了參數(shù)優(yōu)化。

(4)模型的評(píng)估和解釋

為了驗(yàn)證模型的泛化能力,本研究基于驗(yàn)證集或使用10折交叉驗(yàn)證的方式對(duì)模型的預(yù)測(cè)價(jià)值進(jìn)行判斷。在10折交叉驗(yàn)證中,將數(shù)據(jù)集拆分成為10個(gè)不同的子集,每次用9個(gè)作為訓(xùn)練集,剩余1個(gè)作為測(cè)試集對(duì)模型性能進(jìn)行驗(yàn)證,最后獲得10個(gè)模型性能的平均表現(xiàn)作為10折交叉驗(yàn)證的結(jié)果。另外,我們也使用了HL(Hosmer-Lemeshow)擬合優(yōu)度檢驗(yàn)對(duì)模型的預(yù)測(cè)概率的效果進(jìn)行評(píng)估。最后通過SHAP對(duì)用于早期識(shí)別的lightGBM模型進(jìn)行解釋。然而,相較于概率而言,SHAP值對(duì)那些不了解的人來(lái)說也并不友好,于是我們將SHAP值通過一元插值法轉(zhuǎn)換成了概率并加以展示使其更加方便理解[15]。

結(jié) 果

1.特征篩選

(1)人口學(xué)特征及實(shí)驗(yàn)室指標(biāo)

如表1所示,兩組人群在經(jīng)過傾向性評(píng)分匹配之后,無(wú)論是在訓(xùn)練集還是驗(yàn)證集中,年齡、性別、BMI、空腹血糖(FPG)、糖化血紅蛋白等指標(biāo)在兩組間的差異均未達(dá)到顯著性水平(P>0.05),而糖尿病病程無(wú)論是在訓(xùn)練集或驗(yàn)證集中兩組水平均顯示有差異(P<0.05),而高血壓史在訓(xùn)練集中表現(xiàn)出差異(P<0.05)。因此我們?cè)谌丝趯W(xué)特征及實(shí)驗(yàn)室指標(biāo)中選擇糖尿病病程以及高血壓史進(jìn)入模型。

表1 研究對(duì)象的人口學(xué)和臨床特征等一般情況

圖1A顯示了22種脂質(zhì)物質(zhì)在模型中的重要性排序,22種脂質(zhì)的特征重要性是從上到下依次排列的,其排列的方式主要由模型中某一特征的平均SHAP絕對(duì)值所決定。隨后我們按照?qǐng)D1A所示的特征重要性排續(xù)依次將脂質(zhì)物質(zhì)納入模型,并對(duì)不同數(shù)量的脂質(zhì)物質(zhì)所構(gòu)建模型的AUC進(jìn)行評(píng)估,發(fā)現(xiàn)不管是在驗(yàn)證集或交叉驗(yàn)證中,當(dāng)納入前五種脂質(zhì)時(shí)AUC達(dá)到較為穩(wěn)定的狀態(tài)。再考慮到臨床實(shí)際應(yīng)用和奧卡姆剃刀原理,最后我們選取了前5的特征進(jìn)行模型的構(gòu)建,依次為:OxPC_2[OxPC 34∶2+1O(OxPC 16∶0-18∶2+1O)]、LPG[LPG 18∶1]、FA_3[FA 18∶1]、Acar_1[ACar 16∶2]、Acar_2[ACar 8∶0]。

圖1 脂類物質(zhì)的篩選過程

2.模型構(gòu)建及評(píng)價(jià)

根據(jù)上述的選擇標(biāo)準(zhǔn),我們最終選擇了OxPC_2[OxPC 34∶2+1O(OxPC 16∶0-18∶2+1O)]、LPG[LPG 18∶1]、FA_3[FA 18∶1]、Acar_1[ACar 16∶2]、Acar_2[ACar 8∶0]、糖尿病病程以及高血壓史共7種特征在訓(xùn)練集中構(gòu)建lightGBM模型。為了評(píng)價(jià)已構(gòu)建模型的性能,我們使用了各項(xiàng)評(píng)價(jià)標(biāo)準(zhǔn)對(duì)其進(jìn)行評(píng)估。如圖2所示,在訓(xùn)練集中構(gòu)建的lightGBM模型的ROC曲線下面積為0.92,在測(cè)試集中同樣表現(xiàn)出較好的區(qū)分能力,AUC為0.84。除此之外,如表2所示,模型在訓(xùn)練集中的Accuracy、Precision、Recall、F1-score分別為0.823、0.776、0.918、0.841,而在驗(yàn)證集中分別為0.810、0.800、0.800、0.800,在10折交叉驗(yàn)證中分別為0.861、0.769、0.779、0.846,均顯示出了較為優(yōu)秀的分類能力。并且我們對(duì)模型的預(yù)測(cè)概率進(jìn)行了HL擬合優(yōu)度檢驗(yàn),在訓(xùn)練集中卡方值為10.07(P=0.359),測(cè)試集中卡方值為5.54(P=0.854),表示該模型具有較好的校準(zhǔn)度。

圖2 lightGBM在訓(xùn)練集和驗(yàn)證集中的ROC曲線分析

表2 lightGBM模型性能的評(píng)估

3.模型解釋

(1) 全局的解釋

圖3是對(duì)模型全局的解釋,根據(jù)圖3A所示,左側(cè)縱坐標(biāo)表示了各種特征在模型中的重要性排序,而右側(cè)的黑色柱狀圖則根據(jù)平均SHAP值所得出的每種特征的具體重要程度,由此可見在lightGBM中期重要性依次為:LPG[LPG 18∶1]、OxPC_2[OxPC 34∶2+1O(OxPC 16∶0-18∶2+1O)]、Acar_2[ACar 8∶0]、Acar_1[ACar 16∶2]、FA_3[FA 18∶1]、高血壓史[HYPERTENSION]以及糖尿病病程[DM_duration]。橫坐標(biāo)表示了每一個(gè)個(gè)體,相對(duì)應(yīng)上方的f(x)則是表示每個(gè)個(gè)體的SHAP值,由此可見在LPG[LPG 18∶1]的SHAP值較低時(shí)多數(shù)人相較于平均水平患DR的風(fēng)險(xiǎn)下降。而圖3B不僅特征重要性的同時(shí),可以觀察到每個(gè)特征SHAP值的分布而進(jìn)一步了解特征對(duì)結(jié)局的影響是正相關(guān)還是負(fù)相關(guān)。如LPG[LPG 18∶1]的值越小而SHAP值增大,說明了LPG[LPG 18∶1]的增加可能發(fā)生DR的概率下降,Acar_2[ACar 8∶0]、Acar_1[ACar 16∶2]、FA_3[FA 18∶1]與LPG[LPG 18∶1]相類似,而OxPC_2[OxPC 34∶2+1O(OxPC 16∶0-18∶2+1O)]、高血壓史則相反。

圖3 lightGBM模型的全局解釋

(2) 特征的解釋

上述的圖形只能表示特征對(duì)結(jié)局的影響方向,為進(jìn)一步探究特征如何對(duì)結(jié)果產(chǎn)生影響的,我們繪制了SHAP依賴圖對(duì)某一特征的邊際效應(yīng)進(jìn)行描述。然而由于文章篇幅的限制,本文僅選取了FA_3[FA 18∶1]進(jìn)行展示。如圖4A所示,隨著FA_3[FA 18∶1]的增大,SHAP值降低,也就是說,隨著FA_3[FA 18∶1]的上升,發(fā)生DR的概率下降。然而,正如上文中提到的,相對(duì)于概率來(lái)說SHAP值對(duì)于多數(shù)人來(lái)說并不友好,因此本文在此基礎(chǔ)上,將SHAP值轉(zhuǎn)換成概率(如圖4B所示),相較于A圖來(lái)說,我們可以清楚知道,在FA_3[FA 18∶1]的值為75以下,DR發(fā)生的概率增大,并且在FA_3[FA 18∶1]值達(dá)到100以上時(shí),其邊際效應(yīng)保持在了-0.18%左右,繼續(xù)增大也并未起到較好的效果。

圖4 FA_3[FA 18∶1]在模型中的邊際效應(yīng)

同時(shí),為了探究FA_3[FA 18∶1]與高血壓史[HYPERTENSION]是否存在交互作用,通過對(duì)高血壓的分層,我們繪制了特征交互圖。由圖5可知,有高血壓史的人相較于無(wú)高血壓史的人患有DR的風(fēng)險(xiǎn)較高,而黃色和藍(lán)色兩條曲線并無(wú)相交,說明了FA_3[FA 18∶1]與高血壓史[HYPERTENSION]不存在交互作用。

圖5 FA_3[FA 18∶1]和高血壓史[HYPERTENSION]在模型中的交互

(3)個(gè)體的解釋

我們通過SHAP解釋某一特征的同時(shí),想具體了解某一個(gè)體在模型中的預(yù)測(cè)結(jié)果及其被預(yù)測(cè)正確或錯(cuò)誤的原因,因此我們針對(duì)不同的個(gè)體繪制了SHAP解釋力圖。如圖6所示,A、B、C分別顯示了三個(gè)不同的個(gè)體在模型中的表現(xiàn)。其中圖6A表示了被模型預(yù)測(cè)為低風(fēng)險(xiǎn)的個(gè)體,由于其LPG=0.01492導(dǎo)致其預(yù)測(cè)有DR的概率下降的最大,而OxPC_2=0.1431則增加了這個(gè)個(gè)體被模型判定為發(fā)生DR的概率;而圖6B則是被模型預(yù)測(cè)為高風(fēng)險(xiǎn)的個(gè)體,其中導(dǎo)致他被預(yù)測(cè)為高風(fēng)險(xiǎn)的原因是LPG=0.00669;圖6C則是一個(gè)被模型預(yù)測(cè)為中等風(fēng)險(xiǎn)的個(gè)體。而其中值得一提的是,f(x)在下圖中的解釋是模型的輸出概率的對(duì)數(shù),因此其產(chǎn)生了負(fù)值。

圖6 SHAP解釋力圖:不同的個(gè)體在模型中的表現(xiàn)

討 論

首先,本文在這項(xiàng)基于PSM的多中心病例對(duì)照研究中,成功篩選出5種差異內(nèi)源脂類代謝物質(zhì)與DR發(fā)生風(fēng)險(xiǎn)顯著相關(guān),可作為糖尿病病人中早期識(shí)別DR病人的標(biāo)志物。另外,我們將5種差異性特征代謝物結(jié)合糖尿病病程、糖尿病史構(gòu)建DR早期識(shí)別模型,經(jīng)過模型評(píng)估后驗(yàn)證了其良好的分類效果,為糖尿病視網(wǎng)膜病變的早期識(shí)別提供新思路。其次,本文使用了一種可解釋的機(jī)器學(xué)習(xí)框架對(duì)黑盒模型進(jìn)行解鎖,主要從全局解釋,特征的邊際效應(yīng)以及個(gè)體層面的解釋三個(gè)不同的角度解釋模型,完整地介紹了SHAP。最后,本文在SHAP的基礎(chǔ)上,將SHAP值轉(zhuǎn)換成了概率,不僅進(jìn)一步增強(qiáng)了機(jī)器學(xué)習(xí)的可解釋性,同時(shí)也量化了某一特征對(duì)結(jié)局的影響程度。

近年來(lái),機(jī)器學(xué)習(xí)已經(jīng)被廣泛應(yīng)用于生物醫(yī)學(xué)領(lǐng)域,然而由于其較難解釋的特性,在臨床環(huán)境中的應(yīng)用仍然有限,而SHAP能解鎖黑盒子模型的能力將會(huì)極大推動(dòng)機(jī)器學(xué)習(xí)模型在臨床的實(shí)際應(yīng)用。另外,SHAP中針對(duì)個(gè)體可解釋能力不僅能推動(dòng)模型應(yīng)用于臨床實(shí)際工作中,而且將會(huì)成為精準(zhǔn)醫(yī)療和個(gè)體化醫(yī)療的重要決策依據(jù)。

(致謝:感謝所有項(xiàng)目組成員的支持和幫助,感謝溫州醫(yī)科大學(xué)公共衛(wèi)生學(xué)院預(yù)防醫(yī)學(xué)系代謝組學(xué)研究團(tuán)隊(duì)的奉獻(xiàn)和辛勤工作,感謝2021年浙江省大學(xué)生科技創(chuàng)新活動(dòng)計(jì)劃暨新苗人才計(jì)劃、國(guó)家重點(diǎn)研發(fā)計(jì)劃等項(xiàng)目的支持。)