基于融合Lasso的非參數加性分位數回歸模型

付漫俠 周水生

分位數回歸[1]是一種用于估計感興趣的條件分位數的魯棒方法,對重尾分布具有天然的魯棒性,可正確擬合異方差關系,為響應的條件分布提供全面的分析,其稀疏度水平可隨分位數變化而變化.因此,分位數回歸在研究收入分配、股票價格波動、疾病治療效果和教育等方面應用廣泛.

加性模型的優勢在于它可捕捉數據中的復雜關系,如非線性關系.此外,可避免維數困擾的問題.維數困擾是指在高維數據中建立全局模型時,參數估計變得困難或不可靠.通過將模型分解為可加組件,加性模型可更好地處理高維數據,使參數估計更穩定、可靠.

受加性模型和廣義加性模型工作的啟發,de Gooijer等[2]提出基于核方法建模以估計加性非線性條件分位數模型,證明對于固定維數p,加性分位數回歸達到與加性均值回歸相同的收斂速度,從而在理論上減輕維數困擾,但是方法需要對維數p≥5的情況進行偏差校正.Gaillard等[3]通過B樣條擬合廣義加性分位數模型生成溫度情景,將其插入概率預測負荷模型,得到電力負荷預測的最佳效果.然而,此方法需要研究者手動選擇平滑參數.Fasiolo等[4]使用樣條基展開擬合響應和感興趣的分位數之間的平滑關系,給出條件分位數校準公式,快速自動估計平滑參數.模型結構類似于廣義加性模型,只在維數較低時進行,都是使用B樣條基函數逼近分量,樣條函數的性能很大程度上取決于節點的選擇.Sherwood等[5]提出QA-SCAD,利用組懲罰實現超高維加性分位數回歸問題估計和變量選擇的同時進行.雖然將加性分位數回歸模型擴展到超高維問題,但仍基于B樣條函數,隨著維數增加,計算速度越來越慢.此外,樣條函數的性能很大程度上取決于節點的選擇.然而,節點的選擇是主觀的,可能因為選擇的不同導致擬合結果的不同.

現有的均值回歸模型在理論和計算上的結論無法直接應用在分位數回歸上,因此,在理論方面,Padilla等[6]給出分位數回歸預測一致性的損失函數.Ye等[7]在此基礎上研究分位數K-NN Fused Lasso,證明K-NN Fused Lasso接近極小極大收斂速度.在計算方面,Pietrosanu等[8]研究高維分位數回歸問題,使用ADMM(Alternating Direction Method of Multipliers)[9]、MM(Majorize-Minimize) Algorithm[10]以及Coordinate Descent Algorithm[11]進行變量選擇,獲得估計量.Hochbaum等[12]開發O(nlogn)次的分位數融合Lasso快速算法.Brantle等[13]提出Parallelizable ADMM Algorithm,計算一維數據的k階分位數趨勢濾波估計量.

本文提出基于融合Lasso的非參數加性分位數回歸模型(Nonparametric Additive Quantile Regre-ssion Model Based on Fused Lasso, AQFL),是在融合Lasso罰和l2罰之間折衷的可對加性分位數回歸模型進行估計和變量選擇的模型.融合Lasso罰使模型能快速計算,并在局部進行自適應,實現對所需分位數甚至極端分位數的預測.同時結合l2罰,在高維數據中將對響應影響較小的協變量函數值壓縮為0,實現變量的選擇.此外,本文給出保證收斂到全局最優的塊坐標ADMM算法(Block Coordinate Alter-nating Direction Method of Multipliers, BCADMM),證明AQFL的預測一致性.最后,在合成數據和碎豬肉數據上的實驗表明AQFL在預測準確性和魯棒性等方面較優.

1 相關工作

傳統的回歸分析只關注因變量的條件均值,未充分考慮因變量條件分布的完整特征.相比之下,分位數回歸能分析因變量條件分布的完整特征.因此,本文希望將一些在均值回歸中有效的研究方法推廣到加性分位數回歸中.

1.1 均值回歸

