錢 彤,朱永忠

(河海大學 理學院，江蘇 南京 211100)

0 引言

考慮高維稀疏線性回歸問題:

y=Xβ+ε，

(1)

其中y表示響應變量，y∈Rn，X表示回歸矩陣x∈Rn×p，β則表示估計的回歸系數(shù)向量，β∈Rp，其中p為預測變量數(shù)，n為樣本觀測數(shù).另外，ε=(ε1,…,εn)是誤差項，假設εi～N(0,σ2).在此回歸模型中，pn，且假設p個變量中只有少數(shù)變量是“相關變量”，即只有少數(shù)變量的回歸系數(shù)遠離零值.這里記 “相關變量”個數(shù)為pn，本文只考慮pn=o(n)的情況.近年來，大數(shù)據(jù)帶來的此類高維稀疏模型的估計已然成為許多領域(基因?qū)W，計量經(jīng)濟學，臨床醫(yī)學等)的重要問題.在此類問題中，由于樣本協(xié)方差XTX奇異，常見的頻率論方法，諸如最小二乘估計(OLS)，Lasso，Adaptive Lasso算法等很難處理.有效的變量選擇方法很有必要.

在貝葉斯框架下，常見的變量選擇方法包括隨機搜索變量選擇(SSVC)[1]，貝葉斯正則化(包括貝葉斯Lasso[2]，貝葉斯Adaptive Lasso[3]等)，spike-and-slab先驗變量選擇[4]，以及全局-局部(global-local)收縮先驗變量選擇，例如Dirichlet-Laplace(DL)先驗[5]，廣義doubel-Pareto(GDP)先驗[6]，Horseshoe先驗[7]以及Horseshoe+先驗[8]等.這些變量選擇方法都是依賴計算出的后驗包含概率(對于每個變量或每個模型)的閾值化，來選擇最終納入模型中的變量，從而得出最終稀疏模型.

Bondell(2012)[9]提出的貝葉斯懲罰后驗置信區(qū)域變量選擇法，將模型擬合與變量選擇2個步驟分離，依賴于后驗置信區(qū)域來尋找最稀疏模型.該方法只需要設置1個全模型的先驗∏(β)，而避免了對模型空間設置先驗.為了便于分析，該方法在模型擬合步驟使用了一個共軛正態(tài)先驗：

β|σ2,τN～N(0,σ2/τN),

(2)

其中σ2是誤差方差，τN是精確度參數(shù)，表示先驗精度與誤差精度的比率，控制收縮程度.通常給參數(shù)σ2一個逆Gamma超先驗，給參數(shù)τN一個固定值或1個Gamma超先驗.

傳統(tǒng)的全局-局部(global-local)收縮先驗類可以表示成正態(tài)先驗的全局-局部尺度混合[10]：

(3)

其中τ是全局收縮參數(shù)，控制向量β的整體收縮.λj是局部收縮參數(shù)，控制βj的不同程度收縮.π(·)表示適當?shù)南闰灧植?全局-局部先驗類的密度函數(shù)在零點處的“峰值”與其“重尾”的特點，使得在處理稀疏問題時，往往能夠“重”壓縮小系數(shù)，“輕”壓縮大系數(shù).其中，Bhadra(2015)針對“超稀疏信號勘測”提出了Horseshoe+先驗.該先驗是著名的Horseshoe先驗的拓展，并且在Horseshoe先驗的理論優(yōu)勢基礎上，具有更高的后驗集中度，對小系數(shù)的壓縮更 “重”，對大系數(shù)“輕”壓縮甚至不壓縮，使大系數(shù)估計的后驗均方誤差更小.因此在高度稀疏模型中的變量選擇效果更好.

為此，本文改進貝葉斯懲罰置信區(qū)域法，基于正態(tài)先驗在稀疏問題中的表現(xiàn)不佳，將原方法中的正態(tài)先驗替換成Horseshoe+先驗，這樣做的優(yōu)勢在于：一是Horseshoe+先驗的貝葉斯分層模型能夠進行塊更新，從而可以應用Gibbs抽樣;二是提高模型擬合產(chǎn)生的后驗集中度，從而提高該方法的變量選擇效果.

