確定函數(shù)型數(shù)據(jù)主成分個數(shù)的自適應(yīng)加權(quán)截斷法

田 密，羅幼喜

(湖北工業(yè)大學(xué)理學(xué)院，湖北 武漢，430068)

函數(shù)型數(shù)據(jù)[1-2]具有函數(shù)的特征，對其處理不再是把單個數(shù)據(jù)看作樣本，而是將數(shù)據(jù)擬合為一條函數(shù)曲線作為樣本進(jìn)行分析。由于函數(shù)型數(shù)據(jù)理論上隸屬于無窮維空間上的數(shù)據(jù)，因此實際運用中通常需要對其進(jìn)行基函數(shù)展開降維處理。目前常用的基函數(shù)有兩類：第一類是給定的基函數(shù)，如Zhou等[3]、Lian等[4]對函數(shù)型變量使用的B樣條基展開，但是這種基函數(shù)與數(shù)據(jù)是獨立的，不隨數(shù)據(jù)的改變而發(fā)生變化，因此難以確定原始數(shù)據(jù)的主要信息由多少個基函數(shù)決定；第二類是結(jié)合樣本數(shù)據(jù)形成的基函數(shù)，如Shin等[5]、Zhang等[6]使用的函數(shù)型主成分基展開，其中所選取的主成分基是由數(shù)據(jù)所驅(qū)動，它隨著數(shù)據(jù)的不同而相應(yīng)發(fā)生改變，并且能夠解釋函數(shù)型變量變差的百分比，因此這類主成分基展開更具有實際意義。然而，對函數(shù)型數(shù)據(jù)采用主成分基展開也涉及到基函數(shù)的個數(shù)K的確定這一重要問題。K的選取表現(xiàn)為估計的偏差與方差之間的平衡:K越大，估計的偏差越小，而方差越大；反之，K越小，估計的偏差越大，而方差越小。

眾多研究者從不同角度提出了確定函數(shù)型數(shù)據(jù)主成分個數(shù)的相關(guān)準(zhǔn)則。Rice等[7]提出依據(jù)方差解釋百分比(percentage of variance explained, PVE)對主成分?jǐn)?shù)量進(jìn)行截斷；Müller等[8]也采用此方法確定函數(shù)型主成分的個數(shù);但是Zhu等[9]研究指出，盡管基于PVE的主成分基截斷操作簡單，但這種方法僅考慮了特征值的大小，而對于許多復(fù)雜問題，函數(shù)型主成分的個數(shù)對響應(yīng)變量的影響并不僅僅與特征值大小有關(guān)。Su等[10]用PVE法確定主成分基個數(shù)并應(yīng)用于假設(shè)檢驗，證實了PVE并不是用于檢驗的最佳標(biāo)準(zhǔn)，于是提出關(guān)聯(lián)變異指數(shù)(association-variation index，AVI)這一概念，并由此建立關(guān)聯(lián)變異解釋百分比準(zhǔn)則(percentage of association-variation explained，PAVE)替代PVE用于函數(shù)型主成分的排序和選擇。雖然PAVE法能夠考慮到回歸模型中函數(shù)型變量對響應(yīng)變量的影響，但這種方法忽略了由于特征值迅速衰減導(dǎo)致的函數(shù)型協(xié)變量的主要變異。

鑒于以上方法的不足，本文提出一種新的自適應(yīng)加權(quán)截斷法用于函數(shù)型主成分個數(shù)的選擇。該方法首先分別基于PVE和PAVE對函數(shù)型數(shù)據(jù)的主成分基進(jìn)行排序和選擇，然后對兩個方法截斷出的基函數(shù)個數(shù)進(jìn)行加權(quán)，通過求解最小化估計誤差獲得最優(yōu)權(quán)重參數(shù)，從而最終確定主成分基展開項數(shù)。該方法不僅考慮了由于特征值迅速衰減導(dǎo)致的主要變異，還將函數(shù)型協(xié)變量和響應(yīng)變量之間的關(guān)聯(lián)納入選擇標(biāo)準(zhǔn)中，并且權(quán)重的自適應(yīng)選擇使得最后的主成分?jǐn)?shù)量恰當(dāng)，進(jìn)而使估計的偏差與方差之間能達(dá)到相對平衡。

1 模型與算法

1.1 函數(shù)型回歸模型與主成分基展開