在均值回歸模型中,學者們采用不同的方法構建復雜模型,其中一種方法是使用協變量的多項式或樣條函數引入非線性關系.使用非參數建模的方法擬合均值回歸模型也是常見的,此方法可以更靈活地適應數據的特征,不需要對模型做出具體的函數形式假設.Tibshirani等[14]提出的非參數融合Lasso罰模型,通過響應向量y和有序的單變量x估計向量

θ=(θ1,θ2,…,θn),

其中

θi=E[yi|xi]=f(xi) ,
x1

具體模型表達式如下：

其中,λ>0為調優參數,差分矩陣

上述方法將部分系數和系數之間的差收縮為0,得到模型系數的稀疏解和系數差分的稀疏解,使相鄰系數之間更平滑.

Kim等[15]將討論擴展到單變量非參數趨勢濾波罰模型,并將趨勢濾波罰函數定義為輸入數據上的k+ 1階差分算子,將融合Lasso歸納為k= 0的情況.之后,Tibshirani[16]將趨勢濾波罰與樣條方法進行詳細對比.研究發現,相比平滑樣條方法,非參數趨勢濾波方法在模型的擴展性、理論分析、算法實現、局部自適應等方面都具有優勢.

1.2 加性分位數回歸

xi=(xi1,xi2,…,xip)∈Rp

為p維向量.給定分位數τ∈(0,1).加性分位數回歸模型有如下形式：

令

省略τ,則估計的加性分位數回歸模型：

其中

F(yi|xi)表示xi條件下yi的累積分布函數,隨機誤差εi滿足

P(εi≤0|xi)=τ.

此外,令

可提升模型的可識別性.

令

其中,損失函數

非對稱絕對偏差函數

ρτ(t)=(τ-1{t≤0})t.

相比傳統的均值回歸,加性分位數回歸通過對不同分位數的目標函數進行建模,捕捉協變量在不同分位數上的非對稱性,有利于處理偏態數據和非線性關系.加性分位數回歸可靈活適應于不同的數據分布和模型假設,不需要對數據的正態分布和模型的線性關系進行假設,對異常值和離群點具有較好的魯棒性,適用于各種類型的數據和問題.此外,加性分位數回歸還可用于構造預測區間,而均值回歸只能提供不同位置的點估計.

2 基于融合Lasso的非參數加性分位數回歸模型

擬合加性分位數模型的方法通常使用樣條函數逼近分量.但樣條函數的性能很大程度上取決于預先選定的節點,如果節點數量過多或位置選擇不當,可能會出現過擬合現象.此外,樣條函數的計算復雜度相對較高,特別是在高維數據或復雜形狀下,這可能導致計算時間較長或需要更高的計算資源,因此本文提出基于融合Lasso的非參數加性分位數回歸模型(AQFL).

2.1 模型建立

融合Lasso罰是一種用于變量選擇和估計的統計方法,本文將其應用于非參數加性分位數回歸,提出如下模型：

(1)

其中,

λ≥0,Xj=(x1j,x2j,…,xnj)T,

Pj表示將向量Xj從小到大排序的置換矩陣.可將DPj作為一個整體,若DPj的第q行元素

則DPjθj的第q個元素

(DPjθj)q=θlj-θkj.

此外,融合Lasso罰函數通過對模型參數施加稀疏性懲罰,自動選擇對目標函數具有重要影響的變量,從而提高模型的解釋能力和泛化能力.這種方法將模型參數推向0附近,有效減少過擬合的風險.通過稀疏性懲罰,還可減少異常值和噪聲對估計結果的影響,提高模型的魯棒性.

在高維數據中,本文希望對響應影響較小的協變量函數值為0,進而添加一個l2罰,并在兩個罰項間折衷.由此得到如下優化問題：

(2)

其中0≤α≤1.α在鼓勵θj成為分段常數和在整個向量θj上使用組Lasso誘導稀疏性之間提供一種權衡[17],從而實現模型擬合和變量選擇的同時進行.本文稱式(2)為AQFL.

2.2 求解算法

首先考慮單個變量時如何求解估計器.當p= 1時,式(2)改為

(3)

使用ADMM求解,將式(3)重寫為

