李建華, 朱澤陽, 徐禮勝,2, 孫國哲

(1. 東北大學 醫(yī)學與生物信息工程學院, 遼寧 沈陽 110169; 2. 沈陽東軟智能醫(yī)療科技研究院有限公司,遼寧 沈陽 110167; 3. 中國醫(yī)科大學附屬第一醫(yī)院 心血管內(nèi)科, 遼寧 沈陽 110001)

重癥監(jiān)護單元(intensive care unit,ICU)是監(jiān)護和搶救危重癥患者的特殊醫(yī)療單元,被稱為生命的最后一道防線[1].自電子病歷出現(xiàn)以來,重癥領(lǐng)域的研究者們可以借助于回顧性的ICU患者生理數(shù)據(jù)開展相關(guān)的研究,例如死亡風險評估、生存期預(yù)測、器官衰竭預(yù)測等[2].這些研究都是依托于電子病歷,因此數(shù)據(jù)質(zhì)量對研究結(jié)果影響很大.

對ICU患者數(shù)據(jù)而言,存在缺失是常見現(xiàn)象.造成數(shù)據(jù)缺失的原因是多樣的,常見的有人為原因和設(shè)備故障等[3].數(shù)據(jù)缺失會造成樣本信息大量損失,進而影響數(shù)據(jù)分析的結(jié)果,因此處理數(shù)據(jù)缺失是數(shù)據(jù)分析任務(wù)的重中之重.目前相關(guān)研究中對缺失值的處理相對簡單,El-Rashidy等[4]在進行ICU患者死亡預(yù)測模型的研究中,將存缺的樣本直接剔除,使用無缺失的樣本構(gòu)建模型,這種做法雖然簡單,但是浪費了大量的可用信息.在更多的相關(guān)研究中[5-7],均值插補法被用于缺失值插補,這種方法確實能減少缺失對樣本的影響,但是如果缺失的樣本較多,會造成樣本間的差異性被縮小.

針對這些問題,本文提出基于深度嵌入聚類構(gòu)造鄰近度矩陣的缺失值插補算法,本算法可以有效地控制用于替代缺失值的樣本數(shù)量,更適用于ICU患者數(shù)據(jù).

1 實驗數(shù)據(jù)與方法

1.1 實驗數(shù)據(jù)

本研究的所有樣本均來自于MIMIC(medical information mart for intensive care)數(shù)據(jù)庫[8].該數(shù)據(jù)庫由貝斯以色列女執(zhí)事醫(yī)療中心、麻省理工學院、牛津大學和麻省總醫(yī)院的急診科醫(yī)生、重癥科醫(yī)生、計算機科學專家等共同建立.本文使用的MIMIC-III V1.4版本包括了2001—2012年期間在貝斯以色列女執(zhí)事醫(yī)療中心重癥監(jiān)護室內(nèi)接受治療的約40 000名患者的數(shù)據(jù),以關(guān)系型數(shù)據(jù)庫表格形式存儲,共計26張數(shù)據(jù)表,包括了人口統(tǒng)計學信息、生理數(shù)據(jù)、精神狀態(tài)、用藥信息、治療方式、患病史等重要數(shù)據(jù).

此23組變量均存在不同程度的缺失,圖1描述了這些特征的缺失情況,圖中缺失率由式(1)定義:

(1)

其中:M_rate為缺失率;M_sample為存在缺失的樣本數(shù)量;All_sample是樣本總數(shù).

本文選取了數(shù)據(jù)庫中23組特征均不存在缺失的5 260例樣本作為研究對象,首先生成多組缺失率不同的人造缺失數(shù)據(jù),然后使用不同方法生成插補后的數(shù)據(jù),比較這些數(shù)據(jù)與原始數(shù)據(jù)的相似性,進而比較插補方法性能.

圖1 所選特征的缺失情況

本文提出一種基于深度嵌入聚類的K近鄰插值算法.本算法以深度嵌入聚類為核心,通過多次聚類構(gòu)造樣本鄰近度矩陣,再自適應(yīng)地選擇缺失樣本的K個近鄰樣本,用這些近鄰樣本的平均值填補缺失.本算法無需手動選擇K值或聚類的簇數(shù),且能有效控制存缺樣本的鄰居數(shù)量,從而有效解決均值插補方法縮小樣本間差異的問題.

