基于機器學習算法探討糖尿病視網(wǎng)膜病變的風險因素

宋亞男，武惠韜，應 俊，李琬悅，陳 康，劉鐵城，張卯年，張 穎

1 解放軍總醫(yī)院醫(yī)學研究創(chuàng)新部 醫(yī)療大數(shù)據(jù)中心，北京 100853；2 解放軍總醫(yī)院醫(yī)療信息處，北京 100853；3 解放軍醫(yī)學院，北京 100853；4 解放軍總醫(yī)院第一醫(yī)學中心 內(nèi)分泌科，北京 100853；5 解放軍總醫(yī)院第一醫(yī)學中心 眼科，北京 100853

糖尿病視網(wǎng)膜病變(diabetic retinopathy，DR)是糖尿病患者主要并發(fā)癥。大約1/3 糖尿病(diabetic mellitus，DM) 患者并發(fā)有不同程度的DR，并且1/10 的DM 患者會因此而視力損害[1]。隨著生活水平的提高，我國糖尿病患病率顯著上升，絕對患病量不容忽視，其中2 型糖尿病(type 2 diabetic mellitus，T2DM) 占90% 以上[2-3]。由于T2DM 患者的臨床數(shù)據(jù)樣本量大、影響因素多、結構復雜，缺失和冗余信息易造成樣本分布不均衡，單純使用傳統(tǒng)統(tǒng)計學方法略顯不足，如欠擬合、準確度下降、不能很好地處理大量多類特征或變量、對非線性特征需轉(zhuǎn)換等[4-5]。近年來，機器學習應用于DR 的人工智能輔助診斷和預測模型被廣泛關注[6-8]。本研究基于解放軍總醫(yī)院大數(shù)據(jù)平臺的DM 并發(fā)癥預警數(shù)據(jù)集，提取臨床觀察變量指標，運用機器學習的方法篩選危險因子，并分析其在DR 形成中的交互作用，結合臨床實際對DM患者預防、控制和管理DR 做出初步探討。

對象與方法

1 研究對象 回顧性分析“國家人口與健康科學數(shù)據(jù)共享平臺”公布的“解放軍總醫(yī)院糖尿病并發(fā)癥預警數(shù)據(jù)集”3 000 例T2DM 患者的臨床資料(2013 年1 月 -2017 年12 月)，其中并發(fā)DR 者1 500 例。提取指標包括年齡、身高、體質(zhì)量等患者基本信息；血壓、血常規(guī)、凝血功能以及血糖等血生化相關指標；常見系統(tǒng)并發(fā)癥情況。

2 數(shù)據(jù)清洗 本研究對得到的數(shù)據(jù)進行了數(shù)據(jù)清洗，去掉與研究無關以及存在過度缺失的變量，并去掉剩余數(shù)據(jù)中存在明顯錯誤或數(shù)據(jù)缺失的研究個體。清洗后共納入2 329 例T2DM，其中合并DR 患者組1 208 例，無DR 并發(fā)癥(non diabetic retinopathy，NDR) 組1121 例；共提取到58 個觀察指標變量。

3 特征選擇 使用遞歸特征消除(recursive feature elimination，RFE) 進行特征選擇，確定納入變量的個數(shù)最優(yōu)解及最終進入模型的預測變量。RFE 方法以整個數(shù)據(jù)集為起點，以預測精度為評價標準，通過每次迭代消除最小相關變量。在此基礎上進行特征排序。變量越相關，排名越高。RFE 方法根據(jù)上述評價準則生成若干特征子集，并最終選擇最優(yōu)特征子集。在本研究中，隨機森林(random forest，RF) 被確定為RFE 的基本分類器，并在訓練集上進行特征選擇。特征篩選的標準是模型優(yōu)化。因此，本研究不涉及變量間的多重共線性。

4 數(shù)據(jù)集劃分 采用隨機抽樣法生成訓練集和驗證集。其中訓練集為70% 的陽性樣本和陰性樣本，以此對預測模型進行訓練。再利用其余樣本構成驗證集，以評估機器學習模型預測糖尿病患者DR 的能力。

5 預測模型訓練與驗證 本研究利用XGBoost開發(fā)預測模型，該算法是在決策樹算法基礎上不斷累加，使精確率不斷提高，達到效果提升的目的。為了使模型更具說服力，我們還比較了XGBoost 與其他兩種流行的機器學習技術：RF 和logistic 回歸。將對應的準確度、敏感度、特異性以及受試者特征曲線下面積(ROC-AUC)等指標作為比較模型性能的標準指標，進行十折交叉驗證比較，并確定總體最佳性能。各標準指標以真陰性(TN) 值，真陽性(TP)、假陰性(FN) 和假陽性(FP)值，根據(jù)下述公式計算之：