上式的增廣拉格朗日表達式為

(4)

其中R為更新中控制步長的懲罰參數.通過迭代更新下式以求解式(4)：

(5)

(6)

u←u+θ-z,

(7)

其中u為拉格朗日乘子.式(5)有如下形式的閉式解：

(8)

式(6)有如下估計[18]：

(9)

上式可使用現有的求解器,如R語言中的glmgen包進行計算.

當p>1時,使用塊坐標下降法得到全局最優解,步驟如算法1所示.

算法1BC-ADMM

step 2 循環

j= 1,2,…,p,1,2,…,p,1,2,…,p,…

step 2.1 計算殘差

step 2.2 由ADMM求解式(3)：

step 2.1.1 初始化

step 2.1.2 循環k=1,2,…,直至

或達到最大迭代次數Niter.

step 2.3 計算截距：

并居中

step 3 重復step 2,直至式(2)收斂.

3 參數分析

調優參數λ的選擇至關重要,本節討論AQFL的解不全為0的充分條件,從而得到參數λ的取值范圍,并給出選擇最優模型的兩種方法.

3.1 λ的取值范圍

在α= 0,1時,分別研究如何選擇變量λ的值,使AQFL的解

成立,從而得到λ的取值范圍.

引理1定義

λ≥‖VTγ(0)‖∞;

當α= 0時,式(2)的解全為0當且僅當

λ≥‖γ(0)‖2.

其中

證明若

為

(10)

當α= 1時,令

βj=DPjθj,

則有

令

得到式(10)的最優性條件

-VTγ(β)+λs(β)=0,

其中

γ(β)=(τ-1{y1<(Vβ)1},

τ-1{y2<(Vβ)2},…,τ-1{yn<(Vβ)n})T.

當βj≠0時,

s(βj)=sign(βj);

當βj=0時,

s(βj)∈[-1,1].

若β=0,有

-VTγ(0)+λs(0)=0

成立當且僅當

λ≥‖VTγ(0)‖∞.

當α= 0時,式(2)的最優性條件

-γ(θj)+λs(θj)=0,j= 1,2,…,p.

當θj≠0時,

當θj=0時,

s(θj)∈{g∶‖g‖2≤1}.

若θj=0,j= 1,2,…,p,

-γ(0)+λs(0)=0

成立當且僅當

λ≥‖γ(0)‖2.

由引理1可得AQFL的解不全為0的充分條件.

推論1?α∈[0,1],若式(2)的解不完全為0,則

推論1的結果給出參數λ的合適范圍.根據推論1的計算結果,最大值記為λmax,這意味著當λ值超過λmax時,整體擬合將達到完全稀疏性(即式(2)的解全為0).

因此,本文將λ設置在0.01×λmax到λmax之間,在這個范圍內選擇合適的稀疏性水平.

3.2 λ的選擇

調優參數λ的選擇是模型估計中一個重要的實際問題,控制估計器的平滑程度.AQFL的λ值可通過K折交叉驗證或基于貝葉斯信息準則(Bayesian Information Criterion, BIC)選擇.分位數回歸的BIC[19]計算公式如下：

其中,σ>0,可經驗地選擇為

v等于罰項的非零元素的個數[20].

4 理論分析

從而將均值回歸上的評估指標推廣到分位數回歸.它類似于Huber損失函數,在最小二乘代價函數和l1范數代價函數之間折衷,因此對數據中的異常值不太敏感.本文利用這一損失函數證明AQFL存在風險上界,具有較快的收斂速率.

為了簡單起見,假設y在本節中的均值為0.令{an}?R,{bn}?R為2個正序列.如果存在常數C>0和n0>0,使得對n≥n0有an≤Cbn,記

an=O(bn).

若對于?ε>0,存在C>0,使得

P(an≥Cbn)<ε

成立,記

an=OP(bn);

若

an=OP(bn),bn=OP(an),

則記

an=Θ(bn).

首先為加性分位數回歸的理論分析做出如下假設1.

假設1令

則

假設2[6]若存在一個常數L>0,使得

‖δ‖∞≤L,

