基于函數型主微分與主成分的嶺回歸模型與應用

高海燕,張 悅

(1.蘭州財經大學 統計學院,甘肅 蘭州 730020;2.甘肅省數字經濟與社會計算科學重點實驗室,甘肅 蘭州 730020)

0 引言

隨著互聯網技術的快速發展以及數據收集技術的不斷進步,大量復雜連續數據應運而生.對于如何高效處理這些高頻連續數據并探究其內在規律性, Ramsay[1]提出函數型數據分析方法(Functional Data Analysis,FDA),以函數曲線的形式表示離散觀測值,運用函數曲線的研究方法剖析數據.FDA包括函數型主成分分析(Functional Principal Component Analysis,FPCA)[2]、主微分分析(Principal Differential Analysis ,PDA)[3]、函數型聚類分析[4]以及函數型回歸分析[5]等方法.

與傳統回歸模型不同,函數型回歸分析模型無過多假定和約束,其參數是關于時間的函數,適用性較為廣泛.為提高函數型線性回歸模型(Functional Linear Regression,FLR)[1]的預測能力和可解釋性,眾多學者對FLR做了拓展[6-10].如胡錫健等[11]針對汾渭平原的空氣質量數據與氣象數據,建立關于氣溫曲線與月均SO2濃度的函數型空間自回歸模型;蘇梽芳等[12]提出基于殘差函數主成分的估計方法,以預測股市開盤價;Oshinubi等[13]基于FPCA分析法國COVID-19數據,并通過建立FLR模型預測死亡人數.函數型回歸模型要求協變量或響應變量具有函數型變量,有以下4種情況[5]:

(1)協變量為函數型變量,響應變量為標量;

(2)協變量為向量,響應變量為函數型變量;

(3)協變量和響應變量都為函數型變量;

(4)協變量為函數型變量和標量的混合變量.

鑒于西班牙COVID-19數據的函數特性,本文考慮第三種情況.

在實際應用中,由于觀測遺漏、數據記錄或錄入錯誤和設備維護等,數據缺失是常見的問題[14-15],其會不同程度地增大統計分析的復雜性和難度,降低統計推斷的精度,最終導致統計分析結果偏誤.例如在創新藥開發臨床試驗中,因受試者失訪、對干預措施不耐受或缺乏療效等原因中途退出試驗等導致數據缺失,嚴重影響臨床試驗結果[16].因此,如何準確插補醫學數據中的缺失數據具有重要意義.在處理缺失數據時,通常采用刪除法、多重插補法和回歸插補法等[17].回歸插補法是最常見且有效的方法,其利用變量間的關系進行插補,使得模型更具解釋性.例如惠嬌嬌[18]構建了函數型線性空間自回歸模型和函數型線性空間誤差模型,并將其應用于西班牙氣象數據,結果表明所提兩種模型效果均優于FLR;Acal等[19]基于西班牙COVID-19數據,提出多元函數型主成分回歸模型,插補住院人數和ICU人數.

本文針對隨機缺失機制[20]下函數型響應變量缺失的數據,結合FPCA和PDA的思想,構建基于函數型主微分與主成分的嶺回歸模型(Ridge Regression Model based on Functional Principal Differential and Principal Component,FPDPCRR).該方法同時考慮原始數據曲線信息和曲線波動特征信息兩個視角,將兩視角的估計聚合為最終預測結果,解決函數型響應變量存在缺失的插補問題.以西班牙COVID-19數據為例,驗證FPDPCRR模型的估算能力,具有較好的實際應用效果.同時,典型相關分析呈現入院率與疾病反應之間存在高度相關關系,進一步說明該模型自變量和響應變量選取的合理性.

1 相關理論基礎

1.1 離散數據的曲線擬合

在離散時間點{tj}觀察一個連續可微的過程,

yj=x(tj)+εj,

(1)

其中:{yj}為觀測值;x(t)為連續可微的擬合函數;εj為誤差成分.通常使用最小二乘準則刻畫x(t)的準確性;另外,采用粗糙懲罰法保證x(t)的勻滑程度.曲線擬合方程為

(2)

其中:Dk為k階導數;λ為修勻參數,可由廣義交叉驗證準則得到.一般地,最小化式(2)可得到x(t).

1.2 多元函數型主成分回歸模型

多元函數型線性回歸模型是利用J個函數型預測變量X=(X1,…,XJ)′來估計函數型響應Y.具體模型為

(3)

其中:α(t)為截距函數;βj(s,t)為J個系數函數;εi(t)是獨立的誤差函數.

由于多元函數型線性回歸模型易受多重共線性的影響,導致參數估計精度下降.因此,Acal[19]提出了多元函數型主成分回歸模型(Multivariate Functional Principal Component Regression Model,MFPCR).