準確度=(TP+TN)/(TP+TN+FP+FN)

敏感度=TP/(TP+FN)

特異性=TN/(TN+FP)

XGBoost 輸出只對變量的重要性進行排序，卻無法衡量變量對結果影響的方向和水平。為更好地解釋機器學習模型的結果，對最終模型的危險因子將使用沙普利可加性特征解釋方法(shapley additive explanation，SHAP) 進行解釋。SHAP 是Lundberg 和Lee[5]在2017 年提出的方法，被廣泛用于各種分類、回歸模型的解釋。對于每個預測樣本，模型都產(chǎn)生一個預測值，SHAP 值就是該樣本中每個特征所分配到的數(shù)值，其絕對值反映特征影響力大小，正負反映其正向或負向作用。當SHAP 值＞0 時，表示特征改善了結果值，具有積極的效果；相反，當SHAP 值 <0 時，表明該特征降低了結果值，并產(chǎn)生了負面影響。該方法可以輸出特征的重要性排序，以及特征與結果之間的關系，從而解釋變量怎樣影響預測結果。

6 統(tǒng)計學方法 對變量進行基線描述分析，其中正態(tài)分布的數(shù)值變量以表示，非正態(tài)分布的數(shù)值變量使用中位數(shù)(第一四分位數(shù)，第三四分位數(shù))表示；而分類變量使用頻率(百分比)的形式進行描述。對變量進行差異性假設檢驗，其中分類變量使用χ2檢驗；對呈正態(tài)分布、方差齊的數(shù)值變量使用t檢驗；對不同時滿足正態(tài)分布和方差齊的變量使用了秩和檢驗。以α=0.05 為檢驗水準。本研究在Oracle 數(shù)據(jù)庫上使用程序語言/SQL 進行數(shù)據(jù)檢索。統(tǒng)計分析采用R 程序語言(3.6.1 版)和Python 語言(3.7.7 版)。

結果

1 納入病例兩組間臨床資料比較 糖尿病合并DR 組的收縮壓、舒張壓、糖化血紅蛋白(glycosylated hemoglobin，HbA1c)、總膽固醇、低密度脂蛋白、纖維蛋白原、血尿素、血肌酐、血清尿酸、乳酸脫氫酶和凝血酶原活動度指標均高于NDR 組(P＜ 0.05)；而年齡、血紅蛋白、血細胞比容、球蛋白、總蛋白、血清白蛋白、總膽紅素、直接和間接膽紅素、谷丙轉(zhuǎn)氨酶、谷草轉(zhuǎn)氨酶、谷氨酰胺轉(zhuǎn)移酶、堿性磷酸酶、凝血酶原時間和部分活化凝血酶原時間低于NDR 組(P＜0.05)。見表1。

表1 兩組糖尿病患者臨床資料比較分析Tab.1 Comparison of baseline variables between the two groups

2 兩組基礎疾病比較 DR 組高血壓、腦卒中、脂肪肝、膽道疾病、腎衰、腎病、下肢動脈病變、血液病和其他內(nèi)分泌疾病等的并發(fā)比例高于NDR 組(P＜ 0.05)；而高脂血癥、動脈粥樣硬化、冠心病、心肌梗死和風濕免疫疾病等低于NDR 組(P＜ 0.05)。DR 組排名前6 位的系統(tǒng)病癥依次為高血壓、腎病、動脈粥樣硬化、脂肪肝、下肢動脈病變、冠心??；NDR 組排名前6 位的系統(tǒng)病癥依次為高血壓、動脈粥樣硬化、冠心病、高脂血癥、脂肪肝、腎病。見表2。

表2 兩組糖尿病患者合并其他系統(tǒng)病癥情況比較(n,%)Tab.2 Comparison of patients with other systemic diseases between the two groups (n,%)

3 建立模型 與其他模型相比，XGBoost 算法的準確率(0.961)、F1 分值(0.783)、精確度(0.783)、召回率(0.783)和ROC-AUC(0.842)表現(xiàn)良好(表3)。對其產(chǎn)生的特征重要性(特征權重值) 進行排序，前20 位的變量依次為腎病、冠心病、下肢動脈病變、身高、其他腫瘤、HbA1c、血尿素、血清白蛋白、腎衰、高脂血癥、消化系統(tǒng)腫瘤、血肌酐、年齡、神經(jīng)系統(tǒng)疾病、谷氨酰胺轉(zhuǎn)移酶、直接膽紅素、谷草轉(zhuǎn)氨酶、谷丙轉(zhuǎn)氨酶、其他內(nèi)分泌疾病、總蛋白等(圖1)。這些特征對最終的模型影響較大，但是無法準確描述特征與DR 最終預測結果的關系。