其中fyi|xi為給定xi時yi的條件概率密度.

假設2確保yi對應的分位數是唯一的,并且累積分布函數在分位數的鄰域周圍呈一致的線性增長.這是對yi分布的一個普遍假設,適用于大多數常見的分布序列.通過上述假設,有如下定理1成立.

成立.特別地,當α= 1時,若

至少有1-ε/2的概率使得

成立.

定理1表明,當α= 1時,模型(2)至少有n-1/2的收斂速率(忽略對數因子).隨著α的減小,收斂速率越來越慢,但模型仍是收斂的.為了證明定理1,首先引入定義1.

定義1[6]定義經驗損失函數：

其中

令

使用定義1中的符號,考慮M-估計量：

且

假設

θ*∈S?Rnp,

其中S表示貫穿本節的一般約束集.

定義2[7]對集合S?Rnp,S的次高斯寬度定義為

其中s1,s2,…,snp為獨立的1-次高斯隨機變量.

令

Dj∶=DPj,j=1,2,…,p,

將式(2)改寫為

其中

因此,構造矩陣

D=diag{D1,D2,…,Dp}(n-1)p×np,

其廣義逆

Ξ=diag{Π1,Π2,…,Πp}np×np,

Πj表示Ker(Dj)上的投影矩陣,則

D?D=Inp-Ξ.

引理2若假設1和假設2成立,令δ∈Rnp,滿足

Δ2(Wδ)≤t2,

則

SGW({δ∶δ∈S,Δ2(Wδ)≤t2})≤

有了上述定義和引理作為基礎,現可給出定理1的證明.

證明固定μ∈[0,1],令

(11)

其中,

因此

其中,

(12)

式(12)可由文獻[7]中引理20得到,c為正常數.由文獻[21]知,至少有1-2ε的概率使得

其中c1、c2均為正常數.因此,令λ=2R2,存在常數C> 0,使得

(13)

成立的概率至少為1-2ε.

令

設函數

顯然,g為連續函數,且

g(0)=0,g(1)=q2.

有

由不等式(11)可知,

因此

故

(14)

令ε∈(0,1),假設事件

發生的概率至少為1-ε/2,則

因此根據式(14)、 馬爾可夫不等式和三角不等式,

定義

H(η)={δ∈A∶Δ(Wδ)≤η}.

若δ∈H(η)且I成立,則由式(13)可定義

則由文獻[7]的引理15、引理16和本文引理2可知,

成立,其中常數cγ>1.根據文中的符號定義,

特別地,當α=1時,選擇λ=2R1,取

有

5 實驗及結果分析

5.1 評估標準

1)均方誤差(Mean Square Error, MSE)：

2)平均絕對誤差(Mean Absolute Error, MAE)：

3)平均檢驗誤差(Mean Check Error, MCE)：

4)真變量(True Vector, TV)：正確包含在模型中的非零協變量的數目.

6)最小比(Proportion Smaller, PS)：測試集上響應真實值的最小值和其預測值的比例,PS的值應接近τ.

5.2 合成數據實驗

本節選擇如下對比方法：QGAM(Quantile Gene-

ralized Additive Models)[4]、QA-SCAD[5]和FLAM(Fu-

sed Lasso Additive Model)[18],驗證AQFL的優越性.QGAM基于秩為30的三次回歸樣條基數,利用qgam包實現.QA-SCAD使用rqPen包實現,參數設置為a= 3.7,罰函數為SCAD(Smoothly Clipped Absolute Deviation)罰,算法選擇QICD(Quantile Iterative Coordinate Descent).FLAM使用flam包實現,采用K折交叉驗證選取參數λ.

σjk=0.5|j-k|.

其次,在向量a=(a1,a2,…,ap)T∈Rp上定義

Φ(a)={Φ(a1),Φ(a2),…,Φ(ap)}T∈Rp,

5.2.1 低維合成數據實驗

從如下3種問題的模型中分別生成500個訓練樣本,并從同一模型生成1 000個測試樣本.

1)問題1：

2)問題2：

3)問題3：

其中,