函數型預測因子和函數型響應的主成分分解分別為

(4)

其中,權重函數fl(t)是樣本協方差算子的特征函數,ξil(t)為主成分得分.由式(3)和式(4)得

(5)

通過截斷每個主成分分解,得MFPCR模型為

(6)

1.3 主微分分析

主微分分析是通過微分方程來擬合噪聲數據,以捕獲單個曲線特征或曲線的變化特征.假設對于t∈R,函數型數據x(t)有如下線性微分算子L:

Lx(t)=β0(t)x(t)+…+
βm-1(t)Dm-1x(t)+Dmx(t)=f(t),

其中:βj(t)為權重函數;f(t)為強迫函數.

由逐點最小化法[21],得β(t)的最小二乘解為

(7)

其中,Z(t)為N×m階矩陣,其第i行為Zi(t)={-xi(t),…,-Dm-1xi(t),fi(t)}.

2 基于函數型主微分與主成分的嶺回歸模型

考慮到函數型主成分代表函數曲線的絕大部分信息,而主微分旨在刻畫函數曲線的主要變化特征,如曲線的變化趨勢、梯度以及曲率等特征.同時考慮主微分與主成分兩個視角,將有助于提高函數型回歸模型的預測精度.因此,構建FPDPCRR模型,從兩個互補視角中估計預測值,并自加權的調節視角權重,獲得最終預測結果.

類似于FPCA中主成分得分,定義xi(k)的第k個PDA得分[22],

(8)

函數型預測因子和函數型響應的主微分分解為

(9)

結合式(3)和式(9),有

(10)

受MFPCR模型[19]啟發,引入PDA結果,整合兩個視角的估計值,FPDPCRR模型為

(11)

其中:γ是正則化參數,l2是正則化項,有助于提高模型的泛化能力.

特別地,當ω2=0時,FPDPCRR退化為MFPCR模型;而當ω1=0時,FPDPCRR退化為多元函數型主微分回歸模型(Multivariate Functional Principal Differential Regression Model,MFPDR):

然而,由于MFPCR模型僅考慮函數曲線的絕大多數信息,而MFPDR模型僅考慮函數曲線波動信息,理論上它們的預測插補性能次于FPDPCRR模型.為此,下面運用實例來證明上述推斷.

3 實例檢驗

3.1 數據說明

本文以西班牙17個地區的COVID-19數據為例,運用所提FPDPCRR方法對其進行預測插補.17個地區包括:安達盧西亞、阿拉貢、阿斯圖里亞斯、巴利阿里群島、加那利群島、坎塔布里亞、卡斯蒂利亞拉曼查、卡斯蒂利亞萊昂、加泰羅尼亞、埃斯特雷馬杜拉、加利西亞、馬德里、穆爾西亞、納瓦拉、派斯瓦斯科、拉里奧哈和巴倫西亞.為描述方便起見依次記為:AC1-AC17.進一步,為驗證FPDPCRR模型中自變量和響應變量選擇的合理性,對插補后的完整數據進行典型相關分析.

3.2 曲線擬合

首先,對確診人數、死亡人數、康復人數、住院人數和ICU人數5個變量進行曲線擬合并修勻,分別表示為X1(t),X2(t),X3(t),Y1(t)和Y2(t).本文選用 4階B-樣條基函數擬合函數曲線.修勻參數λ=1e-3.分別繪制上述變量的擬合曲線,如圖1所示.

3.3 FPDPCRR模型應用

基于FPDPCRR模型處理Y1(t)和Y2(t)的缺失插補預測問題.具體地,首先運用13個具有完整數據的地區估計FPDPCRR模型;其次對4個缺失地區的數據進行插補預測.

估計5個函數型變量的每個函數型主成分,結果表明第一個主成分解釋了5個變量中的絕大部分信息,即X1(t),X2(t),X3(t),Y1(t)和Y2(t)分別為99.2%、98.3%、97.9%、91.9%、93.9%.

在得到5個變量的函數型主成分之后,基于PDA思想,選取二階微分方程分別對X1(t),X2(t),X3(t),Y1(t),Y2(t)進行PDA,并利用式(8)計算得到各變量的兩個主微分得分.然后,將13個完整數據的地區看作預測樣本,分別考慮預測變量和響應變量之間主微分得分、第一主成分的相關性.因此,根據每個預測變量的主微分得分和第一主成分,將函數型線性回歸模型簡化為主微分回歸模型:

(a) X1(t)

(b) X2(t)

(c) X3(t)

(d) Y1(t)

(e) Y2(t)圖1 5個函數型變量擬合曲線圖

以及第一主成分回歸模型

