周 霖，羅幼喜

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

隨著數(shù)據(jù)收集和存儲(chǔ)能力的不斷進(jìn)步，數(shù)據(jù)獲取的方式和途徑越來(lái)越多，數(shù)據(jù)之間的關(guān)系也越來(lái)越復(fù)雜.如何針對(duì)高維復(fù)雜數(shù)據(jù)進(jìn)行統(tǒng)計(jì)建模分析是一個(gè)亟待解決的問(wèn)題.原始的變量選擇方法由于其方法本身的局限性，在面對(duì)高維復(fù)雜數(shù)據(jù)時(shí)顯得力不從心.近些年，隨著對(duì)變量施加正則化懲罰理論的提出，不同的懲罰方法被運(yùn)用到變量選擇的領(lǐng)域內(nèi).Ridge回歸采用L2正則化懲罰，通過(guò)嶺跡的圖像來(lái)判斷冗余變量，從而降低變量之間共線(xiàn)性的影響.Tibshurani[1]提出了LASSO (least absolute shrinkage and selection operation)回歸，采用L1正則化懲罰，在估計(jì)回歸系數(shù)時(shí)自動(dòng)對(duì)變量進(jìn)行選擇，將冗余變量的系數(shù)壓縮為0.LASSO懲罰函數(shù)為凸函數(shù)，且是有偏估計(jì)，通過(guò)減小預(yù)測(cè)模型的方差達(dá)到降低預(yù)測(cè)總誤差的目的，但容易存在過(guò)度壓縮系數(shù)的缺點(diǎn).為了進(jìn)一步降低預(yù)測(cè)總誤差，F(xiàn)an和Li[2]提出了SCAD懲罰，它不僅對(duì)于很小的回歸系數(shù)壓縮到0，而且是一種近似無(wú)偏估計(jì).Zhang[3]提出的MCP懲罰(minimax concave penalty)，采用對(duì)回歸系數(shù)有差別的懲罰，在保留SCAD懲罰優(yōu)點(diǎn)的同時(shí)得到更加精確的估計(jì).基于此，越來(lái)越多的學(xué)者將MCP懲罰運(yùn)用在各個(gè)領(lǐng)域研究中，在此基礎(chǔ)上進(jìn)行優(yōu)化，并在圖像去噪[4]和圖像重建[5]中有了較大的進(jìn)展.

在混合效應(yīng)模型中，考慮到樣本的數(shù)據(jù)有一定的關(guān)聯(lián)性，在模型中添加了隨機(jī)效應(yīng)項(xiàng).如果僅考慮固定效應(yīng)而忽視隨機(jī)效應(yīng)的影響，則會(huì)給最后的估計(jì)帶來(lái)偏差，但是無(wú)關(guān)的隨機(jī)效應(yīng)太多則會(huì)導(dǎo)致求解過(guò)程中隨機(jī)效應(yīng)的協(xié)方差陣不可逆，不利于估計(jì)，因此如何解決個(gè)體的隨機(jī)效應(yīng)所帶來(lái)的誤差是高維混合效應(yīng)模型進(jìn)行變量選擇的關(guān)鍵.

將混合效應(yīng)模型和分位回歸聯(lián)合起來(lái)，在不同的分位點(diǎn)進(jìn)行變量選擇是一個(gè)新興的研究領(lǐng)域.現(xiàn)有的文獻(xiàn)中，多數(shù)學(xué)者都是對(duì)隨機(jī)效應(yīng)進(jìn)行L1正則化[6-8]，該方法只能剔除冗余的隨機(jī)效應(yīng)，減少隨機(jī)效應(yīng)的數(shù)量，但是無(wú)法選擇出固定效應(yīng)中的重要預(yù)測(cè)變量.羅幼喜[9]針對(duì)此問(wèn)題提出來(lái)一種雙Lasso正則化分位回歸法(DLQR)，對(duì)模型中的固定效應(yīng)和隨機(jī)效應(yīng)同時(shí)施加L1懲罰，使得模型的冗余度大大降低.但是L1正則化得到的結(jié)果為有偏估計(jì)，本文以MCP正則化替代L1正則化，提出雙MCP正則化分位回歸法(DMQR).

1 模型與方法

傳統(tǒng)的線(xiàn)性混合效應(yīng)模型為：