1.2 K近鄰插值法

在死亡率預(yù)測的相關(guān)研究中,研究者大多使用均值插補缺失值,因為均值插補計算簡單,效果尚可.均值插補很容易降低樣本間的差異性,而K近鄰數(shù)據(jù)填充法可以彌補這一不足.K近鄰數(shù)據(jù)填充法的核心思想是利用與含缺失值樣本近似的其他樣本的值來填補缺失值,常見方法是先根據(jù)某種距離度量計算得到樣本的“鄰居”,然后用K個近鄰樣本的均值來替代缺失值.求解K近鄰的過程也可以視作聚類的過程,先對數(shù)據(jù)聚類,再在類內(nèi)做插補.

1.3 棧式自編碼器

在解決聚類任務(wù)時,一般步驟為先將高維特征映射到低維空間,再對低維特征使用聚類算法.因為在高維空間中,樣本間的距離難以衡量,即使距離較遠的兩個樣本在某個平面上也可能是近鄰[10],所以降維是聚類前的必要步驟.主成分分析(principal component analysis,PCA)是一種常用的降維方法,其基本思想是在特征空間中找到一條軸,使特征空間的點映射到這條軸上后的方差最大化.PCA作為一種基本的線性降維方法,沒有參數(shù)的限制,大大降低了計算成本,但是學習到的特征較為簡單.

自動編碼器(auto encoder,AE)作為一種深度神經(jīng)網(wǎng)絡(luò),在限制了隱含層的維度后,可以學習到比PCA更全面的特征.由于AE學習到的特征是原特征空間在連續(xù)非交叉曲面上的投影,相比于PCA學習到的低維超平面投影,AE的隱含層表達包含原特征的更多信息[11].在這種考慮下,誕生了深度聚類網(wǎng)絡(luò),即先使用AE學習特征的低維表征,再使用聚類器完成聚類.深度聚類網(wǎng)絡(luò)(deep clustering net,DCN)本質(zhì)上是一個分步模型,降維和聚類是兩個獨立的步驟.深度嵌入聚類則是在DCN的基礎(chǔ)上,使用聚類模型的損失函數(shù)訓練AE的編碼過程,使AE能夠?qū)W習到對聚類任務(wù)友好的低維特征[12].圖2為深度嵌入聚類模型的結(jié)構(gòu),其中Xi表示輸入特征,f(x)表示編碼函數(shù),g(x)表示解碼函數(shù),l(X,Y)是自編碼器輸入和輸出的重構(gòu)誤差函數(shù),f(Xi)是原特征在自編碼器隱含層的表征,也可以看作是降維后的特征,聚類器為K-means.

AE是一種無監(jiān)督學習技術(shù),其訓練過程可以理解為通過編碼函數(shù)對輸入X進行表征學習得到原始特征的編碼,再使用解碼函數(shù)將編碼映射到原特征空間得到重構(gòu)的X′,并使X≈X′.AE相當于重構(gòu)了輸入,其損失函數(shù)為重構(gòu)誤差,如式(2)所示:

(2)

其中:N為樣本個數(shù);Xi表示輸入特征;f(x)表示編碼函數(shù);g(x)表示解碼函數(shù).

因為AE的訓練是最小化重構(gòu)誤差的過程,為了避免其簡單地將輸入復制給輸出,考慮限制隱含層的輸出維數(shù),即限制隱含層的神經(jīng)元個數(shù),使隱含層的輸出維度小于輸入維度,限制各層之間能夠傳遞的信息量,以此強制AE學習到有效的信息[13].而棧式自編碼器(stacked auto encoder,SAE)則是增加隱含層的個數(shù),擴大網(wǎng)絡(luò)容量,使編碼更復雜.

圖2 深度嵌入聚類模型結(jié)構(gòu)圖

1.4 K-means聚類