1 模型及方法

1.1 懲罰置信區(qū)域法

Bondell等于2012年給出了懲罰置信區(qū)域法的詳細介紹.

(4)

(5)

(6)

其中cα是與α有關的值.對τ設置一個先驗就不能產(chǎn)生“橢圓型”后驗分布了.但是，我們?nèi)匀豢梢允褂谩皺E圓型”輪廓來產(chǎn)生置信區(qū)域，盡管對應的后驗置信區(qū)域不再是最高后驗密度區(qū)域，但是依然是合理的.因為“橢圓型輪廓”也是最高后驗密度區(qū)域的一個合理的近似.

使用一個L0與L1之間的光滑同倫[11]以及局部線性近似[12]來解決上述目標函數(shù)求解時間復雜度高以及可能存在不唯一解的問題，得出最終的最優(yōu)化問題：

(7)

其中λα是與α,cα一一對應的值.隨著置信水平α的變化，將會產(chǎn)生一個預測變量的有序集，其中，越“重要”的變量越排在集合的前面.

1.2 Horseshoe+分層模型

將Bhadra(2015)提出的Horseshoe+先驗應用到線性回歸模型中，為了使估計過程更加穩(wěn)健，考慮響應變量y=(y1,y2,…,yn)與回歸系數(shù)β=(β1,β2,…βp)的抽樣分布模型如下：

(8)

其中π(σ2)表示對誤差方差構(gòu)造的先驗，C+(0,κ)是一個標準半柯西分布，概率密度函數(shù)為：

(9)

由于半柯西先驗的后驗分布本身缺乏標準的形式，Enes(2015)針對Horseshoe先驗提出了一種直觀的抽樣方法，通過添加一個輔助變量φj，半柯西分布可表示成逆Gamma分布的尺度混合[14]：

(10)

利用逆Gamma分布的共軛性，可以直接通過Gibbs抽樣計算后驗分布，大大減少了計算復雜度.為了給每一個參數(shù)構(gòu)建一個具有封閉形式的條件后驗分布，沿用上述方法，進一步添加2個潛在變量γ,ωj，可對參數(shù)λj,τ采用如下的全條件分層結(jié)構(gòu)：

(11)

1.3 Gibbs抽樣

考慮到在所有參數(shù)的先驗分布已知的情況下，各參數(shù)的全條件分布容易抽樣，于是可通過Gibbs抽樣來推導出各參數(shù)的全條件后驗分布：

全局收縮參數(shù)τ與局部收縮參數(shù)λj以及四個潛在變量ηj,φj,γ,ωj的全條件后驗分布根據(jù)共軛性服從逆Gamma函數(shù)：

5)抽樣γ|τ2，其全條件后驗分布為逆Gamma分布IG(a5,b5)，其中a5=1,b5=1+1/τ2.

2 數(shù)據(jù)仿真

Bondell(2012)已經(jīng)證明了基于正態(tài)先驗的懲罰置信區(qū)域方法在變量選擇效果上優(yōu)于SSVS，Lasso，adaptive Lasso，Dantzig Selector以及SCAD等方法.在此基礎上，本文主要比較基于不同先驗的懲罰置信區(qū)域法的變量選擇效果以及預測精度.通過在不同稀疏程度的高維線性回歸的情況下比較基于Horseshoe+先驗，Normal先驗的懲罰置信區(qū)域方法的變量選擇效果和預測精度，為增加對比，我們同時還考慮了Bayesian Lasso與Lasso的情況.其中Bayesian Lasso相當于給系數(shù)β一個Laplace先驗：β|σ,τL～DE(σ/τL)，式中DE(a)表示一個均值為0的Laplace分布.為便于分析，以上所有先驗中的誤差方差σ2統(tǒng)一分配一個逆Gamma先驗IG(0.001,0.001).與Bondell(2012)實驗相同，給Normal先驗中的超參數(shù)τN分配一個Gamma先驗Gamma(0.01,0.01)，且Bondell(2012)通過實驗模擬證明超先驗Gamma(a,a)中的參數(shù)a對實驗結(jié)果不敏感.同理，Laplace先驗給超參數(shù)τL一個”無”信息先驗Gamma先驗Ga(1,1).另外，Lasso算法選用Efron提出的5折交叉檢驗的LARS算法.利用R軟件來實現(xiàn)后驗計算中的Gibbs抽樣.

