南京醫(yī)科大學公共衛(wèi)生學院生物統(tǒng)計學系（211166） 董學思 林麗娟 趙 楊 魏永越 戴俊程 陳 峰

多組學聯(lián)合缺失數(shù)據(jù)填補方法的評價*

南京醫(yī)科大學公共衛(wèi)生學院生物統(tǒng)計學系（211166） 董學思 林麗娟 趙 楊 魏永越 戴俊程 陳 峰△

目的 本研究旨在評價不同平臺間“塊缺失”數(shù)據(jù)的填補方法。如何在保證方差-協(xié)方差結(jié)構(gòu)相對穩(wěn)定的前提下提高多組學數(shù)據(jù)填補的精確度，對于后期數(shù)據(jù)挖掘有重要的意義。方法 利用癌癥基因組圖譜（TCGA）數(shù)據(jù)庫的肺癌數(shù)據(jù)（甲基化數(shù)據(jù)、基因表達數(shù)據(jù)），構(gòu)建不同缺失比例的數(shù)據(jù)集（缺失比例分別為5%、20%、35%、50%和65%）。采用統(tǒng)計學填補方法均值法，馬爾科夫蒙特卡洛法（MCMC）和機器學習填補法［鄰近法（kNN），隨機森林法（RF），多層感知機法（MLP）］對缺失數(shù)據(jù)進行填補，填補后數(shù)據(jù)集與原數(shù)據(jù)集進行比較。評價指標包括估計偏差和矩陣-2-范數(shù)。根據(jù)評價指標和填補時間，比較出填補效果最優(yōu)、填補時間較短的方法。結(jié)果 MLP和kNN算法在各種缺失比例下均比其他填補方法有更優(yōu)的效果，填補時間也相對較短。均值法的時間最短，在數(shù)據(jù)集缺失比例較小時（≤5%），填補效果與其他填補方法相當，但在高比例缺失情況下表現(xiàn)較差。在數(shù)據(jù)集高比例缺失情況下，RF和MCMC的填補效果優(yōu)于均值法，但填補時間過長，不適用于實際工作。結(jié)論 綜合比較，機器學習填補方法中的MLP和kNN兩法適合于甲基化數(shù)據(jù)和表達數(shù)據(jù)的填補。

多組學數(shù)據(jù) 塊缺失 統(tǒng)計學填補 機器學習填補 效果評價

生物數(shù)據(jù)的獲取受限于現(xiàn)階段技術(shù)手段所存在的不足（測序過程中對比基因組測序誤差、芯片劃痕、圖像污染等），缺失數(shù)據(jù)的產(chǎn)生不可避免［1］。高維生物數(shù)據(jù)研究中，由于樣本數(shù)據(jù)往往存在不同平臺的測量信息，在樣本信息匹配時，經(jīng)常存在有部分樣本缺失某平臺數(shù)據(jù)的情況——“塊缺失”。

傳統(tǒng)上，研究者在探索缺失值填補方面更側(cè)重于統(tǒng)計填補法，該類方法雖簡單易行，卻難以挖掘數(shù)據(jù)之間的深層關(guān)系，填補效果并不理想［2］；國外近幾年在機器學習領(lǐng)域的填補方法研究較多，該類方法對數(shù)據(jù)分布類型不敏感，通過計算機模擬的形式深入挖掘數(shù)據(jù)結(jié)構(gòu)關(guān)系，填補效果較優(yōu)［4］?，F(xiàn)階段，數(shù)據(jù)填補的主要手段均基于以上兩大類方法（統(tǒng)計填補、機器學習填補），國內(nèi)研究者進行過一系列統(tǒng)計填補方法的效果評價，但對統(tǒng)計填補方法和機器學習填補方法的比較并不多見。

本研究通過構(gòu)建不同缺失比例的數(shù)據(jù)集，采用統(tǒng)計填補法和機器學習填補法進行填補，采用估計偏差和矩陣-2-范數(shù)作為評價指標，比較填補效果，選出優(yōu)勢方法。不同平臺數(shù)據(jù)之間的填補將會提高信息利用率，提高檢驗效能，有助于研究者得到更可靠的結(jié)果。

