縱向數據分位回歸模型的降維算法模擬研究

羅幼喜 ，李翰芳

（1.湖北工業(yè)大學a.理學院；b.產品質量工程研究院，武漢430068；2.華中師范大學 數學與統(tǒng)計學學院，武漢430079）

0 引言

近些年來，分位數回歸[1]已經成為一種極為廣泛的數據建模方法，它能全面刻畫給定一組解釋變量時響應變量的條件分布。由于分位數提供了比平均值更完整的關于響應變量的描述，分位回歸也可以獲得給定協(xié)變量時響應變量的整個條件分布，并在不同分位點處得到響應變量關于協(xié)變量效應的總體評估[2]。對于在實際生活中經常碰到的縱向數據，如臨床試驗、小組研究、流行病學調查以及計量經濟中的面板數據等，Koenker[3]最先給出了相關討論研究，并指出對于這類數據應該充分考慮到由于對相同個體測量而產生的個體變異性，以避免給參數估計帶來的偏差。

目前關于縱向數據的分位回歸模型多是使用非拉普拉斯分布（ALD）。這類方法的關鍵在于它將分位回歸損失函數的最小化與基于ALD分布的似然函數最大化相等價關聯(lián)，使得分位回歸估計可以在似然函數的框架下進行。Koenker和Machado（1999）[4]首次提出了基于誤差項為ALD分布的分位回歸模型擬合優(yōu)度檢驗。利用這個結論，Yu[5，6]提出了基于貝葉斯框架的分位回歸方法。Luo等（2012）[7]通過ALD和分位數回歸之間的聯(lián)系，重點探索了包含隨機效應的縱向數據貝葉斯分位回歸模型的應用。然而，這些基于貝葉斯框架的分位回歸法一方面對誤差分布假設過于苛刻，不利于復雜數據分布的捕捉，另一方面后驗分布密度函數較為復雜，只能通過MCMC算法求解，計算量較大。Koenker（2004）[3]考慮了含隨機截距的縱向數據分位回歸模型，通過在分位回歸損失函數基礎上對隨機截距實施L1懲罰來估計固定效應參數，這種L1懲罰也稱為Lasso[8]。由于L1懲罰的性質，Lasso能給予施加懲罰的參數連續(xù)收縮和自動變量選擇，從而控制由大量隨機截距引入帶來的變異性。Zou（2006）[9]在Lasso基礎上引入了一種自適應Lasso方法，該方法對回歸系數使用自適應加權L1懲罰，與Lasso方法對所有系數實施相同的懲罰不同，自適應Lasso對較大系數選擇較小權重，而對較小系數選擇較大權重，所得估計具有oracle性質。最近，眾多學者利用Lasso懲罰解決了線性和非線性分位回歸模型中的變量選擇問題，如Wang和Song（2011）[10]考慮了加速失效模型的自適應Lasso程序；Yang和Liu（2016）[11]對具有缺失協(xié)變量的線性模型，提供了在分位數回歸檢驗函數基礎上實施加權自適應Lasso懲罰的復合分位回歸法等。但對于復雜縱向數據分位回歸模型中的變量選擇問題則研究還較少，本文擬將Lasso懲罰和自適應Lasso懲罰應用于縱向數據分位回歸模型的變量選擇問題中，并對方法進行模擬比較研究。

1 模型與方法

Koenker（2004）[3]考慮了含隨機截距的縱向數據分位回歸模型：

在檢驗損失函數基礎上提出了對個體效應施加L1懲罰的分位回歸方法，即極小化如下目標函數：

但上述方法是在固定模型中解釋變量基礎上進行討論的，然而，在實際問題中，人們往往需要對初始建模時給定的一系列相關解釋變量進行挑選，只保留較為重要的變量在模型之中，這樣不僅能減少模型的冗余度，也可以提高模型的預測能力。

本文將對含多重隨機效應的縱向數據分位回歸模型考慮其自變量選擇問題，記響應變量的條件τ分位函數如下：

其中 β 為 k維 τ分位回歸系數，αi，i=1，2，…，n為個體隨機效應系數。與Koenker(2004)[3]不同的是，本文考慮對固定效應系數β分別實施Lasso和自適應Lasso的懲罰分位回歸法，從而將非重要自變量權重系數壓縮至0，起到變量選擇的作用。即分別極小化：

式（5）中相合估計，l=1，2，…，k。

式（4）所得估計稱為Lasso懲罰分位回歸估計，式（5）所得稱為自適應Lasso懲罰分位回歸估計。

事實上，Lasso方法也可以從貝葉斯角度獲得一個很好的解釋[8]，即在對回歸系數 β進行先驗假設時，可以根據高維數據稀疏性的特點，假設絕大多數回歸系數值都集中在0左右，而少數非零系數則以較大概率出現(xiàn)，即先驗為比正態(tài)分布更具有尖峰厚尾的條件拉普拉斯分布[12,13]。從而可以通過賦予回歸系數條件拉普拉斯先驗構造與Lasso方法等價的Bayesian Lasso法。

