數據缺失情況下函數型數據聚類方法及應用＊

高明慧，易丹輝，彭 錦，胡鏡清，楊 燕

（1.中國人民大學應用統計科學研究中心 北京 100872；2.中國人民大學統計學院 北京100872；3.中國中醫科學院中醫臨床基礎醫學研究所 北京 100700;4.中國中醫科學院中醫基礎理論研究所北京 100700;5.湖北中醫藥大學 湖北 430065)

1 引言

隨著現代科學的發展，越來越多的時間序列數據對于分析現實問題及預測有重要意義。為充分利用時間觀測整體的信息，將時間觀測作為連續函數關系進行深入分析，James Ramsay和Bernard Silverman在1982年提出函數型數據的分析方法[1]（Ramsay，1982）。函數型數據分析認為離散的時間序列觀測由連續的時間上的函數產生，即離散觀測背后存在著關于時間的連續函數。函數型數據的優勢在于，一方面基函數擬合的思想提供了一種非參數的曲線擬合方法，不再拘泥于傳統已知函數形式的參數方法；另一方面，將離散的觀測連續化，通過求導運算可以更簡便的分析個體在時間上的變化趨勢[2]（Tokushige，2007）。

函數型數據分析的基本思想認為，時間序列觀測背后真實存在連續時間上的函數影響觀測的在不同時間的變化，因而尋找這一真實存在的函數關系是進行深入分析的關鍵。由于時間的連續性，函數曲線本質上是無窮維的，在一段時間內有無窮多個時間點，因此函數型數據在這一時間段內有無窮多個取值。但無窮維的數據無法使用傳統的多元統計方法進行分析，因此將無窮維的連續曲線降維到有限維的空間非常重要。

原始的函數型聚類方法通過觀測數據直接聚類[3]（Abraham，2003）。將每一個時間點看作一個變量，將離散時間點記錄的數據看作多變量數據進行多元聚類分析。Shuichi Tokushige教授在2007年提出用改進的kmeans聚類方法進行函數型數據的聚類分析，即對每一個時間點上的樣本觀測進行kmeans聚類。樣本在不同時間點所屬的類別不同，反映了每一類樣本不同的特點。這兩種方法并不適用于含有缺失的函數型數據聚類分析。Catherine Sugar教授提出基于模型的自適應聚類方法，通過極大似然的方法估計參數，利用所有樣本的觀測值信息進行聚類，解決了缺失數據難以直接進行基函數擬合的問題。

事實上，自適應方法在含有缺失數據的函數型聚類問題中應用廣泛。科羅拉多州立大學教授在2005年提出稀疏的函數型數據的自適應分析方法，其中提出利用函數型主成分分析法，在數據存在缺失情況下的函數曲線估計[4]。繼而有學者在2008年提出針對稀疏數據的聯合模型聚類方法，將函數曲線分解為固定效應和隨機效應兩部分，利用對數似然函數進行參數估計，從而實現函數型數據的聚類[5]。Peter Hall教授在2008年提出響應變量不服從連續的正態分布時的函數型數據分析方法，運用函數型主成分分析法對數據進行處理后，再根據主成分得分進行聚類分析[5]。其后，黃輝教授于2014年提出二分類變量函數型數據的聯合模型聚類方法，并對于不同類別的方差進行了差異化設計，提高了聚類結果的準確性[6]。另外，李浩成教授在2014和2015年提出函數型數據的分層模型分析方法，將函數型數據分成樣本個體、每天、每周的層次進行固定效應和隨機效應的分析[7]。常微分方程[8]和生存分析[3]也可結合到函數型數據的分析當中。

本文第二部分介紹傳統函數型數據聚類方法的局限性，第三部分介紹存在數據缺失情況下的自適應聚類方法，包括單變量時間觀測和多變量時間觀測兩種情況。第四部分通過模擬比較處理缺失數據的函數型聚類方法，第五部分通過中醫科學院提供的老年人宗氣數據進行自適應聚類方法的實證分析。第六部分回顧了方法的介紹和實證結果，并對未來的研究提出展望。

2 傳統函數型聚類方法的局限