在實(shí)際生活中，數(shù)據(jù)往往存在異方差或者尖峰、后尾等情況，使得傳統(tǒng)的均值回歸效果不好，且一條回歸線(xiàn)所反映的信息也是有限的，分位回歸相對(duì)于均值回歸，其應(yīng)用條件更加寬松，依據(jù)不同的分位數(shù)，對(duì)自變量進(jìn)行回歸，能更好的描述自變量對(duì)因變量的變化范圍以及條件分布形狀的影響.Koenker和Bassett[10]提出線(xiàn)性分位回歸：

Q(τ)=inf(y：F(y)≥τ)，

其中，inf表示最大下界，Koenker[6]考慮了含隨機(jī)截距的條件分位回歸模型：

i=1，2，…，n，

(1)

并提出對(duì)個(gè)體波動(dòng)施加L1正則化壓縮的分位回歸方法，數(shù)學(xué)表達(dá)式如下：

(2)

在實(shí)際問(wèn)題中，有時(shí)個(gè)體效應(yīng)不僅對(duì)模型截距有影響，也可能對(duì)模型的斜率有影響，將上述回歸模型推廣到既包含隨機(jī)截距又包含隨機(jī)斜率的條件分位回歸模型，即在給定個(gè)體隨機(jī)效應(yīng)αi條件下，響應(yīng)變量的τ分位回歸函數(shù)為：

j=1，2，…，ni；i=1，2，…，n，

(3)

羅幼喜[9]提出了DLQR方法，對(duì)個(gè)體效應(yīng)施加L1懲罰以防止模型的過(guò)擬合，同時(shí)對(duì)固定效應(yīng)施加L1懲罰來(lái)剔除冗余變量，選出重要的自變量，減少模型的冗余度，即極小化下式：

(4)

其中，λα為對(duì)個(gè)體效應(yīng)施加的懲罰參數(shù)，λβ為對(duì)固定效應(yīng)施加的懲罰參數(shù).ρτ(·)為分位回歸函數(shù)，其表達(dá)式如下：

Lasso懲罰函數(shù)為凸函數(shù)，且是有偏估計(jì)，通過(guò)減小預(yù)測(cè)模型的方差達(dá)到降低預(yù)測(cè)總誤差的目的，但容易存在過(guò)度壓縮系數(shù)的缺點(diǎn).為了進(jìn)一步降低預(yù)測(cè)總誤差所提出的SCAD懲罰，它不僅對(duì)于很小的回歸系數(shù)壓縮到0，而且是一種近似無(wú)偏估計(jì).之后提出的MCP懲罰，采用對(duì)回歸系數(shù)有差別的懲罰，在保留SCAD懲罰優(yōu)點(diǎn)的同時(shí)得到更加精確的估計(jì).MCP懲罰函數(shù)的數(shù)學(xué)表達(dá)式如下：

(5)

上式中，λ為正則參數(shù)，用來(lái)調(diào)整懲罰力度；a為調(diào)整參數(shù)，用來(lái)控制懲罰范圍.其懲罰力度為：

(6)

當(dāng)|β|≤aλ時(shí)，MCP的懲罰力度會(huì)隨著參數(shù)絕對(duì)值的增大而減小，當(dāng)|β|=0時(shí)其懲罰力度最大為λ，當(dāng)|βi|=aλ時(shí)，其懲罰力度變?yōu)?，當(dāng)回歸系數(shù)繼續(xù)增大時(shí)其懲罰力度保持為0不變.將Lasso與MCP懲罰的懲罰力度可視化如下圖1所示.

圖1 Lasso與MCP懲罰力度圖Fig.1 Lasso and MCP penalty intensity diagram

從上圖可以看出，Lasso對(duì)回歸系數(shù)采用恒定的懲罰力度，而MCP對(duì)回歸系數(shù)采取有差別的懲罰，當(dāng)|βi|≥aλ，其懲罰力度為0，其參數(shù)的估計(jì)類(lèi)似于最小二乘估計(jì)，其擬合的系數(shù)結(jié)果更精確.基于此，本文將DLQR法中L1懲罰替換為MCP懲罰，提出了雙MCP正則化的分位回歸方法(DMQR)，即極小化下式：

(7)

其中，Pλβ(|βl|)表示對(duì)固定效應(yīng)施加MCP懲罰，Pλα(|αit|)表示對(duì)個(gè)體效應(yīng)施加MCP懲罰.

2 回歸系數(shù)估計(jì)方法與正則化參數(shù)的選取準(zhǔn)則

2.1 回歸系數(shù)估計(jì)方法