假設在給定 αi條件下 Yij，j=1，…，ni，i=1，…,n 相互獨立且服從分布ALD(+σ，τ)，即：

再對回歸系數 βl，l=1，…，k賦予條件拉普拉斯先驗信息：

根據Bayes定理，可以推得回歸系數β的后驗密度函數：

不難看出極大化式（8）以得到回歸系數的貝葉斯估計，這等價于極小化式（9）：

記 λ=則求解式（9）與求解式（4）是等價的，從而β本文獲得了Lasso估計法在貝葉斯角度的一個解釋，而不難看出其中的懲罰調整參數則與“信噪比”有關。

事實上，如果假定 βl， l=1，…，k有如下獨立的條件Laplace先驗：

則類似可以得到與自適應Lasso懲罰分位回歸估計相應的貝葉斯解釋。

2 算法與懲罰參數選取準則

2.1 迭代算法

本文以Lasso懲罰分位回歸估計為例來說明如何求解式（4），算法如下：

（1）給定初始值=0，i=1，2，…，n，求解一般Lasso分位回歸估計可得=argminL(β，0)。

（2）按照下面兩步交替迭代，m=0， 1，…

①=argminαL(，α)，此步可看成是調整殘差為r(m)=yij-的一般分位回歸估計求解問題。

②=argminβL(β，)，該步驟可以先將被解釋變量yij調整為yij-，然后再求解普通的帶Lasso懲罰的分位回歸估計。

l體迭代過程中設定的誤差最大容忍上限。

類似的，本文也可以利用上述迭代算法求解式（5）獲得自適應Lasso分位回歸估計，只需要將其中求解一般Lasso分位回歸估計改為求解一般自適應Lasso分位回歸估計即可，而自適應Lasso估計可以利用線性規(guī)劃算法來解決[14]，考慮到線性規(guī)劃算法可以很容易給出不同懲罰參數下的估計路徑，從而本文只需根據給定的懲罰參數選取準則在每步迭代中挑選最優(yōu)的參數值即可。而對于縱向數據一般分位回歸估計，Koenker（2005）[2]早已給出了較為成熟的求解算法。值得一提的是，在自適應Lasso求解過程中，是未知的，本文可以取均值回歸模型下的最小二乘估計絕對值的倒數來代替。

2.2 SIC參數調優(yōu)準則

本文將采用SIC準則（Schwarz Information Criterion）來解決式（5）求解過程中的最優(yōu)懲罰參數λβ選取問題，該準則既考慮到了盡量使模型擬合度高，又不至于讓模型太復雜。SIC準則計算方法如下：

其中為檢驗函數殘差和，衡量了模型的擬合度；是樣本總量，| M |是模型中非零回歸系數的總數，衡量了模型的簡潔度。不難看到當λβ使得式（11）達到最小時，模型的擬合度和簡潔度可以取得一個平衡。在具體的計算過程中，可以先確定一個搜索區(qū)間并將區(qū)間進行等份，然后在所有的等份點上利用交叉驗證的辦法得到最優(yōu)的懲罰參數。

3 蒙特卡羅模擬比較

3.1 生成模型設定

為了檢驗本文提出的Lasso懲罰分位回歸估計和自適應Lasso懲罰分位回歸估計的效果，本文利用下面的縱向數據模型來模擬實際數據：

在上述生成數據模型中，本文設置各變量和參數如下：

（1）固定效應部分：8個自變量 X1，…，X8均來自N(0，1)分布，且任兩個自變量之間的相關性隨著其下標差絕對值增大而減少，具體的 ρXlXk=ρ|l-k|，且 ρ=0.5；

（2）隨機效應部分：=(1，xij1，xij2，…，xij5) ，αi=

（3）隨機誤差部分：考慮隨機誤差εij分別來自N(0，1)，t(3)及Cauchy(0，1)分布，其中 t(3)，Cauchy(0，1)為厚尾分布，尤其是Cauchy(0，1)分布，容易產生絕對值較大的異常點，從而可以檢驗方法各種方法對異常點的穩(wěn)健性。

（4）樣本量大?。簄=30， m=10；

（5）模型稀疏度：模型截距β0=0，考慮各個自變量前系數分別稠密、稀疏和高度稀疏三種情形，其中:

①稠密情形：

β=(0.85， 0.85， 0.85， 0.85， 0.85， 0.85， 0.85， 0.85)T

②稀疏情形：β=(3， 1.5， 0， 0， 2， 0， 0， 0)T

③高度稀疏情形：β=(5， 0， 0， 0， 0， 0， 0， 0)T

3.2 比較指標設定