數(shù)據(jù)與方法

1．數(shù)據(jù)

選擇癌癥基因組圖譜（TCGA）公共數(shù)據(jù)庫的肺癌甲基化、基因表達數(shù)據(jù)，將含有缺失變量的樣本剔除，保留完整樣本，便于填補效果比對。選擇肺癌經(jīng)典生物通路WNT通路中基因表達變量141個，甲基化位點3962個，樣本782例，作為第一個研究數(shù)據(jù)集［5］；再按照同樣的變量類型和變量數(shù)，從全基因組中隨機抽取141個基因表達變量和3962個甲基化位點，保留完整樣本880例，作為第二個數(shù)據(jù)集。將兩個數(shù)據(jù)集按5%、20%、35%、50%、65%的樣本缺失比例構(gòu)造缺失，缺失部分樣本的選擇為完全隨機，避免偶然性。其中，缺失的樣本中50%樣本的表達數(shù)據(jù)缺失，50%樣本的甲基化數(shù)據(jù)缺失，即在不完整觀測中，缺失甲基化數(shù)據(jù)的樣本擁有完整的基因表達數(shù)據(jù)；缺失基因表達數(shù)據(jù)的樣本擁有完整的甲基化數(shù)據(jù)。這樣的缺失比例更加符合“塊缺失”的數(shù)據(jù)情況。在接下來的填補工作中，筆者將用甲基化數(shù)據(jù)對表達數(shù)據(jù)進行填補，用表達數(shù)據(jù)對甲基化數(shù)據(jù)進行填補。

本研究利用這兩個數(shù)據(jù)集，保證在相關(guān)結(jié)構(gòu)和非相關(guān)結(jié)構(gòu)數(shù)據(jù)集中填補方法的穩(wěn)定性，這樣構(gòu)建的相關(guān)結(jié)構(gòu)更貼近實際情況。將各種缺失比例的數(shù)據(jù)集進行5種填補方法的填補，每種方法填補100次，將填補后數(shù)據(jù)集與原數(shù)據(jù)集進行對比，并計算綜合評價指標，評價填補效果。模擬試驗流程見圖1。

圖1 模擬試驗流程

2．方法

此次研究采用統(tǒng)計填補法和機器學習填補法：分別是均數(shù)法（mean）、馬爾科夫蒙特卡洛法（Markov Chain Monte Carlo，MCMC）、隨機森林法（random forest，RF）、K鄰近法（k-nearest neighbor，kNN）和多層感知機法（multi-layer perceptron，MLP）。均值法是最經(jīng)典的填補方法之一，簡便易行。MCMC填補在之前研究中被廣泛提及，因其能充分利用完整數(shù)據(jù)部分作為先驗，因此填補的效果往往高于常見的統(tǒng)計方法填補。機器學習的三種填補方法在單一數(shù)據(jù)集填補的研究中有過報道，但尚未應(yīng)用于多組學數(shù)據(jù)的填補中。

（1）統(tǒng)計填補

①均值法（mean）

采用缺失變量非缺失部分的均值對缺失部分進行填補，填補方法簡單，填補時間極短，但未考慮到數(shù)據(jù)本身變異性，降低填補后數(shù)據(jù)方差，破壞原有數(shù)據(jù)結(jié)構(gòu)。其步驟可以整理為：計算缺失變量中非缺失部分的均值；所計算均值代替缺失數(shù)據(jù)。

②馬爾科夫蒙特卡洛填補（MCMC）

MCMC利用變量均值向量和方差-協(xié)方差陣作為先驗信息，構(gòu)建馬爾科夫鏈，保證其元素的分布可以收斂到一個平穩(wěn)分布，通過抽樣反復模擬該馬爾科夫鏈，得到平穩(wěn)的后驗分布，產(chǎn)生缺失數(shù)據(jù)的估計［6］。其步驟可以整理為：

a將數(shù)據(jù)集拆分為完整觀測部分Xfull和不完整觀測部分Xmiss。

