基于函數(shù)型主微分與主成分的嶺回歸模型與應(yīng)用

高海燕,張 悅

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

0 引言

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

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

(1)協(xié)變量為函數(shù)型變量,響應(yīng)變量為標量;

(2)協(xié)變量為向量,響應(yīng)變量為函數(shù)型變量;

(3)協(xié)變量和響應(yīng)變量都為函數(shù)型變量;

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

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

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

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

1 相關(guān)理論基礎(chǔ)

1.1 離散數(shù)據(jù)的曲線擬合

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

yj=x(tj)+εj,

(1)

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

(2)

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

1.2 多元函數(shù)型主成分回歸模型

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

(3)

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

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

函數(shù)型預測因子和函數(shù)型響應(yīng)的主成分分解分別為

(4)

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

(5)

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

(6)

1.3 主微分分析

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

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

其中:βj(t)為權(quán)重函數(shù);f(t)為強迫函數(shù).

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

(7)

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

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

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

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

(8)

函數(shù)型預測因子和函數(shù)型響應(yīng)的主微分分解為

(9)

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

(10)

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

(11)

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

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

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

3 實例檢驗

3.1 數(shù)據(jù)說明

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

3.2 曲線擬合

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

3.3 FPDPCRR模型應(yīng)用

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

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

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

(a) X1(t)

(b) X2(t)

(c) X3(t)

(d) Y1(t)

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

以及第一主成分回歸模型

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

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

(12)

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

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

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

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

圖2 住院人數(shù)的觀測曲線、預測曲線及置信區(qū)間

圖3 ICU人數(shù)的觀測曲線、預測曲線及置信區(qū)間

3.4 典型相關(guān)分析

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

表4 典型相關(guān)系數(shù)及顯著性結(jié)果

兩組變量的標準化典型相關(guān)系數(shù),其大小表示各變量對相應(yīng)典型變量的貢獻,如表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) 第一對典型相關(guān)變量

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

表5 HOR的標準化典型系數(shù)

表6 IR的標準化典型系數(shù)

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

表7 HOR變量和典型變量之間的平方相關(guān)性

表8 IR變量和典型變量之間的平方相關(guān)性

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

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

表9 方差解釋率

4 結(jié)語

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