圖1為數據的生成圖,其中fj為分段常數函數或在定義域的某些區域恒定其它區域高度可變的函數,用于說明AQFL的局部自適應性.進而由問題1～問題3得到服從正態分布、重尾分布和具有異方差性質的數據,說明AQFL可靈活適應于不同的數據分布和模型假設.將問題1～問題3稱為低維設置.問題1和問題2是同方差的,第τ分位數

Q(yi|xi)(τ)=f1(xi1)+f2(xi2)+f3(xi3)+f4(xi4)+Φ-1(τ),

主要關注中位數τ= 0.5的情況.問題3是異方差的,第τ條件分位數

Q(yi|xi)(τ)=f1(xi1)+f2(xi2)+f3(xi3)+f4(xi4)Φ-1(τ),

主要關注分位數τ= 0.9,0.1的情況.

分位數回歸方法可直接對分位數進行建模,但均值回歸方法不能提供非中位數的估計值,所以對于問題3生成的數據,本文只對比分位數估計器.

(a)f1(b)f2

(c)f3(d)f4

為了全面評估AQFL的性能,本文對3種不同的問題下生成的數據進行模擬實驗.在實驗中,使用測試集MSE最小時對應的調優參數λ評估AQFL.依據定理1結論,取參數α= 0.5,0.75,1,用于提升模型的收斂性.由于QGAM只適用于低維情況,所以在低維數據上計算測試集的MSE、MCE和PS值,對比AQFL、QGAM和FLAM的性能.在對中位數建模時,MAE是MCE的2倍.因此,FLAM通過MAE值除以2得到MCE值.為了可靠評估3個模型的性能,對每個模型的模擬實驗重復運行50次.

在低維情況下,3個模型在問題1～問題3上的擬合結果如表1所示.由表可知,在低維數據上,QFAL在α=1時效果最優,因此本節主要對比α=1時的實驗結果.在問題1中,FLAM擬合效果更優的原因是它擬合均值模型,而AQFL(τ=0.5)擬合中位數模型.因為正態分布的均值是其位置參數,均值回歸在正態分布的數據上通常能夠更好地捕捉數據的整體趨勢,從而產生較小的MSE.問題2中AQFL的各指標值都優于其它兩個模型,充分說明相比均值回歸模型,分位數模型更適合具有重尾分布的數據.在問題3的異方差特征的數據中,相比QGAM,AQFL在不同分位數水平上的擬合效果更優.因此,AQFL在預測方面做得更優.此外,在區域的一個部分快速變化而在另一個部分恒定的數據上,AQFL比QGAM擬合效果更優,說明AQFL具有局部自適應性.

表1 低維情況下各模型在3個問題上的擬合結果

本文使用問題1中的響應生成方式,即使用參數α= 1和MSE最小時對應的λ值繪制圖2,對比FLAM、AQFL(τ=0.5)和QGAM(τ=0.5)在4個函數上的預測曲線以及它們重復50次后的標準誤差,其中,AQFL和FLAM使用自舉法,逐點計算標準誤差.

(a1)FLAM (a2)QGAM (a3)AQFL

(b1)FLAM (b2)QGAM (b3)AQFL

(c1)FLAM (c2)QGAM (c3)AQFL

(d1)FLAM (d2)QGAM (d3)AQFL

5.2.2 高維合成數據實驗

本節考慮在5.2.1節的問題1～問題3上添加96個噪聲函數(即f5=0,f6=0,…,f100=0)的高維設置.為簡單起見,仍記為問題1～問題3,參數α和λ等其它設置與5.2.1節一致,并在高維數據的計算中增加TV和NZ兩個評估標準,以評估AQFL、QA-SCAD和FLAM預測和變量選擇的性能.為了可靠評估3個模型的性能,對每個模型的模擬實驗重復運行50次.

在高維情況下,3個模型在問題1～問題3上的擬合結果如表2所示.由表可見,在具有噪聲函數的高維數據上,AQFL在α= 0.75時的指標值最優.當α= 1時,AQFL只能將模型參數推向0附近,無法實現與α< 1的AQFL相當的稀疏性,這充分說明模型添加l2罰對變量選擇的重要性.而在α= 0.5時,l2罰函數作用增加,可能會過度懲罰一些模型參數,使原本非零的變量整體為0.