μ＝［μ′1，μ′2］代表Xfull和Xmiss的均值向量；∑11、∑22分別代表Xfull、Xmiss的方差-協(xié)方差矩陣，∑12代表Xfull和Xmiss的方差-協(xié)方差矩陣；

b給定Xfull＝X1，Xmiss的均向量為μ2．1＝μ2+∑′12∑′11（X1-μ1）；相應(yīng)條件協(xié)方差矩陣為∑22．1＝∑22-∑′12∑-1

11∑12；c給定Xfull，從Xmiss的條件分布中隨機抽取數(shù)值，對缺失部分進行填補；

d經(jīng)過填補后，產(chǎn)生完整數(shù)據(jù)集，循環(huán)上述步驟，估算新產(chǎn)生的均向量和協(xié)方差矩陣進行填補，直至收斂。

（2）機器學習填補法

①隨機森林填補（RF）

RF填補應(yīng)用集成決策樹的思維，完全分裂產(chǎn)生回歸樹，每一棵分類樹代表一個多元非線性模型，產(chǎn)生缺失變量的加權(quán)平均值，對缺失數(shù)據(jù)進行填補［3］。其步驟可以整理為：

a將完整數(shù)據(jù)部分作為自變量，缺失數(shù)據(jù)部分作為預測變量；

b在數(shù)據(jù)集中采用Bagging的方法，隨機抽取部分的樣本作為單棵決策樹的訓練集；

c按照完全分裂構(gòu)造決策樹回歸器，每棵樹產(chǎn)生一批填補值，最后將各棵樹的結(jié)果取平均值，作為填補值進行填補。

②K鄰近填補（kNN）

kNN填補在樣本數(shù)據(jù)集的特征空間中，按照馬氏距離選取相近（即特征空間中鄰近）的樣本集，計算對應(yīng)變量的加權(quán)平均值進行填補。相較均值填補，kNN考慮了樣本間的變異，保持了數(shù)據(jù)結(jié)構(gòu)的穩(wěn)健性［7］。其步驟可以整理為：

a構(gòu)建完整樣本集的矩陣結(jié)構(gòu)；

b計算含有缺失變量的樣本集Xmiss與完整樣本集Xfull中各樣本的馬氏距離；

d所計算的均值代替缺失數(shù)據(jù)。

③多層感知機填補（MLP）

MLP是人工神經(jīng)網(wǎng)絡(luò)的重要分支，通過訓練集樣本訓練神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，經(jīng)多次層間映射，產(chǎn)生缺失變量估計值［8］。MLP尤其適用于混合分布數(shù)據(jù)庫，在高維多平臺數(shù)據(jù)中，可以綜合不同平臺信息訓練構(gòu)建人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，進行缺失數(shù)據(jù)填補。本次研究采用標準的三層單向神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)：第一層為輸入層，輸入某樣本完整部分的變量；第二層為隱藏單元；第三層為輸出層，產(chǎn)生缺失數(shù)據(jù)的估計值，見圖2。層間由權(quán)重矩陣連接，輸入層經(jīng)過隱藏單元層映射至輸出層，產(chǎn)生填補值。其步驟可以整理為：

圖2 多層感知機網(wǎng)絡(luò)結(jié)構(gòu)

a構(gòu)建完整數(shù)據(jù)集矩陣Xfull，作為訓練數(shù)據(jù)集，含有缺失變量的樣本作為預測集Xmiss；

b采用剪枝算法，交叉驗證，計算隱藏單元數(shù)目；

c采用共軛梯度法，計算層間權(quán)重向量；

d根據(jù)輸入向量，通過映射函數(shù)映射至隱藏單元：

zh是隱藏單元，h＝1，…，h，hj是第一層權(quán)重矩陣，who是殘差項，f（）為激勵函數(shù)，通常為雙曲正切函數(shù)或logit函數(shù)；

e隱藏單元再經(jīng)過一次激勵函數(shù)轉(zhuǎn)化至輸出層。yk是輸出單元k＝1，…，k，wkh是第二層權(quán)重矩陣，wko是殘差項，g（）為線性激勵函數(shù)：