假設(shè)響應(yīng)變量Y為標(biāo)量，X(t)是定義在某個區(qū)間T上的函數(shù)型協(xié)變量，滿足E(X(t))=0，E(X2(t))<∞?？紤]一元函數(shù)型線性回歸模型：

(1)

式中：t∈T，不失一般性，可取T為[0,1]；β0為未知的截距項；β(t)為未知的光滑斜率函數(shù)；ε為零均值、方差有限的隨機(jī)誤差項，且與X(t)獨立。

本文采用主成分基函數(shù)對函數(shù)型變量X(t)和函數(shù)系數(shù)β(t)進(jìn)行擴(kuò)展。首先，函數(shù)型變量的協(xié)方差矩陣有一個譜分解：

其中，Φ(s,t)表示協(xié)變量的協(xié)方差函數(shù)，λ1≥λ2≥…≥0是算子Φ的有序特征值，滿足

{φk(t)∶k=1,2,…}構(gòu)成均方可積函數(shù)空間L2(T)中的一組標(biāo)準(zhǔn)正交基，即函數(shù)型數(shù)據(jù)X(t)的主成分基函數(shù)，再根據(jù)Karhunen-Loeve定理有

對X(t)與β(t)都用主成分基展開，則式(1)變?yōu)椋?/p>

(2)

用最小二乘法最小化如下目標(biāo)函數(shù)：

(3)

由于實際應(yīng)用中隨機(jī)函數(shù)X(t)是未知的，所以需要通過觀測到的樣本{X1(t),X2(t),…,Xn(t)}來估計未知的特征根λk和主成分得分ξik，故需要估計出協(xié)方差函數(shù)Φ。

于是改變?yōu)樽钚』缦履繕?biāo)函數(shù)：

(4)

在式(4)中，為了達(dá)到降維的目的，需要根據(jù)一些準(zhǔn)則確定主成分基函數(shù)的個數(shù)。

1.2 函數(shù)型主成分截斷個數(shù)的確定

1.2.1 基于PVE和PAVE的主成分選擇

首先通過方差解釋百分比PVE和關(guān)聯(lián)變異解釋百分比PAVE分別截斷出函數(shù)型主成分個數(shù)K1、K2。

1.2.2 權(quán)重α的確定

對上述確定的截斷個數(shù)K1和K2分別取權(quán)重(1-α)、α后相加，得(1-α)K1+αK2，但由于截斷數(shù)應(yīng)為整數(shù)，因此采用以下公式對(1-α)K1+αK2取整，得到截斷個數(shù)K：

(5)

式中：Θ[· ]表示向下取整。

權(quán)重α從(0,1)區(qū)間選擇，每一個權(quán)重依據(jù)式(5)都有一個截斷個數(shù)K，即為對應(yīng)的主成分基函數(shù)個數(shù)。將函數(shù)系數(shù)β(t)和函數(shù)型協(xié)變量X(t)均通過主成分基展開后，再對函數(shù)系數(shù)進(jìn)行估計，得到估計誤差

挑選出使估計誤差最小的權(quán)重α，若得到的權(quán)重不止一個，則取接近于0. 5的值，因為α接近于0.5不會使截斷個數(shù)過大或者過小。過大的截斷個數(shù)可能會使數(shù)據(jù)出現(xiàn)過擬合的現(xiàn)象，而過小的截斷個數(shù)可能會使得估計結(jié)果損失過多的樣本信息，因此權(quán)重表示為

該權(quán)重對應(yīng)的K值即為通過本文方法最終得到的主成分截斷個數(shù)，相應(yīng)地將目標(biāo)函數(shù)式(4)轉(zhuǎn)換為：

(6)

2 蒙特卡羅模擬研究

2.1 模擬數(shù)據(jù)

2.2 權(quán)重

本文方法需要選擇出自適應(yīng)的權(quán)重，不妨令α∈[0.1,0.2,…,0.9]。通過最小二乘法對回歸系數(shù)進(jìn)行估計，用均方誤差(RMSE)評價估計誤差。

通過對以上情形的模擬得到不同權(quán)重對應(yīng)的函數(shù)系數(shù)估計誤差，由于篇幅所限，圖1僅給出閾值為0.95的模擬結(jié)果，圖中實心點表示對應(yīng)的權(quán)重使得均方誤差RMSE達(dá)到最小，如果實心點多于1個，則靠近0.5的權(quán)重即為所求的α值。