圖1 基于XGBoost 算法的糖尿病視網(wǎng)膜病變風險前20 位特征重要性排序Fig.1 The importance of the top 20 features on diabetic-retinopathy-risk based on XGBoost algorithm

表3 三種模型預測指標比較Tab.3 Comparison of model training

4 SHAP 模型分析 在特征重要性排序基礎上，利用SHAP 方法對最終模型進一步分析，各變量的SHAP 絕對值權重排序見圖2。根據(jù)圖3 中SHAP 模型在訓練集中生成的散點圖，可結合臨床危險因子的意義進行分析。其中，糖化血紅蛋白(0.59)、腎病(0.44)、血尿素(0.32)、下肢動脈病變(0.25) 四項觀察指標SHAP 值的絕對值均高，在SHAP 值的分布上呈現(xiàn)明顯的分類狀況，SHAP值 ＞ 0，其特征提升了DR 預測值，即為DR 的顯著危險因素。

圖2 SHAP 方法的特征權重排序Fig.2 Feature weight sorting by SHAP method

圖3 SHAP 分析的變量散點圖Fig.3 SHAP summary plot of XGBoost model

5 HbA1c、血尿素、腎病三個變量分析 SHAP 分析所得出的特征重要性較高的三個變量。與HbA1c 最強相關的因素為血尿素，且隨HbA1c 的值升高(圖4)，該因素對模型呈現(xiàn)陽性結果的作用越來越明顯。當HbA1c 在5%～ 8% 時，該因素對模型結果呈現(xiàn)陰性影響。在HbA1c ＜ 5% 時，血尿素的值均高于正常水平。血尿素與腎病如圖5所示，尿素指標大于5 mmol/L 時，患者患腎病的比例明顯增多，且隨著尿素的增高，對模型呈現(xiàn)陽性結果的影響逐漸增大。對血尿素、HbA1c 與是否患DR 做3D 散點圖，所得結果如圖6 所示，血尿素值過高的病例幾乎全部出現(xiàn)在DR 群體中，而HbA1c 在兩者間的分布趨勢無明顯差異。

圖5 血尿素與腎病交互的SHAP 值分布圖Fig.5 Scatter plot of SHAP value about interaction between serum urea and nephropathy

圖6 糖化血紅蛋白-血尿素-糖尿病視網(wǎng)膜病變 3D 散點圖Fig.6 3D scatter plot of SHAP value about HbA1c-serum urea-diabetic retinopathy (NDR:no diabetic retinopathy;DR:diabetic retinopathy)

討論

回顧既往相關研究，糖尿病病程、血糖控制不佳以及高血壓被認為是DR 發(fā)生的主要風險因素[9-10]。這在理論及臨床實踐中被廣泛認可。但實際工作中會發(fā)現(xiàn)血糖和(或)血壓的異常并不一定導致DR 的發(fā)生，即有些血糖和或血壓控制差的DM 患者并不發(fā)生DR，而部分血糖控制好且無高血壓的患者反而會發(fā)生DR，提示還有其他因素影響了DR 的發(fā)生發(fā)展。研究表明，代表血糖水平的HbA1c 僅能解釋至多11%的DR 風險性，并不是絕對和唯一影響要素，變異性較大的另89%風險歸因于糖尿病患者的其他整體因素[10]。近來研究結果指出，早在臨床可見眼底DR 之前，視網(wǎng)膜血管和神經(jīng)的病理性改變就已長期存在[11]。因此，對于探究DR 發(fā)病的個體差異和原因，分析和明確DM 患者其他更多相關的系統(tǒng)影響因素，對于早期有效預防和治療DR、降低DM 患者視功能喪失風險十分必要。

根據(jù)本研究模型結果，HbA1c 為兩組間呈現(xiàn)顯著差異且與DR 高度相關的變量，與過往研究報道一致[12-14]。另外，本研究顯示DR 組患者的HbA1c 水平高于NDR 組，提示DR 患者的血糖控制水平普遍較差，因此也更易出現(xiàn)糖尿病腎病等其他微血管并發(fā)癥。英國前瞻性DM 研究長期隨訪(≤ 20 年) 顯示，早期血糖控制對大血管并發(fā)癥的有益作用是確認的，HbA1c 目標值 ＜ 7%(＜53 mmol/mol)可減少微血管并發(fā)癥[15]。另有研究指出HbA1c 48 mmol/mol(6.5%)的閾值可用來明確診斷糖尿病和DR[16]。