f產(chǎn)生神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)；

g輸入預測集樣本Xmiss，估計預測集中缺失數(shù)據(jù)。

5種方法的填補效果用估計偏差和矩陣-2-范數(shù)來評價。估計偏差：原數(shù)據(jù)集中變量均數(shù)，與填補后數(shù)據(jù)集中變量的均數(shù)之差的絕對值之和，再取平均值?？啥x為：

其中，x-ori為原數(shù)據(jù)集中各變量的均值，x-imp為填補后的數(shù)據(jù)集中各變量的均值，p為變量數(shù)。該指標反應(yīng)填補的精確度，估計偏差越小，填補的精確度越高。

矩陣-2-范數(shù)：轉(zhuǎn)置矩陣d’與原矩陣d的積的最大特征根的平方根值。幾何意義指空間上兩個矩陣（向量）的距離。待比較的兩個矩陣作差得矩陣d，求得矩陣-2-范數(shù)，反映的是差值矩陣距離原點的距離，即：方差-協(xié)方差矩陣變化幅度，可定義為：

矩陣d為填補后數(shù)據(jù)集的方差-協(xié)方差矩陣與原數(shù)據(jù)集方差協(xié)方差矩陣的差，eigen（）函數(shù)分解矩陣特征根，max（）函數(shù)求得最大特征根。該指標反映的是填補數(shù)據(jù)集與原數(shù)據(jù)集的數(shù)據(jù)結(jié)構(gòu)變化幅度，矩陣-2-范數(shù)越大，數(shù)據(jù)結(jié)構(gòu)的變化越大，反之則小。

本研究模擬試驗采用Linux shell進行數(shù)據(jù)整理，填補過程采用R語言編程實現(xiàn)，主要工具包為“RSNNS”和“missForest”。

結(jié) 果

1．估計偏差

圖3和圖4分別是WNT通路變量數(shù)據(jù)集和隨機變量數(shù)據(jù)集填補之后的估計偏差。由圖3和圖4可知，5種填補方法在填補精度上均隨著缺失數(shù)據(jù)比例的升高而降低，但MLP和kNN的穩(wěn)定性較高，在各種缺失比例情況下均高于其他填補方法。RF和MCMC方法估計偏差接近，二者均高于均值法。均值法隨著缺失比例的升高，對填補精度的損失比較大。

圖3 5種填補方法在WNT通路數(shù)據(jù)集中的估計偏差

圖4 5種填補方法在隨機變量數(shù)據(jù)集中的估計偏差

2．矩陣-2-范數(shù)

由圖5和圖6可知，在WNT通路變量數(shù)據(jù)集和隨機變量數(shù)據(jù)集中，MLP和kNN填補更加穩(wěn)健，均值填補則傾向于破壞原有數(shù)據(jù)結(jié)構(gòu)（變異被低估）。MCMC和RF填補在維持數(shù)據(jù)結(jié)構(gòu)方面亦優(yōu)于均值填補，但是不及MLP和kNN。

圖5 5種填補方法在WNT通路數(shù)據(jù)集中的矩陣-2-范數(shù)

圖6 5種填補方法在隨機變量數(shù)據(jù)集中的矩陣-2-范數(shù)

估計偏差和矩陣-2-范數(shù)從兩個方面評價填補效果，即填補精度和數(shù)據(jù)結(jié)構(gòu)穩(wěn)定性。然而，在實際高維數(shù)據(jù)挖掘工作中，對缺失數(shù)據(jù)的填補，還要考慮到填補效率的高低，見圖7。MCMC和RF填補在5種填補方法中耗時最長，且效果不及MLP和kNN。因此，不推薦應(yīng)用這兩種方法。得益于填補值計算的簡便，均值填補的填補時間最短。但均值填補的填補效果差強人意，亦不做推薦。MLP和kNN填補的填補效率高，以上兩種方法無論在估計偏差還是矩陣-2-范數(shù)方面，均表現(xiàn)優(yōu)異，同時，填補時間僅高于均值填補。

討 論

