李建華, 朱澤陽(yáng), 徐禮勝,2, 孫國(guó)哲

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

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

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

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

1 實(shí)驗(yàn)數(shù)據(jù)與方法

1.1 實(shí)驗(yàn)數(shù)據(jù)

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

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

(1)

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

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

圖1 所選特征的缺失情況

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

1.2 K近鄰插值法

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

1.3 棧式自編碼器

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

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

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

(2)

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

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

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

1.4 K-means聚類

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

(3)

其中,μi表示類Ci的聚類中心點(diǎn).

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

1) 在樣本中隨機(jī)選擇K個(gè)點(diǎn)作為聚類中心記作μi,K為指定的聚類簇?cái)?shù).

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

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

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

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

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

3) 輸出聚類簇.

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

深度嵌入聚類網(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)重和每個(gè)隱含層的偏差.

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

(5)

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

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

(6)

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

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

1.6 K-means++初始化

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

1) 隨機(jī)選擇第一個(gè)點(diǎn)作為聚類中心,記作C1.

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

(7)

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

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

2 實(shí)驗(yàn)過程

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

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

(8)

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

圖3 典型的鄰近度矩陣

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

1) 計(jì)算當(dāng)前特征存在缺失的樣本行號(hào),對(duì)應(yīng)鄰近度矩陣D的行號(hào).

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

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

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

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

2.2 評(píng)估指標(biāo)

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

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

(9)

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

(10)

3 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 樣本間差異性分析

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

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

4 結(jié) 語(yǔ)

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

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

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