傳統函數型聚類方法包括直接利用原始數據聚類和篩選方法（filtering method）聚類。這兩種方法在解決含有缺失的函數型數據聚類問題中均存在局限。

2.1 原始數據聚類的局限性

為解決時間序列觀測的聚類問題，最樸素的聚類方法將每一時間觀測點視為一個變量，離散時間點記錄的數據看作多變量的函數型數據進行多元聚類分析。這種方法在大部分情境中并不適用。一方面，當時間點的取值過于密集或測量時間段較長時，時間的維度數高于樣本數，聚類結果會受到維度過高的影響，尤其當變量數多于一個時，聚類分析的維度等于時間維度和變量維度的乘積，大大增加普通聚類分析的難度；另一方面，當樣本之間的時間點取值不同時，這種聚類方法難以將觀測時間點不同的樣本統一看作多元數據，將離散的時間點轉化為多維度的變量進行分析。另外，直接利用原始數據聚類忽視了時間上的連續性，不考慮時間先后的關系，會遺漏數據隱含的信息，導致分析結果不準確。在數據存在缺失的情況下，將每一個時間觀測點當作一個變量進行多元聚類，會由于每個樣本在多個變量上沒有觀測，無法進行聚類分析。

2.2 篩選法聚類（filtering method）的局限性

2.2.1 篩選方法的原理

Gareth James教授和Catherine Sugar教授首次將基函數擬合進行聚類的方法命名為篩選方法（filter?ing method）。篩選方法先用基函數對時間序列觀測進行擬合和降維，然后再對有限維的函數型數據進行聚類。篩選方法的聚類依據是降維后系數之間的距離。篩選方法第一步利用基函數進行降維。基函數的形式記為

那么xi(t)即可表示為

其中，基函數的選擇可以是傅里葉基函數、樣條基函數、小波基函數等。Ramsay的研究認為，當觀測在時間上呈現明顯的周期性時，可以選用傅里葉基函數；當觀測沒有周期性時，可以選用樣條基函數；當觀測存在頻繁且劇烈的波動時，建議選用小波基函數[1]（Ram?say，2006）。

圖1 某樣本觀測缺失時心律擬合值

當函數曲線表示為xi()t=CiΦi(t)時，為使得基函數能準確表示原函數，Ci的維度實際上是無窮維的，為了用更低的維度代替原始數據進行聚類分析，需要進行降維。

如果基函數表達的連續時間函數形式為，

則可以有，使得

K即利用前K個基函數表示的時間函數x~i()t可以近似真實的時間函數xi(t)。這樣通過犧牲一部分的準確性，完成了函數曲線從無窮維降到K維。由于樣本的基函數選擇相同，因此函數曲線可以由基函數系數代表，函數曲線之間的距離轉化為基函數系數之間的距離，作為聚類的依據。

2.3.2 篩選方法的局限性

篩選方法的關鍵在于用有限個基函數系數表示函數曲線。在分析含有缺失的函數型數據時，每個樣本由于數據缺失，觀測數較少，只能通過較少的基函數擬合函數曲線，而過少的基函數會導致擬合效果不理想，與真實曲線相差較遠的情況。如在宗氣數據中，心率變量的取值一般在60至80之間。由于數據缺失較為嚴重，某樣本在直接通過基函數擬合時，缺失部分的函數曲線的取值出現200以上或0以下的情況，如圖1所示，與人的心率變化完全不符。從圖1還可以看出，該樣本在下午13點至18點的心率變量使用篩選方法得到的心率擬合值在17點之后迅速降低至-1 500以下，原因在于17點之后數據沒有觀測，擬合時缺乏足夠的信息，導致基函數的擬合存在嚴重誤差。當基函數系數無法代表函數曲線時，根據基函數系數進行聚類的結果是不準確的。

3 存在缺失的函數型數據自適應聚類

自適應方法是篩選方法的改進，篩選方法中基函數系數是固定的，在確定函數曲線擬合的基函數形式后，基函數的系數隨之確定；而自適應方法中基函數系數是由類別中心決定的，并且存在隨機效應。類中心由該類別中所有樣本決定，因而，某個樣本的系數不僅僅由這個樣本自身的觀測決定，還受屬于同一類別的其他樣本的影響，故其系數可根據樣本所屬類別不同而變化。一方面保證基函數系數的靈活性，另一方面借助類內不同樣本的信息，使用最大似然法確定樣本所屬類別，在數據存在缺失的情況下可以有效補充缺失部分的信息。