針對(duì)高維非凸懲罰的分位回歸問(wèn)題，Wang和Peng[11]提出了迭代坐標(biāo)下降(QICD).QICD算法將MM算法與坐標(biāo)下降算法相結(jié)合.具體地說(shuō)，首先用一個(gè)優(yōu)化函數(shù)來(lái)替代非凸懲罰函數(shù)，并得到新的目標(biāo)函數(shù).然后，將關(guān)于目標(biāo)參數(shù)的每次和每次循環(huán)替代所有參數(shù)，直到收斂為止.對(duì)于每個(gè)坐標(biāo)下降步驟，只需要計(jì)算單變量加權(quán)中位數(shù)，可確?？焖儆?jì)算.具體計(jì)算流程圖見(jiàn)表1．雙MCP正則化分位回歸估計(jì)可以采用交替迭代算法求解，即每次只求解單個(gè)MCP正則化分位回歸，具體迭代算法求解見(jiàn)表2.

表1 QICD 算法偽碼表Tab.1 QICD algorithm pseudo code table

表2 雙MCP懲罰分位回歸迭代算法Tab.2 Double MCP penalty quantile regression iterative algorithm

2.2 正則化參數(shù)選取準(zhǔn)則

選擇合適的正則化參數(shù)對(duì)于MCP懲罰的效果很關(guān)鍵.在DMQR法中一共有4個(gè)正則化參數(shù)，即Pλβ(|βl|)中的aβ和λβ，Pλα(|αit|)中的aα和λα，為簡(jiǎn)化計(jì)算量和降低模型的復(fù)雜度，一般建議aα=aβ=3.0[12].此時(shí)需要選擇的正則化參數(shù)只有兩個(gè)，即調(diào)整固定效應(yīng)的懲罰力度的λβ和調(diào)整個(gè)體效應(yīng)懲罰力度的λα.運(yùn)用上述迭代算法，每次迭代過(guò)程中，只用考慮單個(gè)正則化參數(shù)，使得參數(shù)選取的復(fù)雜度大大降低.

Wang[13]在Fan和Li[14]的基礎(chǔ)上提出了具有一致性的BIC法用于正則化參數(shù)選取.其數(shù)學(xué)表達(dá)式為：

(8)

(9)

其中，

(10)

3 計(jì)算機(jī)模擬比較研究

3.1 稀疏模型下兩種懲罰方法的比較

數(shù)據(jù)生成的公式如下：

(11)

表3給出了τ=0.25，0.5，0.75三個(gè)不同分位點(diǎn)下的模擬結(jié)果，重復(fù)模擬次數(shù)為100次，衡量模型精度的均方誤差(MSE)定義如下：

(12)

C1=重要預(yù)測(cè)變量被正確選擇個(gè)數(shù)/真實(shí)重要預(yù)測(cè)變量總個(gè)數(shù)，

C2=冗余預(yù)測(cè)變量被正確排除個(gè)數(shù)/真實(shí)冗余預(yù)測(cè)變量總個(gè)數(shù).

C1和C2取值介于0和1之間，C1越接近1 表示其挑選重要預(yù)測(cè)變量的能力越強(qiáng)，反之越差；C2越接近1表示其排除冗余預(yù)測(cè)變量的能力越強(qiáng)，反之越差.表3 還給出了100次重復(fù)模擬中3 個(gè)重要預(yù)測(cè)變量X1、X2、X5被正確選擇和5 個(gè)冗余預(yù)測(cè)變量X3、X4、X6、X7、X8被正確排除的總次數(shù).其中，判斷是否選擇了某變量是根據(jù)對(duì)該變量系數(shù)估計(jì)是否為0 得到的(本次模擬基于R 3.6.2，platform為x86＿64-w64-mingw32條件下實(shí)驗(yàn)運(yùn)行).

表3 不同條件下兩種懲罰分位回歸結(jié)果Tab.3 Regression results of two penalty quantiles under different conditions

從變量選擇的角度分析，在不同的σ取值下，DLQR和DMQR都能準(zhǔn)確的選擇出模型中的重要預(yù)測(cè)變量，但是在排除冗余變量方面，DLQR最多只能準(zhǔn)確排除75%左右的冗余變量，而DMQR則能正確排除90%以上的冗余變量，使模型中的干擾變量更少，具體效果見(jiàn)圖2.

圖2 變量選擇效果Fig.2 Variable selection effect

從估計(jì)的精度來(lái)看，在τ=0.5；σ=1時(shí)，DMQR的估計(jì)精度最優(yōu).且在不同分位點(diǎn)，不同的σ取值時(shí)均比同等情況下的DLQR估計(jì)精度要高，且方差更小，具體效果見(jiàn)圖3.