本研究中，根據(jù)評價指標和填補時間，讀者可以發(fā)現(xiàn)：kNN和MLP方法的填補效果要優(yōu)于均值法、MCMC以及RF。究其原因，kNN和MLP方法對于數(shù)據(jù)的分布類型并不敏感，穩(wěn)定性較好［8-9］，因此保證了在數(shù)據(jù)為多平臺來源時，填補效果較為可靠?；趯嶋H應(yīng)用考慮，雖然RF和MCMC的填補效果優(yōu)于均值法，但由于填補時間過長，不推薦RF和MCMC方法用于填補。均值法在數(shù)據(jù)缺失比例很低的情況下（≤5%）亦可以起到比較好的填補效果，同樣值得應(yīng)用。

在RF填補中，RF的每一棵樹代表一個獨立的非線性模型，多個模型組成的隨機森林在抗過擬合方面較單個模型更加可靠，但在一些噪音較大的回歸問題上，RF也會陷入過擬合，Weiss曾深入討論過此問題［10］。MCMC過程依賴于反復的模擬，形成足夠長的馬爾科夫鏈，對于高維的非線性數(shù)據(jù)，其效果有待商榷［11］。筆者推薦的兩種機器學習方法對于多平臺生物數(shù)據(jù)，有較好的穩(wěn)健性和容錯性，能夠在不依賴數(shù)據(jù)分布類型的情況下映射高維非線性復雜數(shù)據(jù)。MLP算法在層間權(quán)重矩陣的設(shè)定方面采用共軛梯度法，該方法不但克服了最速下降法收斂慢的弊端，也規(guī)避了牛頓法大量的矩陣運算過程，同時又不需要任何外來參數(shù)，是處理非線性高維數(shù)據(jù)最有效的方法之一。kNN算法通過屬性空間距離的最相近，構(gòu)造目標變量的候選集合，該方法對異常值不敏感，計算時間也相對較短。MLP和kNN兩種方法不僅適用于數(shù)值擬合，在判別分類方面亦有穩(wěn)定可靠的表現(xiàn)。鑒于以上兩類方法分別在隱藏單元個數(shù)和k值方面不能自適應(yīng)，需要預先設(shè)定，亦有研究者在MLP中采用隱藏單元個數(shù)等于訓練集樣本數(shù)，在kNN中k值等于訓練集樣本數(shù)的平方根的方法［12］。

圖7 5種填補方法填補時間對比

根據(jù)研究結(jié)果，不難發(fā)現(xiàn)，無論何種方法，在缺失比例升高的情況下，填補效果必然下降。當缺失比例高于70%時，填補效果較差，因此本研究未作更高缺失比例情況下的討論?！皦K缺失”數(shù)據(jù)可分為兩種情況：①測有某平臺的變量數(shù)據(jù)的樣本中，部分樣本不含有另一平臺的數(shù)據(jù)；②測有某平臺的變量數(shù)據(jù)的樣本中，所有樣本完全不含有另一平臺的數(shù)據(jù)。本研究主要討論了第一種情況。在第二種情況中，幾乎難以從現(xiàn)有數(shù)據(jù)中挖掘信息用于填補，此種情況下，只能從其他數(shù)據(jù)庫獲取先驗信息，當獲得信息足夠充分時，填補效果才能可靠，此情況有待進一步研究。

多組學數(shù)據(jù)的“塊缺失”也是完全隨機缺失的一種形式，對于其他缺失機制，并未作討論。矩陣-2-范數(shù)作為評價填補后數(shù)據(jù)集與原數(shù)據(jù)集中數(shù)據(jù)結(jié)構(gòu)變化的指標，考慮了矩陣之間的差異度，但是矩陣-2-范數(shù)僅用到最大特征根，也是本指標的一個局限，有待后續(xù)繼續(xù)研究。

［1］邱浪波，王廣云，王正志，等．基因表達缺失值的加權(quán)回歸估計算法．國防科技大學學報，2007，29（1）：111-115．