從模擬結(jié)果可知：在大約55.56%的情形中得到的權(quán)重小于0.5，尤其是樣本量n=50的情形,α全部小于0.5，權(quán)重小于0.5說明此時在本文方法中方差解釋百分比所起的作用較大；在大約27.78%的情形中權(quán)重大于0.5，其中以樣本量n=500、σ2=1的情形為主，說明此時關(guān)聯(lián)變異解釋百分比所起的作用更大；在大約16.66%的情形中權(quán)重為0.5，此時方差解釋百分比和關(guān)聯(lián)變異百分比在本文方法中起著同樣的作用。

2.3 三種方法的比較

將分別采用三種方法的模擬結(jié)果進(jìn)行對比，主成分截斷個數(shù)的平均值如表1所示。可以看出，除去個別情形外，通過方差解釋百分比截斷的主成分個數(shù)都多于基于關(guān)聯(lián)變異解釋百分比的截斷個數(shù)，而采用本文方法得到的截斷個數(shù)處于二者之間，主成分?jǐn)?shù)量比較恰當(dāng)，對于函數(shù)協(xié)變量可以避免出現(xiàn)過擬合或者損失過多樣本信息的情況。

對自適應(yīng)加權(quán)截斷法的模擬結(jié)果進(jìn)行統(tǒng)計：當(dāng)閾值分別為0.95、0.97、0.99時，對應(yīng)的主成分基截斷個數(shù)總平均值分別為5.828、6.101、7.287，即隨著閾值的增大，主成分基個數(shù)也隨之增多；當(dāng)σ2分別為0.5和1時，對應(yīng)的主成分基截斷個數(shù)總平均值分別為6.353和6.457，即隨著誤差的增大，主成分基個數(shù)也有增大的趨勢；然而，當(dāng)樣本量n分別為50、100、500時，對應(yīng)的平均截斷個數(shù)為7.114、6.419、5.684，即隨著樣本量的增多，通過本文方法截斷的主成分個數(shù)在減少。

(a)n=50,σ2=0.5 (b)n=100,σ2=0.5 (c)n=500,σ2=0.5

表1 三種方法的主成分截斷個數(shù)平均值

為了進(jìn)一步說明自適應(yīng)加權(quán)截斷法的優(yōu)越性，比較三種方法對函數(shù)系數(shù)β(t)的估計，得到表2的結(jié)果。從表2可以看出，本文方法的估計誤差均小于其他兩種方法的對應(yīng)值，表明其對函數(shù)系數(shù)的估計精確度有所提高。對自適應(yīng)加權(quán)截斷法的結(jié)果進(jìn)行統(tǒng)計：隨著閾值由0.95逐漸增至0.99，β(t)的平均估計誤差也隨之增大，從0.3409到0.3759再到0.5032；當(dāng)σ2分別等于0.5和1時，β(t)的平均估計誤差分別為0.2665和0.5466，即平均估計誤差隨著隨機(jī)誤差的增大而增大；當(dāng)樣本量n分別為50、100、500時，β(t)的平均估計誤差分別為0.7974、0.3656、0.0567，即隨著樣本量的增多，函數(shù)系數(shù)的估計誤差在減少。

表2 三種方法對函數(shù)系數(shù)β(t)的估計誤差RMSE

2.4 自適應(yīng)加權(quán)截斷法的擬合結(jié)果

將本文方法擬合的曲線與函數(shù)型回歸曲線β(t)進(jìn)行比較，得到閾值分別為0.95、0.97、0.99的結(jié)果。限于篇幅，圖2僅給出v=0.95時的擬合效果，由對比結(jié)果來看，在各種情形下，本文方法都能很好地擬合函數(shù)型回歸曲線β(t)。

(a)n=50,σ2=0.5 (b)n=100,σ2=0.5 (c)n=500,σ2=0.5

3 實際數(shù)據(jù)分析

本節(jié)通過彌散張量成像(diffusion tensor imaging, DTI)數(shù)據(jù)對新方法進(jìn)行展示和比較分析。Goldsmith等[11]、田茂再等[12]都曾對該數(shù)據(jù)進(jìn)行過研究，原始數(shù)據(jù)可以在Huang等[13]編寫的R軟件包“Refund”中下載。