針對(1)式線性模型，

2.1 變量選擇效果

為了衡量上述方法得出的變量排序的準確性，考慮變量排序每一步的模型準確度.對排序的每一步，引入指標TP(True Positive)、FP(False Positive)，TN(True Negative)，F(xiàn)N(False Negative)分別表示正確選擇進模型的真相關變量、錯誤選擇進模型的假相關變量、正確排除在模型外的真不相關變量、錯誤排除在模型外的假不相關變量.進而引入指標TPR(True Positive Rate)或稱specificity，F(xiàn)PR(False Positive Rate)或稱1-specificity，其中TPR等于TP占真實相關變量的比例，F(xiàn)PR等于FP占真實不相關變量的比例.以FPR為x軸，TPR為y軸，繪制ROC(Receiver Operating Characteristic)曲線.根據(jù)ROC曲線下的平均AUC(Area Under ROC Curve)面積值進行比較.通常，AUC的值介于0.5與1之間，較大的AUC值表示變量選擇效果更好.對于仿真一和仿真二的AUC值見表1～表5.

表1 p=50對應的2個仿真實驗中的AUC值

表2 p=200對應的2個仿真實驗中的AUC值

表3 p=500對應的2個仿真實驗中的AUC值

表4 p=1000對應的2個仿真實驗中的AUC值

表5 p=2000對應的2個仿真實驗中的AUC值

根據(jù)表格中的數(shù)據(jù)，可以得出結(jié)論：當真實線性回歸模型中的相關變量回歸系數(shù)產(chǎn)生自均勻分布U(0,1)，即大多接近0時，隨著維度的增加，Horseshoe+先驗的AUC值與Normal先驗的以及Laplace的AUC值相差不大.且在相關系數(shù)較高的情況下Horseshoe+先驗表現(xiàn)出略微的劣勢.當然，Lasso算法仍然存在處理相關性高的數(shù)據(jù)時效果差的弊端；當真實回歸模型中的重要回歸系數(shù)產(chǎn)生自均勻分布U(0,5)，即包含接近0的系數(shù)以及遠離0的系數(shù)，這種情況也更接近于現(xiàn)實數(shù)據(jù).在這種情況下，Horseshoe+先驗的優(yōu)勢更加突出，即使在相關系數(shù)高的情況下，依然表現(xiàn)出絕對的優(yōu)勢.這一結(jié)論是符合預期的，因為在理論上，Horseshoe+先驗面對高維稀疏回歸問題，對大系數(shù)“輕”壓縮，對小系數(shù)“重”壓縮的特點，使得在該類問題中產(chǎn)生更加集中的后驗分布，對大系數(shù)估計的后驗均方誤差更小.尤其是在仿真二實驗中，存在較大的回歸系數(shù)的情況下，該優(yōu)勢更加突出.從而在懲罰置信區(qū)域的方法中更易找出準確的最稀疏解，變量選擇效果更好.總的來說，基于Horseshoe+先驗的懲罰置信區(qū)域法在高維稀疏回歸問題中的變量選擇效果，尤其是在真實回歸系數(shù)同時包含小系數(shù)與大系數(shù)的回歸問題中，優(yōu)于Normal先驗與Laplace先驗.

2.2 預測精度

為了驗證該變量選擇方法產(chǎn)生的相關變量排序應用到預測中的效果，以及變量排序每一步的準確性，我們考慮仿真一以及仿真二實驗中，p={50,500,200 0}，ρ={0.5,0.9}的情況，獨立產(chǎn)生滿足仿真條件的20個測試集，根據(jù)50個訓練集產(chǎn)生的相關變量排序，依次選取不同數(shù)量的預測變量納入預測模型內(nèi)，計算出20個測試集的平均預測均方根誤差RMSEP，并繪制出以上模擬實驗中20個測試集的平均預測均方根誤差.見圖1～圖12.

(12)

圖1 仿真一p=50,ρ=0.5平均預測均方根誤差 圖2 仿真一p=50,ρ=0.9平均預測均方根誤差

