從SHAP到概率*
——可解釋性機器學習在糖尿病視網膜病變靶向脂質組學研究中的應用

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

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

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

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

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

yi:表示模型的預測;f(xik):表示第i個樣本第k個特征對最終預測的貢獻,即對應的SHAP值。

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

本文在課題組前期研究的基礎上,將機器學習方法應用于糖尿病視網膜病變的靶向脂質組學數據分析中,在構建早期識別模型的同時,使用SHAP解鎖模型并進一步將SHAP轉換成概率,以獲得其生物學和醫學的解釋,為DR早期診斷提供科學依據。

資料與方法

1.數據來源

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

2.統計分析

統計分析主要包括數據預處理、特征選擇、模型的構建及評估、模型的解釋四個方面。

(1)數據預處理

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

(2)特征選擇

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

(3)模型的構建

依據上述方式選擇的一般特征以及脂質特征,基于訓練集構建了DR的早期識別模型——lightGBM,并進行了參數優化。

(4)模型的評估和解釋

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

結 果

1.特征篩選

(1)人口學特征及實驗室指標

如表1所示,兩組人群在經過傾向性評分匹配之后,無論是在訓練集還是驗證集中,年齡、性別、BMI、空腹血糖(FPG)、糖化血紅蛋白等指標在兩組間的差異均未達到顯著性水平(P>0.05),而糖尿病病程無論是在訓練集或驗證集中兩組水平均顯示有差異(P<0.05),而高血壓史在訓練集中表現出差異(P<0.05)。因此我們在人口學特征及實驗室指標中選擇糖尿病病程以及高血壓史進入模型。

表1 研究對象的人口學和臨床特征等一般情況

圖1A顯示了22種脂質物質在模型中的重要性排序,22種脂質的特征重要性是從上到下依次排列的,其排列的方式主要由模型中某一特征的平均SHAP絕對值所決定。隨后我們按照圖1A所示的特征重要性排續依次將脂質物質納入模型,并對不同數量的脂質物質所構建模型的AUC進行評估,發現不管是在驗證集或交叉驗證中,當納入前五種脂質時AUC達到較為穩定的狀態。再考慮到臨床實際應用和奧卡姆剃刀原理,最后我們選取了前5的特征進行模型的構建,依次為: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 脂類物質的篩選過程

2.模型構建及評價

根據上述的選擇標準,我們最終選擇了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種特征在訓練集中構建lightGBM模型。為了評價已構建模型的性能,我們使用了各項評價標準對其進行評估。如圖2所示,在訓練集中構建的lightGBM模型的ROC曲線下面積為0.92,在測試集中同樣表現出較好的區分能力,AUC為0.84。除此之外,如表2所示,模型在訓練集中的Accuracy、Precision、Recall、F1-score分別為0.823、0.776、0.918、0.841,而在驗證集中分別為0.810、0.800、0.800、0.800,在10折交叉驗證中分別為0.861、0.769、0.779、0.846,均顯示出了較為優秀的分類能力。并且我們對模型的預測概率進行了HL擬合優度檢驗,在訓練集中卡方值為10.07(P=0.359),測試集中卡方值為5.54(P=0.854),表示該模型具有較好的校準度。

圖2 lightGBM在訓練集和驗證集中的ROC曲線分析

表2 lightGBM模型性能的評估

3.模型解釋

(1) 全局的解釋

圖3是對模型全局的解釋,根據圖3A所示,左側縱坐標表示了各種特征在模型中的重要性排序,而右側的黑色柱狀圖則根據平均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]。橫坐標表示了每一個個體,相對應上方的f(x)則是表示每個個體的SHAP值,由此可見在LPG[LPG 18∶1]的SHAP值較低時多數人相較于平均水平患DR的風險下降。而圖3B不僅特征重要性的同時,可以觀察到每個特征SHAP值的分布而進一步了解特征對結局的影響是正相關還是負相關。如LPG[LPG 18∶1]的值越小而SHAP值增大,說明了LPG[LPG 18∶1]的增加可能發生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) 特征的解釋

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

圖4 FA_3[FA 18∶1]在模型中的邊際效應

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

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

(3)個體的解釋

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

圖6 SHAP解釋力圖:不同的個體在模型中的表現

討 論

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

近年來,機器學習已經被廣泛應用于生物醫學領域,然而由于其較難解釋的特性,在臨床環境中的應用仍然有限,而SHAP能解鎖黑盒子模型的能力將會極大推動機器學習模型在臨床的實際應用。另外,SHAP中針對個體可解釋能力不僅能推動模型應用于臨床實際工作中,而且將會成為精準醫療和個體化醫療的重要決策依據。

(致謝:感謝所有項目組成員的支持和幫助,感謝溫州醫科大學公共衛生學院預防醫學系代謝組學研究團隊的奉獻和辛勤工作,感謝2021年浙江省大學生科技創新活動計劃暨新苗人才計劃、國家重點研發計劃等項目的支持。)