對于蒙特卡羅模擬結果和進行比較的方法，本文作如下設定：

（1）相互比較的幾種方法：①普通分位回歸估計,記為QRE；②只考慮隨機截距的Lasso懲罰分位回歸估計，記為IQRE；③對所有固定效應系數帶Lasso懲罰的分位回歸估計，記為LQRE；④對所有固定效應系數帶自適應Lasso懲罰的分位回歸估計，記為ALQRE。

（2）各種估計法估計精確度衡量指標：

Σ=為第 s次模擬過程中 β 的系數估計值，另外本文還給出MSE在N次模擬中的標準差SD。

（3）各種估計方法對自變量選擇準確度衡量指標：

易見，上述三個準確度衡量指標的取值范圍均為[0，1]，其值越接近1準確度越高。且三個指標分別反映不同的側面，其中Acc1反映估計方法保留重要自變量的準確性，Acc2反映估計方法排除非重要自變量的準確性，Acc為前兩者的平均，反映估計方法的綜合選擇準確性。

3.3 比較結果分析

表1至表3給出了4種方法在不同模型下重復100次的模擬結果，其中LQRE和ALQRE兩種方法在迭代過程中最大誤差容忍上限值ε=10-4，所有方法對于系數的估計絕對值若小于10-6則計為估計值為0，也即該方法將對應的自變量排除在模型之外。另外，本文對高中低分位點均進行了模擬，考慮到不同分位點結果變化不大，故此處只對τ=0.5時的結果進行詳細分析討論。

表1 稠密模型下四種方法參數估計結果

首先從表1稠密模型下的估計精度MSE來看，此時IQRE法和QRE法占優(yōu)，尤其是IQRE法，在 N(0，1)和t(3)分布下MSE值均遠遠小于其他方法。這一點也不難理解，首先，由于本模型是稠密模型，所以從理論上來說，所有固定系數均不需要進行壓縮，由于LQRE和ALQRE均對固定效應系數實施了壓縮，從而導致其估計精度較沒有壓縮的IQRE估計和QRE估計有所降低。而對比LQRE和ALQRE，在各種分布下，LQRE均優(yōu)于ALQRE，原因在于本模型中設定的所有固定效應系數均一致為0.85，而在ALQRE中，初始給定的權重估計很難完全一致，從而對每個變量的壓縮程度也會有所不同，所以此時將所有系數壓縮權重進行同等對待的LQRE法較ALQRE法占優(yōu)，然而在實際問題中，這種所有解釋變量系數均相等的情況是較為少見的。另外，由于本模型中所有變量均為重要解釋變量，所以只需比較各方法的重要解釋變量保留能力指數Acc1即可。從Acc1的值來看，顯然完全沒有進行變量挑選的QRE法和IQRE法始終將所有變量保留在模型之中，從而Acc1指數能夠達到最大值1。而本文提出的LQRE和ALQRE法由于對變量進行了選擇，故導致Acc1能力指數有所下降，但可以看到LQRE法能在各種分布下均保持90%以上的準確率，ALQRE法同樣由于采用了不等權重壓縮較LQRE法稍差，但其準確率也都在80%以上。另外總體上當分布逐漸由正態(tài)分布N(0，1)變?yōu)楹裎卜植糃auchy(0，1)時，各種方法估計精度均有所下降，但下降幅度不大，尤其是本文的LQRE法，精度還略有提升，說明方法對異常點具有較強的穩(wěn)健性。

表2 稀疏模型下四種方法參數估計結果

從表2稀疏模型下的結果來看，首先各種方法估計精度MSE的值較之前的稠密模型更為接近了，其中在N(0，1)分布下QRE法最優(yōu)，ALQRE法次優(yōu)；t(3)分布下IQRE法最優(yōu)，ALQRE法次優(yōu)；Cauchy(0，1)分布下ALQRE法最優(yōu)，IQRE法次優(yōu)。所以綜合來看，本文提出的ALQRE法能夠在各種分布下均有較好的表現(xiàn)。再從模型選擇能力指數上來看，雖然QRE和IQRE法的Acc1指數均為1，但Acc2指數均為0，也即它們沒有進行任何變量選擇，而是直接將所有變量納入到模型之中。本文提出的LQRE和ALQRE法則顯然既有重要變量選擇功能，也有冗余變量剔除功能，尤其是ALQRE法，保留重要解釋變量能力的Acc1指數均在0.99以上，且在3種分布下變量選擇能力綜合指數Acc均是最優(yōu)的。對比稠密模型，這一點也不難理解，因為此時固定效應系數不再完全相等，故如果對其進行壓縮懲罰，其懲罰權重也應該不同，LQRE法雖然進行了壓縮，但其權重系數均一樣，沒有如同ALQRE法那樣具有自適應性的功能。另外，當分布逐漸從N(0，1)變?yōu)楹裎卜植糃auchy(0，1)時，雖然ALQRE估計精度MSE有所降低，但其變量選擇能力幾乎不受任何影響，綜合能力指數Acc甚至還有提高。所以，對于這種在實際問題中可能最常見的模型，本文的ALQRE法無論是在估計精度上還是在變量選擇能力上均有著極好的表現(xiàn)。