圖3 仿真二p=50,ρ=0.5平均預測均方根誤差 圖4 仿真二p=50,ρ=0.9平均預測均方根誤差

圖5 仿真一p=500,ρ=0.5平均預測均方根誤差 圖6 仿真一p=500,ρ=0.9平均預測均方根誤差

3種先驗對應的曲線都呈現(xiàn)出一個遞減趨勢，且在區(qū)間(5,10)間遞減速率開始下降，在區(qū)間(10,20)之間遞減速率繼續(xù)下降，趨于平穩(wěn).這是符合預期的，因為相關變量只有10個，因此將10個變量納入最終模型時，大約已包含了所有相關變量，且在納入大于10個變量時，預測均方根誤差基本沒有太大變化.遞減速率逐漸下降說明產(chǎn)生的相關變量排序是按照重要性的大小，越重要的變量排在越前面.

當p=50時，Horseshoe+先驗對應的預測均方根誤差曲線與Normal先驗與Laplace先驗差異不大.隨著p逐漸增大，即稀疏度越來越高，Horseshoe+先驗對應的預測均方根誤差曲線逐漸低于其他兩個先驗的誤差曲線，預測效果更好.當然，在仿真一實驗中，當相關系數(shù)ρ=0.9時，Horseshoe+先驗的預測誤差曲線只在前半部分低于其他兩個先驗的誤差曲線，說明基于Horseshoe+先驗的懲罰置信區(qū)域法在處理相關系數(shù)高的數(shù)據(jù)時，選擇出全部相關變量的效果不夠好,但在前半部分選擇出較重要變量的能力更好.一個合理的解釋正如Juho(2016)指出[15]，當樣本數(shù)據(jù)之間存在強相關性時，Horseshoe+先驗中的全局參數(shù)τ通過樣本數(shù)據(jù)被強識別的能力較弱.但是，在仿真二實驗中，基于Horseshoe+先驗的預測誤差總是低于其他2個先驗的誤差，且當選擇納入模型的變量數(shù)接近10時，預測誤差達到最低，之后趨于平穩(wěn)，說明幾乎準確選擇出所有相關變量.總的來說，基于Horseshoe+先驗的懲罰置信區(qū)域法的預測效果優(yōu)于Normal先驗與Laplace先驗.

圖7 仿真二p=500,ρ=0.5平均預測均方根誤差 圖8 仿真二p=500,ρ=0.9平均預測均方根誤差

圖9 仿真一p=2000,ρ=0.5平均預測均方根誤差 圖10 仿真一p=2000,ρ=0.9平均預測均方根誤差

圖11 仿真二p=2000,ρ=0.5平均預測均方根誤差 圖12 仿真二p=2000,ρ=0.9平均預測均方根誤差

3 結(jié)語

在高維稀疏線線性回歸問題中，變量選擇這一步驟很重要.本文針對這一問題，在懲罰置信區(qū)域變量選擇方法的基礎上，利用該方法在高維線性回歸問題中的優(yōu)異表現(xiàn)，以及Horseshoe+先驗在高維稀疏變量選擇問題中的應用，將兩者有效結(jié)合.將Horseshoe+先驗應用到該方法中，取代原方法使用的共軛正態(tài)先驗.利用Horseshoe+先驗在高維稀疏問題中對大系數(shù)“輕”壓縮，對小系數(shù)“重”壓縮的特點，產(chǎn)生更好的后驗集中度與后驗估計，從而提高變量選擇準確性.通過數(shù)據(jù)模擬實驗結(jié)果可以得出結(jié)論，在高維稀疏回歸問題中，基于Horseshoe+先驗的懲罰置信區(qū)域法相比于正態(tài)先驗與Laplace先驗，變量選擇效果與預測效果總體上更好.

當然也存在一些不足，本文只將基于Horseshoe+先驗的懲罰置信區(qū)域變量選擇方法與原始方法中的正態(tài)先驗以及Laplace先驗進行了對比，實際上，懲罰置信區(qū)域方法可以與多種先驗結(jié)合，為此，可進一步研究基于不同先驗的懲罰置信區(qū)域方法，從而提高原方法的變量選擇效果及預測效果.