多組學聯合缺失數據填補方法的評價*

南京醫科大學公共衛生學院生物統計學系（211166） 董學思 林麗娟 趙 楊 魏永越 戴俊程 陳 峰

多組學聯合缺失數據填補方法的評價*

南京醫科大學公共衛生學院生物統計學系（211166） 董學思 林麗娟 趙 楊 魏永越 戴俊程 陳 峰△

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

多組學數據 塊缺失 統計學填補 機器學習填補 效果評價

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

傳統上，研究者在探索缺失值填補方面更側重于統計填補法，該類方法雖簡單易行，卻難以挖掘數據之間的深層關系，填補效果并不理想［2］；國外近幾年在機器學習領域的填補方法研究較多，該類方法對數據分布類型不敏感，通過計算機模擬的形式深入挖掘數據結構關系，填補效果較優［4］。現階段，數據填補的主要手段均基于以上兩大類方法（統計填補、機器學習填補），國內研究者進行過一系列統計填補方法的效果評價，但對統計填補方法和機器學習填補方法的比較并不多見。

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

數據與方法

1．數據

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

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

圖1 模擬試驗流程

2．方法

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

（1）統計填補

①均值法（mean）

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

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

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

a將數據集拆分為完整觀測部分Xfull和不完整觀測部分Xmiss。

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

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

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

d經過填補后，產生完整數據集，循環上述步驟，估算新產生的均向量和協方差矩陣進行填補，直至收斂。

（2）機器學習填補法

①隨機森林填補（RF）

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

a將完整數據部分作為自變量，缺失數據部分作為預測變量；

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

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

②K鄰近填補（kNN）

kNN填補在樣本數據集的特征空間中，按照馬氏距離選取相近（即特征空間中鄰近）的樣本集，計算對應變量的加權平均值進行填補。相較均值填補，kNN考慮了樣本間的變異，保持了數據結構的穩健性［7］。其步驟可以整理為：

a構建完整樣本集的矩陣結構；

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

d所計算的均值代替缺失數據。

③多層感知機填補（MLP）

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

圖2 多層感知機網絡結構

a構建完整數據集矩陣Xfull，作為訓練數據集，含有缺失變量的樣本作為預測集Xmiss；

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

c采用共軛梯度法，計算層間權重向量；

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

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

e隱藏單元再經過一次激勵函數轉化至輸出層。yk是輸出單元k＝1，…，k，wkh是第二層權重矩陣，wko是殘差項，g（）為線性激勵函數：

f產生神經網絡結構；

g輸入預測集樣本Xmiss，估計預測集中缺失數據。

5種方法的填補效果用估計偏差和矩陣-2-范數來評價。估計偏差：原數據集中變量均數，與填補后數據集中變量的均數之差的絕對值之和，再取平均值。可定義為：

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

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

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

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

結 果

1．估計偏差

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

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

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

2．矩陣-2-范數

由圖5和圖6可知，在WNT通路變量數據集和隨機變量數據集中，MLP和kNN填補更加穩健，均值填補則傾向于破壞原有數據結構（變異被低估）。MCMC和RF填補在維持數據結構方面亦優于均值填補，但是不及MLP和kNN。

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

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

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

討 論

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

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

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

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

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

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

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

［3］吳俊杰，趙鵬．非線性噪聲數據集上基于隨機森林的空缺值填補算法．計算機應用與軟件，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語言．李洪成，許金煒，李艦譯．機械工業出版社，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）以及江蘇省高校優勢學科資助

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