彌散張量成像數(shù)據(jù)是醫(yī)學(xué)研究數(shù)據(jù)，該數(shù)據(jù)集包含了382個樣本。通過核磁共振成像記錄大腦中水分子沿胼胝體方向彌散的軌跡，再根據(jù)部分各向異形圖在白質(zhì)束上的位點，得到93個數(shù)據(jù)。同時，實驗還記錄水分子沿著右皮質(zhì)脊髓束方向彌散的軌跡，得到55個數(shù)據(jù)，由此可以得到兩個函數(shù)型變量。另外，每位患者都做了一個有節(jié)奏的聽覺系列加法測試(paced auditory serial addition test, PASAT)，測試結(jié)果為0到60不等的得分。本文分析的目的是量化各向異形圖沿胼胝體方向?qū)憫?yīng)變量PASAT分?jǐn)?shù)的影響。

將382個樣本數(shù)據(jù)進(jìn)行預(yù)處理，對余下的236個樣本數(shù)據(jù)建立一元函數(shù)型線性回歸模型，以胼胝體數(shù)據(jù)為函數(shù)型協(xié)變量，以PASAT得分為標(biāo)量型響應(yīng)變量。在閾值分別為0.95、0.97、0.99的情況下，通過自適應(yīng)加權(quán)截斷法獲取使估計誤差RMSE最小時的權(quán)重，如圖3所示。由圖可見，閾值為0.95時，權(quán)重取0.2和0.3可以使RMSE最小；閾值為0.97時，權(quán)重取0.2、0.3及0.4使得RMSE最小；閾值為0.99時，權(quán)重取0.8和0.9使得RMSE最小。因此，對于彌散張量成像數(shù)據(jù)，在本文提出的主成分基個數(shù)的自適應(yīng)加權(quán)截斷法中，閾值為0.95和0.97時，方差解釋百分比的貢獻(xiàn)要大于關(guān)聯(lián)變異解釋百分比，而閾值為0.99時恰好相反。

(a)v=0.95 (b)v=0.97 (c)v=0.99

比較三種方法對彌散張量成像數(shù)據(jù)的預(yù)測效果，同樣采用RMSE指標(biāo)計算估計誤差，結(jié)果見表3。由表3可見，本文方法的估計誤差要稍小于基于PVE和PAVE的估計誤差，表明新方法可以提高預(yù)測精度，并且本文方法截斷的函數(shù)型主成分個數(shù)介于其他兩種方法截斷的主成分個數(shù)之間，不會因主成分過多導(dǎo)致過擬合，也不會因主成分太少導(dǎo)致?lián)p失過多的樣本信息，進(jìn)而使得估計的偏差與方差達(dá)到相對平衡。

表3 三種方法對于DTI數(shù)據(jù)的預(yù)測誤差和截斷結(jié)果

4 結(jié)語

本文針對函數(shù)型數(shù)據(jù)主成分基函數(shù)的個數(shù)選取提出了一種新的截斷方法。首先通過方差解釋百分比和關(guān)聯(lián)變異解釋百分比準(zhǔn)則對函數(shù)型主成分進(jìn)行排序和截斷，分別得到截斷個數(shù)K1、K2，然后對K1、K2進(jìn)行加權(quán)求和，并通過優(yōu)化算法獲得使估計誤差達(dá)到最小的最優(yōu)權(quán)重，進(jìn)而得到最終的主成分展開截斷項數(shù)。新方法不僅考慮了函數(shù)型協(xié)變量由于特征值迅速衰減導(dǎo)致的主要變異，還將協(xié)變量和響應(yīng)變量之間的關(guān)聯(lián)納入選擇標(biāo)準(zhǔn)中，且權(quán)重α的選擇使得最后的截斷個數(shù)不會過大或過小，進(jìn)而使估計偏差與估計方差處于一個相對平衡的狀態(tài)。

蒙特卡羅模擬結(jié)果顯示，在不同的樣本量、隨機(jī)誤差和閾值條件下，本文方法對未知函數(shù)系數(shù)的估計精度和穩(wěn)定性整體上都優(yōu)于只用方差解釋百分比準(zhǔn)則和只用關(guān)聯(lián)變異解釋百分比準(zhǔn)則的主成分截斷方法。同時，隨著閾值或者隨機(jī)誤差的增大，函數(shù)型主成分的截斷個數(shù)也在增多，并且函數(shù)系數(shù)的估計誤差也有隨之增大的趨勢；但函數(shù)型主成分的截斷個數(shù)以及對函數(shù)系數(shù)的估計誤差隨著樣本量的增大而減小。

本文方法在彌散張量成像數(shù)據(jù)分析中的實際應(yīng)用也顯示，其不僅能量化各向異形圖沿著胼胝體方向?qū)ASAT分?jǐn)?shù)的影響，而且在預(yù)測效果上也優(yōu)于其他兩種方法。