K-means選用歐氏距離作為相似性度量,其目標函數(shù)是最小化類內(nèi)樣本到聚類中心的距離[14].對于樣本X,假設(shè)將其分為K類,類別記作{C1,C2,…,CK},K-means的目標函數(shù)表示為

(3)

其中,μi表示類Ci的聚類中心點.

K-means最小化目函數(shù)的過程是通過迭代完成的,其步驟為

1) 在樣本中隨機選擇K個點作為聚類中心記作μi,K為指定的聚類簇數(shù).

2) 對于n=1,2,…,N,N為樣本的總個數(shù).

①初始化聚類簇,使Ci=φ.

②計算樣本Xn到每個簇的中心的歐氏距離,將該樣本歸入距離最小的聚類中心所在的類.

③將所有樣本分類完成后,重新計算每個簇的聚類中心,即該簇中所有樣本的均值向量.

④如果聚類中心發(fā)生變化,重復步驟①到③,直至聚類中心不再發(fā)生變化或達到最大迭代次數(shù),結(jié)束迭代.

3) 輸出聚類簇.

1.5 深度嵌入聚類網(wǎng)絡(luò)訓練

深度嵌入聚類網(wǎng)絡(luò)(deep embedded clustering,DEC)將SAE的編碼函數(shù)加入K-means的損失函數(shù)中,編碼函數(shù)如式(4)所示:

hi=f(Xi,W),f(·;W):RM→RR.

(4)

其中:hi是輸入Xi通過多層編碼后的瓶頸層輸出;在此過程中特征由M維降低到R維;W記錄f的全部參數(shù),參數(shù)包括連接各層神經(jīng)元的權(quán)重和每個隱含層的偏差.

將hi作為輸入,更新K-means的目標函數(shù)如式(5)所示:

(5)

其中,Si表示Xi所屬的聚類.

隨后將此目標函數(shù)作為正則項加入SAE的重構(gòu)誤差函數(shù),如式(6)所示:

(6)

其中:λ取1;g(Xi;Z)是解碼函數(shù);Z是由隱含層變換到輸出層的參數(shù);l是自編碼器的重構(gòu)誤差.

SAE的訓練過程如下,首先訓練一個自編碼器,得到隱含層的輸出,再使用第一個自編碼器的隱含層輸出作為輸入訓練第二個自編碼器,最后將所有自編碼器堆疊,得到最終的棧式自編碼器.如圖2所示,本研究使用的SAE共5個隱含層,編碼器的三個隱含層維數(shù)分別設(shè)置為30,20,10,即瓶頸層輸出維度為10.通過預(yù)訓練得到的瓶頸層輸出訓練K-means,得到初始化的聚類簇和聚類中心點,然后在求最優(yōu)聚類簇和聚類中心的過程中同時優(yōu)化參數(shù)W和Z,最終得到對聚類任務(wù)友好的特征表征hi.

1.6 K-means++初始化

K-means++初始化是K-means初始聚類中心選取的優(yōu)化方案[15],根據(jù)K-means的聚類思想,每個樣本到類中心點的距離要盡量近,并且聚類中心點之間要盡量遠,K-means++算法就是考慮選取距離最遠的點作為初始聚類中心.K-means++算法步驟如下:

1) 隨機選擇第一個點作為聚類中心,記作C1.

2) 計算每個樣本Xn到C1的距離,記作D(X),計算每個樣本點被選為下一個聚類中心的概率p(CX),p由式(7)定義:

(7)

3) 累加p,得到每個點被選為下一個聚類中心的概率區(qū)間,然后隨機生成一個0～1之間的數(shù),其所屬的區(qū)間對應(yīng)的點就被選作下一個聚類中心.

4) 重復步驟2)和3),直至選出所有聚類中心,初始化完成.需要注意,當進行步驟2)時已經(jīng)產(chǎn)生了多個聚類中心,則需要計算樣本到每個聚類中心的距離,并選擇其中最小的值作為D(X).

2 實驗過程

2.1 基于鄰近度矩陣的K近鄰插值法

