數(shù)據(jù)缺失情況下函數(shù)型數(shù)據(jù)聚類方法及應(yīng)用＊

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

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

1 引言

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

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

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

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

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

2 傳統(tǒng)函數(shù)型聚類方法的局限

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

2.1 原始數(shù)據(jù)聚類的局限性

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

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

2.2.1 篩選方法的原理

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

那么xi(t)即可表示為

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

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

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

如果基函數(shù)表達(dá)的連續(xù)時間函數(shù)形式為，

則可以有，使得

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

2.3.2 篩選方法的局限性

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

3 存在缺失的函數(shù)型數(shù)據(jù)自適應(yīng)聚類

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

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

3.1 自適應(yīng)聚類方法的優(yōu)勢

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

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

3.2 單變量函數(shù)型數(shù)據(jù)的聚類

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

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

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

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

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

在自適應(yīng)模型中，假設(shè)?i和γi都服從均值為0的正態(tài)分布。

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

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

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

其中，

另外，通過參數(shù)的估計值可以得出曲線的擬合結(jié)果。借助第k類的信息，可以計算出完整時間區(qū)間內(nèi)基函數(shù)的系數(shù)η^Mi，將基函數(shù)系數(shù)的估計值η^Mi

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

其中，

3.3 多變量函數(shù)型數(shù)據(jù)的聚類

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

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

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

4 模擬設(shè)定

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

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

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

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

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

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

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

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

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

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

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

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

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

5 應(yīng)用

5.1 數(shù)據(jù)形式

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

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

5.2 聚類分析

5.2.1 基函數(shù)個數(shù)的選擇

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

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

圖2 基函數(shù)個數(shù)的交叉驗證結(jié)果

圖3 聚類個數(shù)的選擇碎石圖

5.2.2 聚類個數(shù)的選擇

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

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

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

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

5.2.3 聚類結(jié)果的解釋

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

圖6-1 三類個體在心率變量上的表現(xiàn)

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

考慮到聚類結(jié)果的解釋意義，可以將老年人的宗氣水平劃分為三類。聚成三類的結(jié)果如下所示。圖5橫軸為個體降維后的α值，縱軸代表樣本編號。中間兩條豎線代表兩類的均值。當(dāng)類內(nèi)樣本的α值較為集中，類間樣本的α值距離較遠(yuǎn)時，即說明樣本的聚類效果較好。圖5中圓圈代表第一類樣本在低維空間的投影，三角形代表第二類樣本在低維空間的投影，十字代表第三類樣本在低維空間中的投影?？梢钥闯龅谝活悩颖舅鶎?yīng)的α值在3左右，第二類樣本所對應(yīng)的α值在8左右，第三類樣本所對應(yīng)的α值在-12左右，第三類和另外兩類的α值分隔較為清晰，類中心相距較遠(yuǎn)，第一類和第二類之間也有一定的區(qū)別。因此可以認(rèn)為統(tǒng)計上的聚類效果較好。

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

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

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

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

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

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

6 結(jié)論與討論

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

圖6-3 三類個體在體溫變量上的表現(xiàn)

圖6-4 三類個體在呼吸頻率變量上的表現(xiàn)

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

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

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

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 溫武兵.論宗氣的生理功能.山東中醫(yī)藥大學(xué)學(xué)報,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.