自適應方法不僅適用完整時間觀測的函數型數據聚類，還用于含有隨機缺失的函數型數據聚類[9]（James，2003）。由于自適應方法同時確定樣本的系數和樣本所屬類別，因此需要迭代算法進行計算，本文在求解過程中選用EM算法。

3.1 自適應聚類方法的優勢

利用原始數據進行聚類，在處理含有缺失的函數型數據時，一方面會遇到時間觀測不完整，多數時間點觀測有缺失的問題，無法將n個時間點看作n個多元變量進行多元聚類；另一方面，當時間點過多時，多元聚類的變量維度過高，運算量非常大。篩選方法在對含有缺失的函數型數據進行聚類時，若每個樣本的觀測數較少，會使選用的基函數個數較少，難以對函數曲線進行較為準確的擬合，從而導致基函數系數之間的距離無法很好代表函數曲線之間的距離，致使聚類結果不準確。

自適應方法需事先確定基函數種類、個數和聚類個數。這種方法不直接使用基函數系數進行聚類，而是在已知分布信息的情況下，給定類中心系數等參數的初始值，運用EM算法，得到類中心對應系數的估計值；再根據函數型曲線的系數與類中心之間的馬氏距離進行聚類。

3.2 單變量函數型數據的聚類

令gi(t)是第i個個體在時間t上的真實值，Yi(t)是第i個個體在時間t上的觀測值，?i是第i個個體的測量誤差。真實值等于觀測值與測量誤差的疊加。

其中，s(t)是p維樣條基函數向量，ηi是服從正態分布的基函數系數向量。由于基函數的系數的隨機效應，ηi可以表示成與類別相關的固定效應μzi與隨機效應γi之和。（3）式是真實函數曲線由無限維向p維空間的投影，即第一次投影。

其中，zi表示個體i所屬類別。James和Sugar認為樣條基函數的系數αi根據所屬類別k不同，服從均值為μk，方差為Σ的混合高斯分布。對于第k類的類中心μk，

其中λ0和αk分別為p維和h維向量，Λ是 p×h維矩陣，其中h≤min(p,G-1)。（5）式是p維的類中心 μk向h維空間中αk的投影，即第二次投影。在保證不損失信息的前提下進一步降維，一方面減少了待估參數，另一方面使高維情況下的聚類結果可視化。結合（3）、（4）、（5）、（6）式可以得出函數型聚類的模型。

其中，n是個體數量。E()?i=0，測量誤差之間、測量誤差與gi之前均為獨立。真實的函數可以用基函數進行擬合。

在自適應模型中，假設?i和γi都服從均值為0的正態分布。

為保證模型可以求得唯一解，需要對λ0,Λ,αk加上兩個約束條件。

根據參數的分布假定可以得出，參數的極大似然函數可以寫成

通過EM算法求解（7）式極大似然的最大值可以估計出樣本所屬類別πk，殘差項方差Γ，隨機擾動的方差σ2，截矩項λ0，投影矩陣Λ，投影后的類中心系數αi。α^i與αk之間的馬氏距離決定了曲線i屬于第k類的概率，即可以通過（8）式得出聚類結果。

其中，

另外，通過參數的估計值可以得出曲線的擬合結果。借助第k類的信息，可以計算出完整時間區間內基函數的系數η^Mi，將基函數系數的估計值η^Mi

與基函數Si相乘，即得到曲線的擬合值。

其中，

3.3 多變量函數型數據的聚類