本方法結(jié)合DEC與集成學習思想構(gòu)造鄰近度矩陣.對包含n個樣本的數(shù)據(jù)X,多次使用DEC,因為聚類初始化點是隨機選取的,這樣多次聚類會得到不同的劃分結(jié)果.假設(shè)進行了m次聚類,最終得到m個劃分,定義同類別矩陣N(n×n),Nij=Nji表示樣本i和樣本j在m次聚類中被劃分在同一類中的次數(shù);式(8)定義了鄰近度矩陣D.D度量了一個樣本與其他樣本被多次劃分為同一類的概率,鄰近指數(shù)Dij越接近1,樣本i和j越有可能是近鄰.

(8)

鄰近度矩陣是一個對稱矩陣,圖3給出了矩陣示例,每行表示一個樣本與其他樣本的鄰近程度,數(shù)值越大則兩樣本是近鄰的概率越大,對角線元素為1.

圖3 典型的鄰近度矩陣

將樣本及其特征構(gòu)成的矩陣稱為特征矩陣.對缺失值的插補算法過程如下,對特征矩陣中的每一列特征進行如下操作.

1) 計算當前特征存在缺失的樣本行號,對應(yīng)鄰近度矩陣D的行號.

2) 對D中相應(yīng)的行從高到低排序,取前K個值,計算其對應(yīng)的樣本在特征矩陣中的行號,獲取這些樣本對應(yīng)的特征值.

3) 檢查這K個特征值中是否存在缺失,若存在,用D中第K+1個樣本對應(yīng)的特征值替換,重復此步驟至此K個特征中不含有缺失值.

4) 用此K個特征值的均值代替缺失值.

相比于直接在聚類簇內(nèi)使用均值插補,本方法可以通過參數(shù)K的選取控制用于計算插補值的樣本數(shù)量,相當于控制聚類簇中的樣本數(shù)量,同時保證了每個簇中只有一個缺失值.本方法主要參數(shù)有聚類劃分次數(shù)m,每次劃分的聚類簇數(shù)C以及用于計算缺失特征替代值的近鄰樣本個數(shù)K;根據(jù)集成聚類的相關(guān)研究[16],C取近似于m的值效果較好,一般C=m+1,本文根據(jù)重復實驗選取m值為20.對K值的選取主要依據(jù)每個樣本的高近鄰樣本的個數(shù)來選取,根據(jù)對樣本的觀察,選取鄰近指數(shù)在[0.9, 1]之間的樣本作為鄰近樣本,鄰近樣本個數(shù)K由缺失樣本在規(guī)定區(qū)間內(nèi)的鄰居個數(shù)決定.

2.2 評估指標

本文主要使用插補后數(shù)據(jù)和原數(shù)據(jù)的相似性以及樣本集內(nèi)樣本之間的差異性來衡量插補方法的性能.

使用余弦相似度作為相似性指標,兩個n維向量A和B的余弦相似度Sθ由式(9)定義,Sθ越接近1說明兩個向量越相似,使用特征矩陣行向量余弦相似度的均值表示兩個矩陣的相似度.

(9)

式(10)定義了一個樣本集的總體差異性指標,對一個含有N個樣本集X,每個樣本為Xi(i=1, 2,…,N),其含義為每個樣本與其他樣本的差異性的平均值.

(10)

3 結(jié)果與討論

3.1 數(shù)據(jù)缺失對樣本的影響

為了直觀地展示數(shù)據(jù)缺失對樣本的影響,本文使用PCA方法將原始特征映射到三維空間以便于觀察,將三個維度的特征命名為f1,f2和f3.以此三個特征建立空間坐標系,其中圖4是1 000個無缺失樣本的空間分布,圖5是原始樣本人工加入5%缺失數(shù)據(jù)后在特征空間的分布情況.

圖4 1 000例無缺樣本在三維特征空間的分布圖

如圖5所示,在無缺失原始樣本中人為地添加缺失值后,樣本發(fā)生了肉眼可見的偏移,且樣本更加分散.

圖5 加入5%缺失后1 000例樣本在三維特征空間的分布圖

3.2 插補后與原數(shù)據(jù)的相似度對比