其中k=1,2;i=1,…,17;j=1,2.上述模型可以基于X1(t)、X2(t)和X3(t)的主微分得分和第一主成分分別準確估計Y1(t)和Y2(t)的主微分得分和第一主成分,其決定系數如表1所列.顯然,第二主微分得分的解釋性優于第一主微分得分.因此,選取第二主微分得分和第一主成分進行函數型數據插補.

表1 X1(t),X2(t),X3(t)分別與Y1(t),Y2(t)的決定系數

(12)

為展示所提FPDPCRR模型的插補預測效果,采用均方根誤差(RMSE)、歸一化均方根誤差(NRMSE)以及均方根百分比誤差(RMSPE)作為評估指標.具體公式為

表2 住院人數曲線的預測效果評估對比

由表2、表3可知,MFPDR模型的三個評估指標均大于FPDPCRR模型和MFPCR模型,表明其插補預測能力較差.個別地區的MFPCR模型的預測效果略優于FPDPCRR模型的,但從平均結果來看,FPDPCRR模型三個評估指標的均值均小于MFPCR模型,其預測效果優于MFPCR模型.這意味著引入主微分,從數據曲線及其波動特征信息兩個互補視角中估計,可進一步提升預測性能.總的來說,三個模型的插補預測性能:MFPDR模型

此外,分別繪制由FPDPCRR模型插補預測的4個數據缺失地區 (AC7、AC8、AC11和AC12)的住院人數和ICU人數的預測曲線、觀測曲線以及其置信區間,分別如圖2、圖3所示,圖中黑色實線表示觀測曲線,灰色實線表示預測曲線,黑色虛線表示置信帶.

圖2 住院人數的觀測曲線、預測曲線及置信區間

圖3 ICU人數的觀測曲線、預測曲線及置信區間

3.4 典型相關分析

為進一步說明所提FPDPCRR模型的插補預測精度,以及模型變量選擇的科學性,下面基于典型相關分析研究入院人數與受疾病影響人數之間的關系.由于函數型變量的主微分之間、主成分之間的相關性較高,則變量非線性獨立.在不區分自變量和因變量的情況下,將典型相關分析應用于這些函數型變量的主成分與主微分加權的得分中,以解釋兩組變量之間的相關關系.

表4 典型相關系數及顯著性結果

兩組變量的標準化典型相關系數,其大小表示各變量對相應典型變量的貢獻,如表5,表6所列.典型變量為

U1=0.08019246×ξy1+0.2844273×ξy2,
U2=-2.2768425×ξy1+0.3162396×ξy2,
V1=-0.01284091×ξx1+

0.71636136×ξx2+0.02283154×ξx3,
V2=-0.14690342×ξx1-
1.12059582×ξx2-0.08767666×ξx3.

(a) 第一對典型相關變量

(b) 第二對典型相關變量圖4 典型相關變量的散點圖

表5 HOR的標準化典型系數

表6 IR的標準化典型系數

HOR、IR變量與典型變量之間的平方相關性如表7和表8所列.第二典型變量對(U2,V2)的相關性較小;而第一典型變量對(U1,V1)與除康復人數之外的其他變量的相關性均較高,表明U1與住院人數和ICU人數更相關,V1與確診人數和死亡人數更相關.

表7 HOR變量和典型變量之間的平方相關性

表8 IR變量和典型變量之間的平方相關性

上述結果表明,醫院的入院率特別取決于住院人數;同時,對于大流行病的應對能力主要取決于確診人數和死亡人數.盡管ICU人數和康復人數在典型變量中起著重要作用,但其貢獻較小.

最后進行典型冗余分析,以研究典型變量對于各變量的信息提取量情況,分析結果如表9所列.可以看出,第一典型變量對(U1,V1)提取了原始變量的大多信息,HOR和IR組的解釋方差比例分別為0.909和0.618,而第二典型變量對幾乎沒有貢獻.

表9 方差解釋率

4 結語

本文考慮FPCA能夠代表原始數據的絕大多信息,PDA既能反映曲線的變化趨勢,又能挖掘數據的梯度、曲率等潛在信息,從這兩個視角出發,針對響應變量存在缺失的函數型數據,構建FPDPCRR模型,并利用西班牙COVID-19數據評估該模型.結果證實,與MFPDR模型和MFPCR模型相比,所提FPDPCRR模型插補預測效果最優.最后,為驗證模型中變量的選取問題,基于插補預測的完整數據,運用典型相關分析解釋HOR與IR之間的相關關系,結果表明兩組變量彼此高度相關,住院人數對HOR的影響較大,確診人數和死亡人數對IR的影響較大,ICU人數和康復人數在典型變量中的貢獻較小.