在實際問題中，在建模前往往存在大量的冗余自變量，即高度稀疏模型。從表3的模擬結果來看，LQRE法和ALQRE顯然比傳統(tǒng)的QRE和IQRE法在變量選擇準確性以及系數估計精確性上都有較大優(yōu)勢，尤其是在變量選擇準確性方面，傳統(tǒng)方法幾乎已經失效，而LQRE法和ALQRE都保持著80%以上的綜合準確率，ALQRE更是保持85%以上的綜合準確率。在MSE值方面，LQRE法和ALQRE也在各種誤差分布下一直占優(yōu)，尤其是在Cauchy(0，1)分布之下，ALQRE的MSE值僅有QRE和IQRE估計的一半，說明本文提出的ALQRE估計即使在高度稀疏模型下也對誤差分布具有很強的穩(wěn)健性。

表3 高度稀疏模型下四種方法參數估計結果

4 結論

（1）本文將Lasso和自適應Lasso方法推廣至含多重隨機效應的縱向數據分位回歸模型中來，不僅能夠對模型中重要自變量進行自動選擇，也能夠對冗余變量進行自動剔除，極大地方便了實際工作者在數據初始建模過程中遇到的變量篩選問題。

（2）本文對兩種方法設計的參數估計交替迭代算法充分利用了現(xiàn)有分位回歸估計求解算法，快速有效地解決了目標函數優(yōu)化問題。對于Lasso和自適應Lasso中懲罰參數選取這一難題，給出了能夠平衡模型擬合效果和模型簡潔度的最優(yōu)懲罰參數選取準則，使得實際工作者不需要逐一試探懲罰參數的調節(jié)。

（3）通過在不同稀疏度情況下的模擬比較發(fā)現(xiàn)，當模型需要進行變量選擇時，本文提出的兩種方法不僅能夠對于模型中存在的冗余變量進行一定程度的剔除，而且能夠極為準確地保留模型中的重要解釋變量。在模型高度稀疏時，兩種懲罰分位回歸法變量選擇綜合準確率達到80%以上，尤其是自適應Lasso懲罰法，不僅參數估計精度最高，而且變量選擇準確度也最高，達到85%以上。所以該方法對于存在大量待篩選變量的情形尤為有效。

（4）通過在不同誤差分布情況下的模擬比較發(fā)現(xiàn)，本文提出的兩種方法對誤差分布有較強的穩(wěn)健性，尤其是在重要解釋變量選擇和冗余變量排除方面幾乎不受到影響。

參考文獻：

[1] Koenker R,Bassett,G.Regression Quantiles[J].Econometrica,1978,(46).

[2] Koenker R.Quantile Regression[M].New York：Cambridge University Press,2005.

[3] Koenker R.Quantile Regression for Longitudinal Data[J].Journal of Multivariate Analysis,2004，(91).

[4] Koenker R，Machado J.Goodness of Fit and Related Inference Pro?cesses for Quantile Regression[J].Journal of the American Statistical Association,1999,(94).

[5] Yu K,Moyeed R A.Bayesian Quantile Regression[J].Statistics and Probability Letters,2001,(54).

[6] Yu K，Stander J.Bayesian Analysis of a Tobit Quantile Regression Model[J].Journal of Economics,2007,(137).

[7] Luo Y,Lian H,Tian M.Bayesian Quantile Regression for Longitudi?nal Data Models[J].Journal of Statistical Computation and Simulation,2012，(82).

[8] Tibshirani R J.Regression Shrinkage and Selection via the Lasso[J].Journal of the Royal Statistical Society,Ser.B,1996,(58).

[9] Zou H.The Adaptive Lasso and Its Oracle Properties[J].Journal of the American Statistical Association,2006,(101).

[10] Wang X，Song L.Adaptive Lasso Variable Selection for the Acceler?ated Failure Models[J].Communication in Statistics-Theory and Methods,2011，(40).

[11] Yang H，Liu H.Penalized Weighted Composite Quantile Estima?tors With Missing Covariates[J].Statistical Papers,2016，57(1).

[12] Yuan M,Lin Y.Efficient Empirical Bayes Variable Selection and Es?imation in Linear Models[J].Journal of the American Statistical As?sociation,2005，(100).

[13] Park T，Casella G.The Bayesian Lasso[J].Journal of the American Statistical Association,2008，(103).

[14] Li Y，Zhu J.L1-norm Quantile Regressions[J].Journal of Computa?tional and Graphical Statistics,2008，(17).