在PCA降維后的三維特征空間中,圖6～圖8分別是均值插補、中值插補和本算法插補后的數(shù)據(jù)在三維特征空間內(nèi)的分布情況.

圖6 使用均值插補法補缺后樣本在三維特征空間的分布圖

由圖6～圖8可知,經(jīng)過插補的數(shù)據(jù)與原數(shù)據(jù)分布近乎一致, 證明插補是一種有效的預(yù)處理手段,但是與本文的方法相比,均值插補和中值插補后的數(shù)據(jù)存在細節(jié)上的不足,一些離群點出現(xiàn)失真.

圖7 使用中值插補法補缺的數(shù)據(jù)在三維特征空間的分布圖

以余弦相似度作為度量,表1比較了均值插補、中值插補、后驗分布估算插補、條件均值插補和基于鄰近度矩陣的插值的性能.從表中可見,本文的方法插值得到的數(shù)據(jù)更接近真實數(shù)據(jù),在缺失率較高時這種優(yōu)勢更加明顯.例如,當缺失率從5%升至15%時,作為對照的4種方法中性能最好的條件均值插補法,余弦相似度從0.993 2降至0.975 4,下降了1.79%;本文方法,相應(yīng)從0.994 8降至0.987 5,僅下降0.73%.

圖8 使用鄰近度矩陣插值法補缺后樣本在三維特征空間的分布圖

表1 使用不同插補方法得出的數(shù)據(jù)與原數(shù)據(jù)的余弦相似度

3.3 樣本間差異性分析

計算可得5 260例樣本的總體差異度為0.869 7,表2給出不同插補方法得出的數(shù)據(jù)的總體樣本間差異度.以均值插補法為例,對于缺失率大于10%的樣本,得到的補缺數(shù)據(jù)樣本間差異度明顯小于原始數(shù)據(jù).當大量樣本存在缺失時,均值插補、中值插補等方法為所有的缺失值都給定相同的替代值,這種計算方法會縮小樣本間差異,對后續(xù)的數(shù)據(jù)分析任務(wù)帶來難度.相比之下,提出的方法在低缺失率時,可以很好地還原數(shù)據(jù),樣本間差異度十分接近原始數(shù)據(jù);而在高缺失率時,樣本差異性的縮小并不明顯,說明本算法很好地還原了數(shù)據(jù)的真實情況.例如,當缺失率從5%升至15%時,4種方法中后驗分布估計插補的總體樣本間差異度變化最小,從0.853 9降至0.835 6,下降了2.14%;本文的方法,相應(yīng)地從0.867 0降至0.860 5,僅下降0.75%,較好地保持了樣本間的差異性.

表2 使用不同插補方法得出的數(shù)據(jù)的總體樣本間差異度

4 結(jié) 語

ICU患者的回顧性生理數(shù)據(jù)中存在缺失值是一種常見現(xiàn)象,在使用數(shù)據(jù)進行統(tǒng)計分析時,缺失值會產(chǎn)生不利影響.缺失值插補是使用估計值代替缺失值,因為ICU患者生理數(shù)據(jù)本身具有變異性,生理參數(shù)值很可能不在正常值范圍內(nèi),對其缺失值的插補是一項具有挑戰(zhàn)性的工作.

本文提出了一種基于深度嵌入聚類構(gòu)造鄰近度矩陣的缺失值插補算法.與ICU患者回顧性數(shù)據(jù)分析的相關(guān)研究中常用的插補方法相比,本算法插補后的數(shù)據(jù)與原始數(shù)據(jù)近似程度較高.本算法能夠有效地確定替代樣本的數(shù)量,根據(jù)樣本間差異性的對比,更好地保留了原樣本的數(shù)據(jù)特性.

此外,本研究還可以在數(shù)據(jù)方面做深入工作.因為MIMIC數(shù)據(jù)庫中不包含中國人的數(shù)據(jù),對國內(nèi)的ICU患者數(shù)據(jù)是否效果較好還需要進一步驗證.未來將與醫(yī)院積極合作,進一步驗證本算法的適用性.