合并腎病的T2DM 患者DR 患病風險增高，為本研究風險因素分析的結果，與其他文獻報道一致[13,17-18]。一項前瞻性研究證實DM 合并腎病的患者進行多因素干預具有益效果，即通過對腎病和DM 的聯(lián)合治療可以顯著降低微血管和大血管事件的發(fā)生率[19]。腎與視網(wǎng)膜在起源、發(fā)育、毛細血管網(wǎng)結構、濾過屏障功能方面有高度相似性[20]。DR 與糖尿病腎病同為DM 的微血管并發(fā)癥，兩者之間的發(fā)生發(fā)展具有某種程度的相關性，如均受到遺傳、血流動力學、脂質(zhì)代謝等共同影響，發(fā)病機制都包括了糖基化終末產(chǎn)物的堆積、多元醇通路活化、氧化應激、生長因子、炎癥介質(zhì)、補體激活等[21]。目前，DR 伴微量白蛋白尿即可作為DM 引起腎病的診斷標準，DR 在糖尿病腎病的臨床輔助診斷中的重要價值和意義已達成專家共識[22]。還有研究指出DR 嚴重程度可用來預測T2DM 相關腎病，但DR 與糖尿病腎病的嚴重程度可不一致[23]。

本研究中，高血壓和動脈粥樣硬化在NDR 和DR 兩組均呈現(xiàn)大于50%的發(fā)病率，尤其高血壓在NDR 和DR 組患病率更分別高達66.5%和73.8%；而高血壓、高脂血癥、腦卒中、脂肪肝、冠心病、心肌梗死、血液病等血循環(huán)系統(tǒng)異常、脂質(zhì)代謝疾病在兩組間患病率差異具有統(tǒng)計學意義，均在一定程度上表明這些因素與DM 和DR 病程進展的機制有關。有研究報道稱，隨DM 病程延長，機體病變逐漸由血管內(nèi)物質(zhì)的轉(zhuǎn)變發(fā)展至小血管病變甚至演變成為大血管病變[24]，本研究所示的腎病以及下肢動脈病變患病率在DR 患者中明顯增高且在模型中預測重要性排名靠前，均印證了這一觀點。

血尿素是機體蛋白質(zhì)代謝的主要終末產(chǎn)物，從腎小球濾過而排出體外。有研究表明，T2DM 患者24 h 血尿素與DR 的發(fā)生相關[25]。這與本研究結果一致。血尿素正常值范圍為2.86～ 7.14 mmol/L，器質(zhì)性腎功能損害可使血尿素增高。本研究結果顯示，血尿素指標＞5 mmol/L 時，T2DM 患者患腎病的比例明顯升高，且隨著血尿素水平的增高，對模型呈現(xiàn)陽性結果的影響逐漸增大，意味著DR 的患病風險越高，與以往研究一致[26]。故血尿素水平除作為腎功能的觀察指標，應作為DM 患者除血糖水平之外需重點關注的一項檢驗學觀測指標，通過飲食起居、藥物治療等調(diào)節(jié)和維持血尿素正常水平以預防DM 引發(fā)更多組織并發(fā)癥和功能損傷。

本研究結果顯示，HbA1c、血尿素、腎病這三個SHAP 分析所得出的特征重要性較高的變量對DR 病程影響存在潛在交互關系。血尿素與HbA1c 強相關，意味著在血糖水平控制良好的情況下，較高的血尿素指標會導致模型呈現(xiàn)陽性結果。血尿素值過高的患者幾乎全部出現(xiàn)在DR 組的患者群體中，而HbA1c 在兩組間的分布趨勢無明顯差異，此時兩組間的差異可能是由于血尿素的差異所致。故在對HbA1c 檢測的同時，對血尿素高度關注可能能夠為DM 病程管理提供更全面的建議。此結果亦為在HbA1c 控制較為良好時，部分患者仍患有DR 的病理機制提供新的研究方向。

綜上所述，本研究通過機器學習算法建立優(yōu)選模型，并利用SHAP 分析方法分別解釋了各危險因素對DR 結局的影響及臨床意義，找到T2DM患者并發(fā)DR 的主要風險因素，即HbA1c、血尿素水平和是否伴腎病，且三者對DR 病程影響具有潛在交互關系。由此給予臨床提示，一方面通過監(jiān)測HbA1c 水平控制血糖，另一方面血尿素測量可作為DM 微血管并發(fā)癥一級預防高度關注的提示性指標，并采取干預措施，可為DM 患者的長期病程管理提供更全面的建議，預防和降低DR 患病風險，也可為未來的相關研究提供方向。作為單中心臨床研究，本文數(shù)據(jù)的代表性具有局限性。因此，進一步開展多中心、信息完備的大樣本數(shù)據(jù)采集為探討真實世界糖尿病視網(wǎng)膜病變發(fā)病影響因素的必要途徑和努力方向。

致謝：解放軍總醫(yī)院醫(yī)療大數(shù)據(jù)中心喬屾工程師、楊明工程師為本研究提供的指導和幫助。