基于缺失森林模型的稀疏函數型數據修復方法

摘 要：為解決函數型數據缺失插補問題，在函數型數據分析框架下，以缺失森林模型（MF）為基礎，采用基于條件期望主成分分析的函數型插補方法PACE進行初始插補，并通過K-means聚類借助樣本之間的相關性插補，給出了一種融合類信息的函數型多重插補方法。模擬數據插補實驗結果表明，在不同缺失比例（5%～55%）下，本文方法相較于Hot.deck、MF、均值插補、PACE、MFP、SFI、HFI等7種插補方法，更能保證插補的準確性和有效性。同時，針對股票數據的實例應用驗證了本文方法插補得到的數據符合實際情況和規律。

關鍵詞：函數型數據；缺失森林；多重插補；缺失插補方法

中圖分類號：O212"" 文獻標志碼：A""" 文章編號：1673-5072（2024）05-0481-07

隨著數據采集程度的不斷提高，出現了一種無窮維、具備連續特征的數據，稱為函數型數據。它是一種基于時間或空間的數據類型，描述某個變量隨著時間或空間的推移而變化的規律，例如股票價格、氣溫和心率數據等都是函數型數據。如果將這些具備函數特征的數據看成一個整體進行統計分析，則稱為函數型數據分析 （Functional Data Analysis，FDA）［1］。FDA方法已經成為醫學［2-3］、經濟［4］、環境［5］等領域中處理大規模復雜數據的重要工具。

函數型數據往往不能被直接觀測，實際采集結果通常是具備曲線特征的離散采樣點，然而在函數型數據的收集過程中，數據缺失是較為常見的數據質量問題。例如：在環境監測和交通系統中，由于設備、電源、傳輸以及監測點增設或停運等原因，通常會導致監測數據存在大量的條狀、塊狀缺失［5-6］；在醫學領域中，由于個體縱向數據搜集、測量成本過高等原因，常會導致測量的數據存在缺失或稀疏不規則［7-8］；在投資理財領域，股票數據經常由于停盤等因素導致其中存在缺失［9］等。數據缺失不僅會妨礙現有FDA方法的應用，還會降低樣本的代表性，不同程度增大統計分析的復雜性和難度，降低統計推斷的精度。因此，如何科學有效地處理缺失數據是各領域亟待解決的問題。

近年來，學者們對函數型數據插補方法展開了研究。如Yao等［10］構建了一種基于條件期望主成分分析（Principal Components Analysis through Conditional Expectation，PACE）的稀疏不規則函數型數據插補方法，該方法是一種單一插補方法，其參數估計概念簡單，且易于實現。然而PACE仍存在一些不足之處：第一，PACE沒有考慮后續的模型是否適合，這導致估計模型參數時存在偏差［11］；第二，PACE作為一種單一插補方法，在形成置信區間、預測區間或p值時，沒有考慮插補的不確定性［12］；第三，該方法的預測精度高度依賴數據的分布，在處理非線性數據時存在一定的局限性。為了克服單一插補方法局限，有學者提出了一些針對函數型數據的多重插補方法。例如：對于縱向缺失數據，Twisk等［13］利用多重插補（Multiple Imputation，MI）算法估計缺失值；Ferro［14］提出了橫截面多重插補方法，并提供了完全信息最大似然的相似估計；此外，Adam等［8］將鏈式方程多重插補（Multiple Imputation by Chained Equations，MICE）與函數型回歸模型相結合，提出針對響應變量缺失的fregMICE算法；Rao和Reimherr［15］提出修復函數型缺失數據的多重插補方法，并應用于處理稀疏的血壓監測數據。然而這些函數型多重插補方法在處理缺失數據時并未充分考慮到樣本曲線之間的相關性和規律性。如空氣質量監測站點在時空上具有一定的相關性，距離越近的監測站點在同一時間段觀測到的數據之間相關性越強，所以樣本之間的相關性是插補缺失值的一項重要依據［16］。

因此，針對上述問題，本文提出一種基于缺失森林模型（Missforest，MF）的稀疏函數型數據修復方法（Missforest Combining Class Information and PACE，CMFP），該方法以MF為核心，采用PACE初始插補，結合K-means聚類挖掘樣本曲線之間的潛在變化模式。本文充分利用樣本曲線的函數特征和類信息，給出一種融合類信息的函數型多重插補方法，為稀疏函數型數據提供一種迭代插補方案。

1 基于缺失森林模型的稀疏函數型數據修復方法

考慮L2（τ）中的隨機函數X（t）用于描述樣本曲線的變動軌跡，其中L2（τ）表示封閉時間間隔τ上平方可積函數的Hilbert空間。假設X（t）有均值函數E（X（t）=μ（t））和協方差函數CX（t，s）=cov（X（s），X（t）），PACE［10］通過條件期望估計函數型主成份得分，其主成分展開式如下：

Xi（t）=μX（t）+∑SymboleB@j=1ξijvj（t），（1）

其中vj（t）是CX的特征函數。主成分得分通過式（2）得到，

ξij=〈Xi－μX，vj〉。（2）

MF具體執行步驟如算法1所示。

算法1：MF插補算法

輸入：n×p的函數型缺失數據矩陣X，終止迭代條件γ。

步驟1 對存在缺失值的變量運用該變量的平均值對缺失值進行插補；

步驟2 計算缺失數據集X中各個變量的缺失率，將缺失率從小到大排序，并將對應的變量存入向量m中；

步驟3 判斷是否達到γ和最大迭代次數，若達到則算法停止，輸出插補后的矩陣Ximpnew，記為Xnew，若沒達到，則繼續步驟4—6；

步驟4 存儲先前插補后的矩陣，記為Ximpnew；

步驟5 對于s∈m，依次執行：

①訓練隨機森林模型：y（s）obs～x（s）obs；

②利用①訓練好的隨機森林模型輸入x（s）mis預測y（s）mis值；

③使用y（s）mis值更新插補后的矩陣，記為Ximpnew；

步驟6 更新γ和最大迭代次數，返回步驟3。

輸出：插補后的矩陣Ximp。

其中Xs表示含有缺失的變量；y（s）obs表示Xs的觀測值；y（s）mis表示Xs的缺失值；x（s）obs表示Xs以外的觀測值；x（s）mis表示Xs的缺失值以外的其余觀測值。終止條件γ定義為Ximpnew和Ximpold之間的差異，如果γ增大，則達到終止條件，算法結束，輸出插補結果Ximp。 Ximpnew和Ximpold間的差異定義如下：

ΔN=∑j∈N（Ximpnew－Ximpold）2∑j∈N（Ximpnew）2，

其中，N為各變量的集合。

CMFP方法的執行步驟如算法2所示。

算法2：CMFP插補算法

輸入：n×p的函數型缺失數據矩陣X，終止迭代條件γ。

步驟1 對缺失數據矩陣X利用PACE進行初始插補；

步驟2 計算初始插補后數據的函數型主成分得分，使用K-means對函數型主成分得分進行聚類，將樣本劃分為k（k=1，2，…，N）類，記最終聚類結果為{C1，C2，…，CN}；

步驟3 對{C1，C2，…，CN}分別使用MF方法，得到對應的插補結果{C1，C2，…，CN}；

步驟4 將插補結果{C1，C2，…，CN}按照缺失數據矩陣X中變量的順序進行整合。

輸出：插補后的矩陣C。

2 模擬試驗

2.1 試驗設計

為驗證CMFP方法在插補函數型數據時的有效性，構造模擬數據集，并在不同缺失比例下評價其插補性能。本文通過以下3個步驟完成試驗設計：

步驟1 隨機構造一個100×50的模擬數據集［15］。具體為：模擬50條獨立同分布的曲線{X1（t），…，X50（t）}，這些曲線服從均值為0，協方差為CX（t，s）=σ2Γ（v）2v－12v|t－s|ρvKv2v|t－s|ρ的高斯分布，其中Kv是第二類的修正貝塞爾函數。設置參數ρ=0.5，v=5/2，σ2=1，曲線在區間［0，1］上等距選取100個時間點進行估計。假設每一個觀測點均包含一個均值為0，方差為3的正態測量誤差，最終得到函數型數據離散的觀測矩陣。

步驟2 隨機生成含有缺失的數據集。為驗證CMFP方法在不同缺失比例下均有較好的插補效果，設置缺失率分別為5%、15%、25%、35%、45%、55%。

步驟3 對比方法與評價指標的確定。將CMFP方法與Hot.deck［17］、MF［18］、均值插補［19］、PACE［10］、MFP［15］、SFI［20］、HFI［20］等7種方法進行插補性能的對比。評價指標采用平均絕對誤差（Mean Absolute Error，MAE）和均方根誤差（Root Mean Square Error，RMSE），其計算公式如下：

MAE=1n1∑n1i=1y^i－yi，RMSE=1n1∑n1i=1（y^i－yi）2，

其中，y^i是第i個缺失值的插補值，yi是對應的真實值，n1為不完整變量中缺失值的個數。MAE和RMSE值的大小反映了插補值與真實值之間的誤差，其值越小表示插補值與真實值越接近，插補效果越好，算法性能越高。

2.2 確定最佳聚類數目

在CMFP方法中，聚類數目會對插補結果產生一定的影響，將聚類數目設置為1～5，通過組內誤差平方和（Sum of Squared Error，SSE）確定最佳聚類數目，不同聚類數目下SSE值如圖1所示。隨著聚類數目增多，樣本數量越來越少，距離越來越近，因此SSE值隨著聚類數目增多而減少，當SSE值下降較大存在“肘點”或減少較緩慢時，插補誤差趨于平穩，停止聚類。從圖1可以看出，當聚類數目為2時存在“肘點”，且聚類數目大于等于3時，SSE值下降趨勢平緩。進一步通過RMSE和MAE 2個評價指標確定最佳聚類數目，評價指標結果如表1所示。在不同缺失比例下，CMFP方法聚類數目為2的插補誤差均小于聚類數目為3的誤差。因此，設定最佳聚類數目為2。

2.3 試驗結果

試驗主要通過R語言實現，計算機環境為：Intel（R） Core（TM） i5-5200U CPU2.20 GHz，內存4GB，Windows10 64位操作系統。試驗結果如表2、表3所示：任何缺失比例下，Hot.deck、均值插補、PACE、SFI和HFI的插補誤差均較大，其原因在于Hot.deck和均值插補作為傳統多元統計插補方法，在插補函數型數據時并沒有考慮到樣本的曲線特征，而PACE、SFI和HFI作為單一函數型插補方法，插補有效性較差；MF、MFP和CMFP 3種方法均具有較小的插補誤差，且較為接近。為了直觀地觀察MF、MFP和CMFP 3種方法的插補效果，將評價指標可視化，如圖2和圖3所示：在不同的缺失率下，CMFP方法相較于MF和MFP方法，RMSE分別降低了1.08%～8.53%和1.27%～8.78%，MAE分別降低了0.82%～4.91%和0.61%～4.37%，CMFP方法插補性能在8種方法中表現最優。

3 實例應用

為了進一步驗證CMFP方法的實際應用效果，將該方法應用于樣本曲線變化波動較大的股票交易收盤價數據。由于股票分類中最常用的方法是按照股票行業進行分類，并且曾有學者指出行業分類的結果會影響股票的價格［21］，同時行業因素能解釋股票收益率26%波動率的結論，說明同行業的股票數據之間具有一定的相關性，可按照行業對股票數據進行分類，因此按照2021年第3季度證監會上市公司行業分類結果（http：//www.csrc.gov.cn/csrc/c100103/c1558619/content.shtml），以2022年1月1日—2022年12月31日交易數據中的收盤價作為研究對象來驗證CMFP方法的適用性。現選取其中4類中的24只股票，包括農業類5只、開采輔助活動類4只、貨幣金融服務類8只、綜合類7只（表4）。

為更好地驗證CMFP方法的有效性，選取的股票數據不含缺失值，因此對股票數據進行5%、25%、55%的隨機缺失，并應用CMFP方法估計缺失數據。以農業類股票為例，基于皮爾遜相關系數將農業類5只股票完整數據（缺失率為0）的相關系數和在不同缺失比例下插補后數據的相關系數進行對比，如果完整數據與插補后數據的相關系數比較接近，則說明插補值符合原有數據規律，插補方法有效。從表5可看出，不論缺失比例如何，完整數據與插補后數據的相關關系基本保持一致，例如完整數據中隆平高科與荃銀高科的相關系數是0.8658，不同缺失比例下，利用CMFP方法插補后數據的相關系數分別為0.8666、0.8657和0.8684，插補值符合原有數據的潛在變化規律，故CMFP方法估計的缺失值有效。

此外，為了說明CMFP方法對后續統計分析的影響，對不同缺失率下插補的完整數據進行K-means聚類，結果如表6所示。

結合表4、表6可以看出，缺失率為5%和55%時，第一類股票與農業類股票相比，增加了宏大爆破和金誠信2只股票；第二類股票與開采輔助類活動類股票相比，有3只股票的變動，分別增加了紫金銀行，減少了宏大爆破和金誠信2只股票；第三類股票與貨幣金融服務類股票相比，有2只股票的變動，分別增加了創元科技，減少了紫金銀行；第四類股票與綜合類股票相比，減少了1只股票，為創元科技。缺失率為25%時，第一類股票與農業類股票相比，增加了1只股票，為中曼石油；第二類股票與開采輔助類活動類股票相比，減少了1只股票，為中曼石油；第三類股票與貨幣金融服務類股票相比以及第四類股票與綜合類股票相比，均無變動。

通過對比不同行業股票分類結果與不同缺失率下插補后聚類分析結果可以得到，插補后不同類別內股票數量及名稱沒有明顯變動，說明插補值對聚類分析的影響較小。因此，CMFP方法插補得到的數據符合實際情況，可以在一定程度上降低缺失值對后續統計分析和模型性能的影響。

4 結 論

近年來，隨著數據科學和機器學習技術的發展，對稀疏函數型數據的修復方法也成為了研究的熱點之一。本文在彌補已有插補方法不足的基礎上，提出一種基于缺失森林模型的稀疏函數型數據多重插補方法CMFP，旨在為稀疏函數型數據插補提供一種可行方案。通過在模擬數據和實例數據上實驗的結果表明，CMFP方法相較于其他插補方法具有較小的插補誤差，在一定程度上提高了插補的準確性。同時，通過相關系數矩陣表明，該方法插補得到數據與現實情況相符，可為后續統計分析奠定較為準確的數據基礎。

參考文獻：

［1］ RAMSAY J O，SILVERMAN B W.Functional data analysis［M］.2nd.New York：Springer，2005.

［2］ HAREZLAK J，WU M C，WANG M，et al.Biomarker discovery for arsenic exposure using functional data analysis and feature learning of mass spectrometry proteomic data［J］.Journal of Proteome Research，2008，7（1）：217-224.

［3］ GERTJEISS J，GOLDSMITH J，CRAINICEANU C，et al.Longitudinal scalar-on-functions regression with application to tractography data［J］.Biostatistics，2013，14（3）：447-461.

［4］ RAMSAY J O，RAMSAY J B.Functional data analysis of the dynamics of the monthly index of nondurable goods production［J］.Journal of Econometrics，2002，107（1）：327-344.

［5］ 黃恒君，漆威.海量半結構化數據采集、存儲及分析：基于實時空氣質量數據處理的實踐［J］.統計研究，2014，31（5）：10-16.

［6］ CHIOU J M，ZHANG Y C，CHEN W H.A functional data approach to missing value imputation and outlier detection for traffic flow data［J］.Transportmetrica B：Transport Dynamics，2014，2（2）：106-129.

［7］ JAMES G M，SUGAR C A.Clustering for sparsely sampled functional data［J］.Journal of the American Statistical Association，2003，98（462）：397-408.

［8］ ADAM C，EVA P，OFER H.Elucidating age and sex-dependent association between frontal eeg asymmetry and depression：An application of multiple imputation in functional regression［J］.Journal of the American Statistical Association，2022，117（537）：12-26.

［9］ 鐘宇航.基于廣義Rescal分解的股票缺失值填充［D］.成都：西南財經大學，2022.

［10］YAO F，MULLER H，WANG J.Functional data analysis for sparse longitudinal data［J］.Journal of the American Statistical Association，2005，100（470）：577-590.

［11］PETROVICH J，REIMHERR M，DAVMONT C.Highly irregular functional generalized linear regression with electronic health records［J］.Journal of the Royal Statistical Society Series C：Applied Statistics，2022，71（4）：806-833.

［12］龐新生.缺失數據插補處理方法的比較研究［J］.統計與決策，2012（24）：18-22.

［13］TWISK J，DE B M，DE V W，et al.Multiple imputation of missing values was not necessary before performing a longitudinal mixed-model analysis［J］.Journal of Clinical Epidemiology，2013，66（9）：1022-1028.

［14］FERRO M A.Missing data in longitudinal studies：cross-sectional multiple imputation provides similar estimates to full-information maximum likelihood［J］.Ann Epidemiol，2014，24（1）：75-77.

［15］RAO A R，REIMHERR M.Modern multiple imputation with functional data［J］.Stat，2021，10（1）：e331.

［16］ZHANG H S，ZHANG Y，LI Z H，et al.Spatial-temporal traffic data analysis based on global data management using MAS［J］.IEEE Transactions on Intelligent Transportation Systems，2004，5（4）：267-275.

［17］BERTSIMAS D，PAWOLOWSKI C，ZHOU Y D.From predictive methods to missing data imputation：an optimization approach［J］.The Journal of Machine Learning Research，2017，18（1）：7133-7171.

［18］STEKHOVEN D J，BUHLMANN P.Missforest--non-parametric missing value imputation for mixed-type data［J］.Bioinformatics，2012，28（1）：112-118.

［19］JINNINEN H，NISKA H，TUPPURAINEN K，et al.Methods for imputation of missing values in air quality data sets［J］.Atmospheric Environment，2004，38（18）：2895-2907.

［20］KIDZINSKI L，HASTIE T.Longitudinal data analysis using matrix completion［J］.ResearchGate，2018.

［21］張紫璇，段紅梅.基于混沌游戲表示和自適應仿射傳播聚類的股票板塊分類［J］.財會月刊，2020（19）：152-155.

Restoration Method of Sparse Functional Data Based on Missforest Model

GAO Hai-yanab，LI Wei-xina，MA Wen-juana

（a.School of Statistics and Data Science，b.Key Laboratory of Digital Economy and Social Computing Science，Lanzhou University of Finance and Economics Lanzhou Gansu 730020，China）

Abstract：In the framework of functional data analysis and with the aid of Missforest model （MF），this paper proposes a functional multiple interpolation method combining class information by employing the functional interpolation method PACE based on conditional expectation principal component analysis for initial interpolation，and utilizing the correlation interpolation between samples through K-means clustering in order to solve the missing interpolation problem of functional data.The experiment results of simulation data interpolation show that the proposed method under different missing rates （5%～55%） can ensure the accuracy and effectiveness of interolation in a better way，compared with other seven imputation methods such as Hot.deck，MF，Mean imputation，PACE，MFP，SFI and HFI.At the same time，an application example of stock data has verified that the data imputed by the proposed method conforms to the actual situation and rules.

Keywords：functional data; Missforest; multiple imputation;missing imputation method