因此,本節主要對比α= 0.75時的實驗結果.FLAM在問題1上擬合效果更優,這與低維設置中的原因一致.

下面重點關注TV和NZ兩個反映模型能否正確進行變量選擇的指標.在問題1和問題2中,AQFL、FLAM和QA-SCAD都包含真正的協變量.因為只有4%的變量是真正的非零變量,而AQFL選擇很多實際為0的無關變量,所以QA-SCAD的NZ值略優于AQFL.可通過類比最小二乘設置中的Lasso來理解這種現象.當使用預測精度MSE選擇Lasso的調優參數時,Lasso在模型預測和變量選擇上實現的效果通常是不一致的[22],故本節選擇的調優參數有利于預測,但不利于變量選擇.在異方差數據的問題3中,對于TV標準,AQFL效果更優,可正確選擇非零變量,而QA-SCAD會導致系數過度收縮,不能正確選擇非零的變量,因此AQFL更適合在有維數噪聲的高維數據上實現極端分位數的預測和變量選擇.

表2 高維情況下各模型在3個問題上的擬合結果

在高維情況下,AQFL、FLAM和QA-SCAD在不同問題上生成響應數據,共計訓練50次,并對訓練時間取平均值,結果如圖3所示.針對問題3的異方差數據,FLAM不能提供非中位數的估計值,因此只對比AQFL和QA-SCAD的訓練時間.由圖可見,AQFL在訓練時間上處于中等水平.這是因為FLAM只能處理均值回歸模型,而AQFL基于分位數損失函數,需要更多的迭代步驟找到全局最優解,導致訓練時間增加.在處理加性分位數回歸問題時,相比QA-SCAD,AQFL在3個問題上都具有更高的訓練效率.因此,AQFL在實際應用中仍具有一定的優勢.

圖3 高維情況下各模型在不同問題上的訓練時間對比

5.3 碎豬肉數據實驗

為了體現AQFL優勢,在一組真實數據上進行實驗,并與QA-SCAD[5]進行對比.本文分析來自faraway包中通過測量240個碎豬肉樣品的脂肪含量和100通道吸收光譜獲得的215個樣本,刪除5個具有離群值的觀測值,留下210個樣本,從中隨機抽取200個樣本作為訓練數據,另外10個樣本作為測試數據.由于信道頻譜數據高度相關,所以使用前30個主成分轉換預測因子.主成分居中并按比例排列,使其平均值為0,標準差為1.使用30個主成分作為協變量對脂肪含量的對數進行擬合.

擬合分位數回歸模型需要選擇τ.2個τ值可用于創建一個不需要對誤差進行參數假設的預測區間.例如,τ= 0.1,0.9的模型可用于創建80%的預測區間.如果對整個條件分布感興趣,可使用更大范圍的τ值.為了估計整個條件分布,定義多個τ值檢驗AQFL,取

τ∈Τ=

{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}.

參數λ選擇如5.2節所述,計算過程重復100次.

使用MSE(τ= 0.5)、MAE(τ= 0.5)、MCE(τ∈Τ)和QB(τ∈Τ)這4個標準評估模型.從5.2.2節實驗知α= 0.75時效果最優,所以表3對比AQFL(α= 0.75)和QA-SCAD的平均值,由此說明AQFL能夠適應不同的數據類型.

表3 AQFL和QA-SCAD在對碎豬肉數據的指標值平均值

6 結 束 語

分段型數據在實際數據中很常見,因此準確估計具有這種結構的數據是一個重要的研究課題,然而在多元情況下,缺乏關于分段數據的穩健估計.因此,本文提出基于融合Lasso的非參數加性分位數回歸模型(AQFL),并給出相應算法,有效計算優化問題的解,從理論角度證明模型的預測一致性.在合成數據和真實數據上的實驗評估AQFL的性能,實驗表明,AQFL具有較優的魯棒性和局部適應性,并且研究結果具有一定的推廣價值.