圖3 兩種懲罰精度對(duì)比圖Fig.3 Comparison of two types of penalty accuracy

保存每次DMQR法模擬產(chǎn)生的30個(gè)6維行向量αi(i=0，1，2，3，4，5)，得到一個(gè)30×6的α矩陣，即為隨機(jī)效應(yīng)α的估計(jì)值，計(jì)算每一列的均值和方差作為隨機(jī)效應(yīng)參數(shù)的一個(gè)估計(jì).在不同分位點(diǎn)和不同信噪比下取100次模擬的平均值作為最后隨機(jī)效應(yīng)α均值和方差的估計(jì)值，具體數(shù)據(jù)如下表4(上方數(shù)據(jù)為均值，括號(hào)內(nèi)的數(shù)字為方差).

表4 隨機(jī)效應(yīng)估計(jì)表Tab.4 Random effect estimate table

從均值的角度分析，在不同的分位點(diǎn)和不同的信噪比下，隨機(jī)效應(yīng)系數(shù)α的均值都接近于0，擬合效果較好；從方差的角度分析，在不同分位點(diǎn)和不同信噪比下，α4，α5方差的擬合效果都很好，方差也接近于0，其他隨機(jī)效應(yīng)系數(shù)的方差，隨著信噪比σ的增大，模擬效果逐漸降低，在不同分位點(diǎn)，當(dāng)σ=1時(shí)其模擬效果最好，最接近于1，且在τ=0.5；σ=1，模擬效果最佳.

τ=0.5時(shí)，在不同信噪比下，保存每次模擬最后一步懲罰參數(shù)λα，λβ的值，100次模擬下，懲罰參數(shù)的箱線(xiàn)圖如圖4所示.

圖4 懲罰參數(shù)選取箱線(xiàn)圖Fig.4 Penalty parameter selection box plot

從上圖可以看出，隨機(jī)效應(yīng)的懲罰參數(shù)λα其分布較為密集，其上四分位點(diǎn)與下四分位點(diǎn)包含的區(qū)間更小，且懲罰參數(shù)基本集中在0.09左右，固定效應(yīng)的懲罰參數(shù)λβ其上四分位點(diǎn)和下四分位點(diǎn)所包含的區(qū)間更大，且其最大值與最小值之間的差距相較于λα更大.

3.2 不同誤差分布時(shí)估計(jì)的比較

下文在不同誤差分布情形下比較DMQR同其他幾種分位回歸方法在變量選擇和估計(jì)精度上的表現(xiàn)，數(shù)據(jù)生成的模型為式(11)，自變量系數(shù)設(shè)置為稠密、稀疏、高度稀疏三種不同情況：

1)稠密模型：β=(0.85，0.85，0.85，0.85，0.85，0.85，0.85，0.85)；

2)稀疏模型：β=(3，1.5，0，0，2，0，0，0)；

3)高度稀疏模型：β=(5，0，0，0，0，0，0，0).

取σ=1，τ=0.5，D=diag(2，2，2，2，0，0)，取誤差εij分別來(lái)自N(0，1)，t(3)及Cauchy(0，1)分布.比較的方法有：①普通分位回歸(QR)；②雙Lasso懲罰分位回歸(DLQR)；③雙MCP懲罰分位回歸(DMQR).表4、表5、表6給出了各種情況下重復(fù)100次的模擬結(jié)果.

表5 稠密模型下模擬結(jié)果Tab.5 Simulation results under dense model

表6 稀疏模型下模擬結(jié)果Tab.6 Simulation results under sparse model

在稠密模型下，所有的變量均為重要預(yù)測(cè)變量，所以不存在正確排除冗余變量，故C2不存在.從MSE角度來(lái)看，不論誤差服從何種分布對(duì)系數(shù)進(jìn)行懲罰后的模型在精度上面不如直接進(jìn)行分位回歸的模型，DMQR模型和DLQR模型的擬合精度相差不大，主要原因是在稠密模型下，所有變量都不需要進(jìn)行懲罰，故直接分位回歸(QR)得到的結(jié)果是最好的.

從誤差服從的不同分布來(lái)看，當(dāng)誤差服從正態(tài)分布時(shí)，三種估計(jì)方法的精度最高，當(dāng)誤差服從柯西分布時(shí)，估計(jì)的誤差最大，且對(duì)系數(shù)進(jìn)行懲罰后，C1的值不能完全達(dá)到1.