自適應方法同樣可以處理存在缺失的多變量函數型數據。Yij代表時間點tij1,…,tijnij上，第i個個體在第j個變量上的觀測值，J為變量個數。Yij=(Yijtij1,Yijtij2,…,

與上述過程相似，模型的形式為，

同樣可以用EM算法估計參數，判斷個體所屬類別，畫出擬合曲線的方法與單變量分析相似。

4 模擬設定

本文通過模擬比較不同缺失比例下函數型數據聚類方法的可行性，說明自適應方法在處理缺失數據時的優勢。可行性是指在多次重復的實驗中，該方法可以進行聚類分析的次數占總實驗次數的比例。

在重復多次的情況下，由于每次隨機生成的數據不同，缺失位置不同，聚類方法在相同缺失比例情況下，運行結果并不相同。模擬生成來自3個類別的100個樣本，每個樣本有60個時間觀測。由（3）式，樣本觀測等于誤差?加上基函數s(t)與系數η的乘積。即其中對于所有樣本，s(t)取基函數個數為10的B樣條基函數。第一類有20個樣本，η1為(η1,1,η1,2,…,η1,10)，10個元素取自均值為-10，方差為5的正態分布。第二類有30個樣本，η2η1為(η2,1,η2,2,…,η2,10)(η1,1,η1,2,…,η1,10)，10個元素取自均值為0，方差為5的正態分布。第三類有50個 樣 本 ，η3η2η1為 ( )η3,1,η3,2,…,η3,10(η2,1,η2,2,…,η2,10)(η1,1,η1,2,…,η1,10)，10個元素取自均值為10，方差為5的正態分布。誤差?服從標準正態分布。每個觀測時間點上隨機缺失的比例分別設定幾種情況：1%，2%，5%，10%，20%。模擬設定實驗重復1000次，記錄成功進行聚類的次數，并計算該次數與總重復次數的比例。

Yi=gi+ ∈i,i=1,…,n,?i～N(0 ,R),R=σ2I，gi(t)=s(t)Tηi。

表1 5%缺失下保留m個樣本的概率

表2 不同缺失比例下聚類方法的可行性結果

對于刪除法，即“直接丟棄含缺失數據的記錄”[23]的方法，當剩余完整觀測的樣本數量大于等于30時，認為處理后的數據集足夠代表原數據集的特征，可以進行聚類分析。對于篩選法，數據需要保證最小二乘估計的矩陣滿秩，估計出基函數系數后進行函數型聚類。對于自適應法，只要EM算法可以進行最大似然估計，就可以進行函數型聚類分析。

結合處理缺失數據的刪除法和不同聚類方法，模擬設定以下4種情況比較：刪除法處理缺失數據，得到完整數據后通過原始數據法進行函數型聚類（下記“刪除-原始數據法”）；刪除法處理缺失數據，得到完整數據后通過篩選法進行函數型聚類（下記“刪除-篩選法”）；直接使用篩選法進行函數型聚類（下記“篩選法”）；使用自適應法進行函數型聚類（下記“自適應法”）。其中原始數據法和篩選法由于自身的局限性，需要通過刪除法得到完整的數據以進行聚類分析；而自適應法本身可以處理含有缺失數據的函數型聚類，不需要刪除數據，因此不進行刪除法與自適應法結合的聚類。

一般在缺失比例較低時可以使用刪除法，但是當缺失比例較高，即模擬設定中缺失比例為5%及以上時，被刪除的樣本可能比較多，剩余樣本無法代表原數據集進行后續分析。表1以缺失比例5%的情況為例說明保留m個樣本的概率。令缺失比例為c，未被刪除的樣本個數為m，至少保留m個樣本的概率為

剩余樣本個數在[k,l)之間的概率為Pk,l

當m等于30時，至少保留30個樣本的概率為9.69E-17，十分接近0，說明經過刪除法處理后的樣本個數在30個以上的概率非常低。當m等于6時，保留的樣本個數小于6的概率為0.69，說明在缺失比例為5%的情況下，處理后的樣本個數很可能剩余不到6個。在1 000次的模擬中，每次剩余樣本的個數均不到30個，因此刪除法不可行。

不同缺失比例下聚類方法的可行性結果如表2所示。

從表2可以看出，原始數據法、篩選法在缺失比例為1%時，數據經過刪除法處理剩余較多樣本，可以進行聚類分析；在缺失比例小于等于2%時，方法可運行的比例在50%以上；當缺失比例達到5%時，樣本被刪除的概率很大，“刪除-原始數據法”和“刪除-篩選法”不再可行。缺失比例大于5%，三種方法基本失效。

因此，前三種處理含有缺失數據的函數型數據聚類方法在缺失比例達到10%及以上時均不可行，只有自適應法可以進行缺失情況下的函數型聚類分析。在實例分析中，老年人的宗氣數據在下午13點至18點之間的缺失比例高達70%，因此只能選擇自適應法對該數據進行函數型聚類分析。

5 應用

5.1 數據形式

數據來源于中國中醫科學院提供的老年人宗氣數據。宗氣的概念來自中醫，主要指人體內臟的消化運輸、升清降濁的功能。宗氣主宰著機體心肺樞機的功能活動，是心肺功能交互為用的結果和產物（溫武兵，2000）[10]。人體的心率、動脈血氧飽和度、體溫、呼吸頻率等特征可以反映宗氣足或不足。通過分析這些數據可以對老年人宗氣狀態進行分類。

數據形式為儀器檢測的58位老年人一天24小時內的心率、動脈血氧飽和度、體溫和呼吸頻率。經過數據清洗，截取下午一點至六點的時間序列觀測，每五分鐘取值一次，每個個體有60個觀測值。由于儀器本身的準確性以及老年人使用操作的影響，部分數據存在缺失，刪除缺失比例95%以上的個體，剩余52位老年人的觀測納入分析。根據心率、動脈血氧飽和度、體溫、呼吸頻率四個變量對人群聚類，在含有缺失數據的情況下使用多變量自適應方法進行函數型數據聚類。

5.2 聚類分析

5.2.1 基函數個數的選擇

在Ramsay的2006年的函數型分析書中提到，函數曲線不存在周期性時一般選用B樣條基函數進行擬合，而基函數個數的選擇可以借鑒BIC和CV等方法，也可以主觀選擇，基函數個數只要在合理的區間內，對于分析結果的影響并不明顯。周教授11在成對的稀疏函數型數據聯合模型分析一文中也提到，基函數主要起到平滑函數曲線的作用。在宗氣數據的研究中，本文將函數型數據進行了兩次投影，最終將無窮維的函數型數據降維至二維平面中，因此基函數主要起到了將離散的時間數據轉化為函數曲線這一作用，降維后基函數個數對于聚類結果的影響較小。另一方面，由于宗氣數據的稀疏性，當所有樣本需要用相同個數的基函數進行擬合時，最多只可以選用10個基函數。

本文通過十折交叉驗證(ten-fold cross validation)的方法確定基函數的個數，結果如圖2所示，其中橫軸為基函數個數，縱軸為交叉驗證的殘差平方和。由于樣本觀測的稀疏性，52個個體中在下午一點至六點之間觀測數最少為10個。Ramsay在函數型數據分析中提到，為保證基函數二階導的連續性，基函數的個數最少取4個。另外，基函數個數不應大于個體的觀測數，因此本文選擇基函數個數的交叉驗證的范圍為4至10個基函數。評判基函數的擬合效果的標準為，殘差平方和（Sum of the Squared Errors,SSE）越小擬合效果越好。從圖2看，適宜選擇10個基函數進行擬合。

圖2 基函數個數的交叉驗證結果

圖3 聚類個數的選擇碎石圖

5.2.2 聚類個數的選擇

聚類個數的選擇一方面需要考慮分類的效果，使類內距離更小，類間距離更大。參考James和Sugar提出的原則，使用平均類內馬氏距離dK確定聚類個數。

根據（9）式類內馬氏距離的變化選擇聚類個數。圖3中縱軸為（8）式求得的類內距離之差jump，橫軸為聚類個數。Jump越大，說明類內距離的變化越大，分隔不同個體的效果越明顯。從圖3中可以看出，聚成三類時jump最高，說明統計學上聚成三類是合理的。另外，結合宗氣水平的中醫理論，人群可分為宗氣充足、宗氣水平一般、宗氣不足三種，因此聚成三類在實際應用中有現實意義。

圖4 某個體的觀測值與擬合值

圖5 聚成三類的樣本在低維空間中的α值

5.2.3 聚類結果的解釋

圖4中粗實線代表一個體的含有缺失數據的觀測值，細實線代表個體的擬合值，虛線代表三個類別的平均水平。從圖4中可以看出，個體的觀測值存在缺失因而觀測曲線有多處間隔，而通過自適應方法的處理，個體的擬合值借助了三類平均水平的信息，因此可以得到連續的平滑的擬合值。這是自適應方法的優勢之一。

圖6-1 三類個體在心率變量上的表現

圖6-2 三類個體在動脈血氧飽和度變量上的表現

考慮到聚類結果的解釋意義，可以將老年人的宗氣水平劃分為三類。聚成三類的結果如下所示。圖5橫軸為個體降維后的α值，縱軸代表樣本編號。中間兩條豎線代表兩類的均值。當類內樣本的α值較為集中，類間樣本的α值距離較遠時，即說明樣本的聚類效果較好。圖5中圓圈代表第一類樣本在低維空間的投影，三角形代表第二類樣本在低維空間的投影，十字代表第三類樣本在低維空間中的投影。可以看出第一類樣本所對應的α值在3左右，第二類樣本所對應的α值在8左右，第三類樣本所對應的α值在-12左右，第三類和另外兩類的α值分隔較為清晰，類中心相距較遠，第一類和第二類之間也有一定的區別。因此可以認為統計上的聚類效果較好。

三類人群在心率、動脈血氧飽和度、體溫、呼吸頻率四個變量上的表現不同。在心率方面，第三類人群的心率較快，且波動較大，在13點至14點之間出現了明顯的峰值；第一類人群和第二類人群的心率均較為平穩，其中第二類人群的心率相比于第一類人群更慢。

在動脈血氧飽和度方面，第一類人群的動脈血氧飽和度較為平穩，呈現緩慢上升的趨勢。第二類人群的動脈血氧飽和度稍有波動，第三類人群的波動最為明顯，其平均水平也低于前兩類人群。

在體溫方面，第一類和第二類人群的體溫較為相近，在13點到14點之間，和15到16點之間有小幅度上升，整體波動不大。第三類人群的體溫波動較為明顯，尤其在16點之后，出現明顯的下降和上升。總體上看，第三類人群的體溫低于前兩類人群。

在呼吸頻率變量上，第二類人群的呼吸頻率變化最為平緩，第一類人群的平均呼吸頻率略高于第二類人群且略有波動。第三類人群的呼吸頻率明顯高于前兩類人群，并且在14點、16點等時間點周圍出現大幅波動。

綜合四個變量的曲線圖中信息可以看出，在下午13點至18點之間，第一類人群和第二類人群的特征較為相似，心率較慢且較為平穩；動脈血氧飽和度變化幅度小，且呈現緩慢上升的趨勢；體溫呈現小幅波動；呼吸頻率變化平緩。其中第二類人群的心率和呼吸頻率低于第一類人群。第三類人群在心率、動脈血氧飽和度、體溫、呼吸頻率四個變量上的波動幅度相對較明顯，并且在心率和呼吸頻率變量上水平明顯高于第一類和第二類人群，在動脈血氧飽和度和體溫變量上的水平略低于前兩類人群。

中醫理論發現，宗氣不足的主要表現有心率快、血氧低、體溫低、呼吸頻率高，并且在每個變量上的波動幅度均比較大。因此可以認為，聚類所得三類人群中，第三類人群與第一類和第二類人群明顯不同，有心率和呼吸頻率較高、血氧含量和體溫較低、四個指標均波動較大的特點，屬于宗氣不足人群。第二類人群相比于第一類人群的心率和呼吸頻率更低，且波動幅度更加平緩，可以將第二類人群定義為宗氣充足人群。第一類人群在各變量的水平介于第二類和第三類人群之間，與第二類宗氣充足人群更為相近，因此可以將第一類人群定義為宗氣水平一般的人群。

6 結論與討論

本文主要介紹了函數型數據的特征，函數型數據的聚類方法以及中醫科學院來年人宗氣數據的實證分析。函數型數據分析的關鍵點在于降維，常用的降維方法是基函數法，用有限個基函數與系數的乘積表示連續的函數曲線。函數型聚類的方法有原始數據聚類法、篩選方法和自適應方法。原始數據聚類法當某些時間點上的樣本觀測存在較多缺失時，聚類方法無法計算樣本所屬類別，缺失比例增大，會導致無法進行聚類。篩選方法通過基函數法對函數曲線降維，再對基函數的系數進行多元聚類分析，當數據存在缺失時，擬合效果不理想，聚類效果不好，缺失比例增大也會出現無法聚類的情況。自適應方法在篩選方法的基礎上，假定基函數系數根據樣本所屬類別不同服從不同的分布，這種方法提高了擬合的靈活性，并且適合處理稀疏數據的聚類問題。自適應方法也有其自身的局限性，一方面，當數據缺失過多，某些時間點上的觀測值小于兩個甚至沒有觀測時，自適應方法有可能無法運行。另一方面，本文僅討論了自適應方法在缺失情況下進行聚類的可行性，其聚類效果有待進一步研究。本文討論數據缺失的處理時，僅考慮了刪除法，沒有涉及插補法，因為當數據缺失較多時，插補的結果很可能是不準確的，基于插補數據的后續分析也很可能出現問題。函數型數據如何進行插補還有待研究。

圖6-3 三類個體在體溫變量上的表現

圖6-4 三類個體在呼吸頻率變量上的表現

應用部分主要分析了老年人在下午一點至六點的時間序列觀測，通過心率、動脈血氧飽和度、體溫、呼吸頻率四個變量的聚類結果分析老年人宗氣狀況。利用針對含缺失數據的自適應聚類方法，將老年人分為宗氣充足、宗氣水平一般、宗氣不足三類人群，其中宗氣不足人群與前兩類人群有明顯差別。

本文中選取基函數個數的方法為交叉驗證，使殘差平方和最小的基函數個數為最優。由于擬合的函數曲線非線性，殘差平方和難以反映擬合的效果，因此在未來的研究中可以考慮其他反應擬合效果的統計量代替殘差平方和。另外，基函數個數的選擇與樣本觀測數之間的關系也有待探究。自適應方法中，擬合每一條樣本觀測的基函數個數受樣本觀測稀疏性的限制，而聚類后每一類內的總觀測數相對較充足，可以考慮擴大基函數個數的選擇范圍。

一般情況下，二次投影后的維度h選取1或2以保證可視化的效果。本文選擇將無限維的函數曲線二次投影到一維空間上，在未來的研究中可以討論將函數曲線投影到二維空間或更高維空間中的情況。

1 Ramsay J O.Functional data analysis.John Wiley&Sons,Inc.,2006.

2 Tokushige S,Yadohisa H,Inada K.Crisp and fuzzy k-means clustering algorithms for multivariate functional data.Comput Stat,2007,22(1):1-16.

3 Abraham C,Cornillon P A,Matzner-L?ber E,et al.Unsuper?vised curve clustering using B-splines.Scandi J stat,2003,30(3):581-595.

4 Yao F,Müller H G,Wang J L.Functional data analysis for sparse longitudinal data.J Ame Stat Associ,2005,100(470):577-590.

5 Hall P,Müller H G,Yao F.Modelling sparse generalized longi?tudinal observations with latent Gaussian processes.J Royal Stat So?ci:Series B(Statistical Methodology),2008,70(4):703-723.

6 Huang H,Li Y,Guan Y.Joint modeling and clustering paired generalized longitudinaltrajectories with application to cocaine abuse treatment data.J Ame Stat Associ,2014,109(508):1412-1424.

7 Li H,Kozey Keadle S,Staudenmayer J,et al.Methods to as?sess an exercise intervention trial based on 3-level functional da?ta.Biostatistics,2015,16(4):754-771.

8 Little R J A,Schenker N.Missing data[M]//Handbook of statisti?cal modeling for the social and behavioral sciences.Springer US,1995:39-75.

9 James G M,Sugar C A.Clustering for sparsely sampled func?tional data.J Ame Stat Associ,2003,98(462):397-408.

10 溫武兵.論宗氣的生理功能.山東中醫藥大學學報,2000,24(4):247-250.

11 Zhou L,Huang J Z,Carroll R J.Joint modelling of paired sparse functionaldata using principalcomponents.Biometrika,2008,95(3):601-619.