［2］張橋，李寧，張秋菊，等．任意缺失模式缺失數(shù)據(jù)不同填補方法效果比較．中國衛(wèi)生統(tǒng)計，2013，30（5）：690-692．

［3］吳俊杰，趙鵬．非線性噪聲數(shù)據(jù)集上基于隨機森林的空缺值填補算法．計算機應(yīng)用與軟件，2013，30（7）：51-53．

［4］W illiam S，Chad E，Herman T．Machine learning data imputation and classification in amulticohorthypertension clinical study．BioinformBiol Insights，2016，9（3）：43-54．

［5］Han D，Cao C，Su Y，etal．Ginkgo biloba exocarp extracts inhibits angiogenesis and its effects onWnt／beta-catenin-VEGF signaling pathway in Lewis lung cancer．J Ethnopharmacol，2016，192（1）：406-412．

［6］M ikhchi A，Honarvar M，Kashan NE，et al．Assessing and comparison of different machine learning methods in parent-offspring trios for genotype imputation．JTheorBiol，2016，399（2）：148-158．

［7］Beretta L，Santaniello A．Nearest neighbor imputation algorithms：a critical evaluation．BMC Med Inform DecisMak，2016，16（3）：63-74．

［8］Jerez JM，Molina I，García-Laencina PJ，et al．M issing data imputation using statistical and machine learning methods in a real breast cancer problem．ArtifIntell Med，2015，50（2）：105-115．

［9］Bibault JE，Giraud P，Burgun A．Big Data andmachine learning in radiation oncology：State of the art and future prospects．Cancer Lett，2016，382（1）：110-117．

［10］Weiss GM．M ining with rarity：A unifying framework．JSIGKDD Explorations，2004，6（1）：7-19．

［11］RiviereMK，Ueckert S，Mentre F．An MCMC method for the evaluation of the Fisher informationmatrix for non-linearmixed effectmodels．Biostatistics，2016，17（4）：737-750．

［12］BrettL．機器學習與R語言．李洪成，許金煒，李艦譯．機械工業(yè)出版社，2015：45-50．

（責任編輯：郭海強）

Evaluations on Several Im putation Approaches of Integrated Omics Data

Dong Xuesi，Lin Lijuan，Zhao Yang，et al（Department of Biostatistics，School of Public Health，Nanjing Medical University（211166），Nanjing）

Objective In post-GWAS era，integrated data from various platforms has become increasingly popular．Because of the complexity of data sources，many new challenges arise，which inevitably include how to treat“block missing data”．Ensuring the imputation accuracy and precision as well asmaintain the variance-covariance structure of the original data is of great importance to missing data imputation．In this project，we aimed to evaluate the effect of several imputationmethods based on both statistical techniques and machine learning techniques，on the integrated data from different data-platforms．Methods We go tlung cancer data-set（DNA methylation and gene expression）from The Cancer Genome Atlas（TCGA），and constructed m issing data-setw ith differentm issing proportions at5%，20%，35%，50%and 65%．The statisticalmethods（Mean imputation method，MCMC）and machine learningmethods（kNN，MLP，RF）were applied．Evaluation indicators included estimation bias and matrix 2-norms．At last，we considered imputation time and finding out a time-saving and efficientmethod．Results MLP and kNN showed high quality imputation effectand less time consuming from differentmissing ratio．Mean imputation had shortest filling time，and the imputation quality was high whenm issing ratio was low（≤5%）．However，whenmissing ratio increasing，the imputation effect decreased．When them issing ratio increasing，RF and MCMCmethod exceled in Mean approach．Nevertheless，RF and MCMC were time-killer．Conclusion After comprehensive comparative analysis，MLP and kNN imputation from machine learningmethods turned out to be suitable approaches in joint imputation process（DNA methylation，gene expression）．

Integrated omics data；Block m issing data；Statistical imputation；Machine learning imputation；Evaluation

本課題受國家自然科學基金重點項目（81530088）、面上項目（81473070，81373102）、國家自然科學青年基金（81402764）以及江蘇省高校優(yōu)勢學科資助

△通信作者：陳峰，E-mail：fengchen＠njmu．edu．cn