在稀疏模型下，普通的分位回歸(QR)在排除模型中冗余變量的能力為0，因?yàn)槠胀ǚ治换貧w不能對(duì)變量進(jìn)行選擇，只能將所有的變量均視為重要變量保留在模型中.而對(duì)系數(shù)施加懲罰后的方法則能排除模型中的冗余變量，但是不同懲罰所得到的結(jié)果也不同.

從排除模型中冗余變量的能力來(lái)看，DMQR剔除冗余變量的效果最好，能剔除90%以上的冗余變量，且在誤差服從標(biāo)準(zhǔn)正態(tài)分布時(shí)的效果最好，達(dá)到了94%.而DLQR只能剔除約70%以上冗余變量，從估計(jì)的精度來(lái)看，DMQR在3種方法中的估計(jì)精度最高，且在誤差為正態(tài)分布時(shí)的精度最高，當(dāng)誤差為柯西分布時(shí)的精度最低，具體效果如圖5.

圖5 稀疏模型下擬合估計(jì)效果圖Fig.5 Fitting estimation effect diagram under sparse model

在高度稀疏模型下，無(wú)論誤差服從什么分布，DLQR法能剔除約80%以上的冗余變量，DMQR法能剔除約90%的冗余變量，雖然DMQR法的效果仍比DLQR法的結(jié)果好，但兩者在剔除冗余變量上的差距相較于稠密模型和稀疏模型縮小了許多.在模型的精度上面，DMQR法要優(yōu)于DLQR法，DMQR法在正態(tài)分布時(shí)精度最高，誤差為0.188，說(shuō)明DLQR法雖然在高度稀疏模型下能夠很好的排除冗余變量，但是對(duì)重要變量的擬合精度不高，擬合出來(lái)的結(jié)果與真實(shí)結(jié)果差距較大，擬合效果如圖6.

圖6 高度稀疏模型下擬合效果圖Fig.6 Fitting effect diagram under highly sparse model

3.3 高維情況比較

本節(jié)在高維情形下考察雙正則化分位回歸DLQR和DMQR的表現(xiàn).數(shù)據(jù)的生成模型仍為式(7)，但減少樣本量至n=10，m=10，即總樣本量為100.在模擬(2)的稀疏模型中額外添加102個(gè)獨(dú)立噪聲變量X9，X10，…，X110，所有變量均獨(dú)立同分布于N(0，0.52)，在添加噪聲變量后，總的變量個(gè)數(shù)為110個(gè)，大于總樣本量.其中重要的預(yù)測(cè)變量為3個(gè)，冗余變量為107個(gè).另外設(shè)σ=0.5 ，D=diag(1，1，1，1，0，0).表7給出了兩種方法在τ=0.5 、0.9 兩個(gè)分位點(diǎn)下進(jìn)行100次模擬估計(jì)的結(jié)果.

表7 高度稀疏模型下模擬結(jié)果Tab.7 Simulation results under highly sparse model

在高維情況下，DLQR法與DMQR法兩者的模擬差距非常明顯.從變量選擇的角度看，DLQR法只能正確剔除約81%的冗余變量，而DMQR法幾乎能剔除冗余變量的能力約為90%，且在τ=0.5 時(shí)能達(dá)到93%.在估計(jì)的精度上面，DMQR則要遠(yuǎn)遠(yuǎn)好于DLQR，DLQR法在不同分位點(diǎn)的平均誤差為5.94，而DMQR法的平均誤差僅為2.82.兩種方法在τ=0.5 時(shí)估計(jì)的精度均優(yōu)于τ=0.9時(shí)的精度.

表8 高維情形不同分位點(diǎn)下估計(jì)結(jié)果Tab.8 Estimation results under different quantiles in the high-dimensional case

4 結(jié)論

本文在羅幼喜[9]提出的雙Lasso的基礎(chǔ)上改進(jìn)了懲罰方法，提出了雙MCP正則化分位回歸(DMQR).通過(guò)模擬結(jié)果可以發(fā)現(xiàn)，不論模型是稠密模型、稀疏模型還是高度稀疏模型，改進(jìn)后的方法不論是在正確選擇模型中重要變量方面還是剔除冗余變量方面均比原來(lái)的方法要好，且模型的精度更高.DMQR法對(duì)誤差有很好的穩(wěn)健性，在正態(tài)分布下擬合的效果最好.在高維情況下，DMQR相對(duì)于DLQR法的優(yōu)勢(shì)更大，基本上能完全選中模型中的重要預(yù)測(cè)變量，且排除冗余變量的能力也能達(dá)到90%左右，模型估計(jì)的精度更是遠(yuǎn)遠(yuǎn)大于DLQR法．