含順序類別自變量的中位數懲罰回歸及應用研究

吉洋瑩，潘雨辰，黃 磊

(西南交通大學 數學學院， 成都 611756)

0 引言

統計研究時，數據往往存在多個變量。近幾十年，關于變量選擇方法的研究非常多。Akaike[1]提出AIC準則，有助于減少模型過擬合的風險。Schwarz[2]考慮了樣本數量，提出了懲罰因子為log(n)的BIC準則。另外，在自變量越來越多的大數據時代下，其維度p甚至大于樣本量n。在這種高維情況下，基于AIC、BIC及其擴展的子集選擇方法將導致計算的負擔和不可靠的結果。因此，統計學家又提出了許多針對高維環境下有用的變量選擇方法。例如，Fan等[3]提出了SCAD(smoothly clipped absolute deviation)懲罰，SCAD既能連續地壓縮系數，也能在系數較大時取得漸近無偏估計。Tibshirani[4]提出了最小絕對收縮選擇算子(least absolute shrinkage and selection operator，LASSO)，在最小二乘基礎上添加l1懲罰得到稀疏模型。由于LASSO對回歸系數使用相同程度的壓縮，可能會造成估計不完全有效，模型選擇結果不一致的現象。因此，Zou[5]改進LASSO并提出自適應LASSO方法，使回歸系數有不同程度的壓縮。自適應LASSO能夠一致地識別真實模型且所得估計具有oracle性質。 除此以外，對于LASSO統計學家還做了不少拓展，針對具有組特征的數據結構，Yuan等[6]提出可從組角度進行變量選擇的group-LASSO方法。Group-LASSO對一組系數向量添加約束，克服了LASSO無法從組水平進行變量選擇的缺點。 Wang等[7]將Zou[5]和Yuan等[6]的思想加以結合修訂，又提出了針對不同組系數施加不同懲罰的自適應group-LASSO方法。

統計建模中，回歸用來量化因變量和自變量之間的關系。其中，最常見的參數估計方法是普通最小二乘法(ordinary least square，OLS)。最小二乘估計是一種基于均值回歸的估計方法，在滿足正態誤差假設時具有良好的性質，如估計量具有無偏性和相合性。然而，當誤差不滿足高斯-馬爾科夫假設條件，就有可能產生有偏的估計結果。若樣本數據中存在異常值，模型估計值可能存在較大偏差[8]。另外，有時候不僅希望研究因變量的期望均值，還希望能探索因變量的全局分布。于是，Koenker和Bassett[9]提出分位數回歸(quantile regression，QR)。分位數回歸應用條件相對寬松，可以描述因變量的全局特征，也可以挖掘更為豐富的信息，另外分位數回歸估計采用最小化加權誤差絕對值和進行估計，通常不受異常值的影響，結果更為穩健。Yu等[10]總結了一些分位數回歸典型的應用領域，分位數回歸在醫藥與生存分析、金融與經濟統計、環境分析等諸多領域得到應用。 Engle等[11]在經濟中對于計算VaR，提出了一種類似廣義自回歸條件異方差形式的計算分位數的方法。劉軍躍等[12]利用分位數回歸法，從碳減排角度實證分析了長江經濟帶在不同碳排放水平下產業結構升級與碳排放的關系。此外，分位數回歸還拓展到各種類型的數據中。例如，具有異方差誤差項的數據、時間序列數據、刪失數據等。

近年來，建立分位數回歸和LASSO類型相結合的分位數懲罰回歸的研究越來越廣泛。Ciuperca[13]提出了含有組特征變量的分位數自適應group-LASSO。Wang等[14]結合中位數(least absolute deviation，LAD)回歸和LASSO構建了LAD-LASSO。 LAD是分位數回歸中τ=0.5的回歸，LAD不受異常值的影響，結合正則化估計方法對高維數據進行變量選擇，可以得到穩健的結果。最近，許多分位數回歸方法的相關研究大多都集中在維數超過觀測值的高維數據中，如Wu等[15]，Wang等[16]研究了l1懲罰分位數回歸模型，證明了LASSO懲罰分位數回歸估計量的oracle性質。

旨在提出一種含有順序類別自變量的高維回歸方法。教育研究中，受教育水平可以分為小學、中學、大學、研究生教育[17]。患者的身體情況表現為嚴重損傷、輕微損傷、正常、非常好[18]。 這種具有順序類別的自變量，通常按照被觀測對象特征高低或大小依次賦值，將順序類別自變量轉換為啞變量形式。針對這類含有順序類別自變量的模型，利用前面提到的分位數懲罰回歸模型中的LAD-LASSO對模型進行選擇，之后若直接對模型進行擬合，不可保證能消除掉組變量中的偽分類，即組內相鄰變量系數相等。 直接進行模型擬合可能會造成過擬合的現象。為了避免發生過擬合， Walter等[19]提出順序類別變量的編碼方案，能夠很好地解釋系數估計，但依然存在過擬合和估計值不存在等問題。Gertheiss等[20]運用嶺回歸能夠得到穩健的參數估計，但過擬合現象還是沒有解決。之后，Tian等[21]提出啞變量線性變換方法(transformed dummies，TD)，將TD與BIC方法結合提出TD-BIC方法檢測偽分類。 Huang等[22]提出TD-LASSO方法，高維環境下采用LASSO較BIC更為容易。周曉霞等[23]將TD與機器學習相結合識別偽分類。晁越等[24]則是提出一種漸近相合的探測logistic回歸模型中順序多分類解釋變量的LTD-BIC偽分類識別方法。又由于LASSO的無差別懲罰，采用自適應LASSO進行改進。因此，將TD與自適應LASSO相結合提出TD-ADLASSO，再與LAD回歸相結合提出LAD-TD-ADLASSO識別偽分類。LAD-TD-ADLASSO既可以做參數估計還能識別偽分類，可以有效避免模型的過擬合風險，提高預測精度。

剩余部分安排如下：第1節詳細介紹了TD和LAD-TD-ADLASSO偽分類識別與融合技術，并整理出一個系統性的可行算法；第2節通過2個實例，根據所提出的LAD-TD-ADLASSO偽分類識別與融合技術算法展示該方法的實用性和有效性；第3節進行總結和展望。

1 模型與方法

在本節中，介紹了具有順序類別自變量的高維模型，并提出LAD-TD-ADLASSO算法。為了描述簡便，在以下模型推導中僅考慮包含一個順序類別自變量的回歸模型，而實際問題中的模型可能包含多個順序類別的自變量。考慮含有一個m分類的順序類別自變量的線性模型：

Yi=α+β2X2,i+…+βmXm,i+

α1Z1,i+…+αkZk,i+εi,

i=1,2,…,N

(1)

式中：εi,i=1,…,N為隨機誤差。實際研究中，不光有順序類別自變量，還有連續型變量，用{Z1,…,Zk}表示k個連續型變量。{X2,…,Xm}為相應的啞變量，m類別的順序類別自變量將生成m-1個啞變量，啞變量定義如下：

(2)

將模型(1)改寫成矩陣形式：

Y=Xβ+ε

(3)

式中：X為設計矩陣，X=(1,X2,…,Xm,Z1,…,Zk)，并假定列是滿秩的，Xj=(Xj,1,…,Xj,N)T,j=2,…,m。 連續型變量Zj=(Zj,1,…,Zj,N)T,j=1,…,k，因變量Y=(Y1,…,YN)T，回歸系數β=(α,β2,…,βm,α1,…,αk)T。誤差向量ε滿足高斯-馬爾科夫條件：

(ⅰ)E(ε)=0；

(ⅱ) Var(ε)=δ2I， Var(·)表示一個隨機向量的協方差矩陣。

使用最小二乘來估計回歸參數向量β，使得誤差向量ε=Y-Xβ的平方和達到最小，得到參數向量β的無偏估計：

(4)

當誤差項不滿足高斯-馬爾科夫條件時，使用最小二乘估計可能造成估計有偏。同時，最小二乘極易受異常值的影響，當樣本中存在異常值，使用最小二乘的誤差平方和可能產生較大偏差。而且，最小二乘確定的是自變量X取值x時，因變量Y的E(y|x)，無法給出因變量的一個全局分布。因此，Koenker等[9]提出分位數回歸，引入分位數回歸替換最小二乘估計，分位數回歸不對誤差分布做任何假定，可以得到一種穩健性的參數估計方法。分位數回歸使用最小化加權誤差絕對值和進行參數估計，對異常值的影響較小。同時，分位數回歸可以給出自變量X取值x時，因變量Y取值的各種τ分位數Qτ(y|x)下的回歸擬合結果，即相應Y取值的各種分位點的估計。

給定X=x，因變量Y的累積條件分布函數為FY(y|X=x)=P(Y≤y|X=x)。因變量Y的τth分位數的分位數回歸模型為：

Qτ(Y|X=x)=g(X,β)=xTβ(τ)

(5)

(6)

分位數回歸不能像最小二乘估計(4)一樣得到顯式解，通常基于損失函數ρτ(μ)，用數值算法對式(6)進行求解：

(7)

主要研究當τ=0.5時的LAD，LAD與分位數回歸一樣都具有對異常值不敏感，不對誤差分布做假設等特點，并且能得到穩健的估計結果，其他τ分位點同理。

主要研究對含有順序類別自變量的樣本進行建模，對于這類樣本建模時需要考慮偽分類的情況，避免模型因為存在偽分類發生過擬合現象。對于解決過擬合現象，可以通過正則化進行處理。Gertheiss等[20]曾提出對不含有連續型變量的模型，運用嶺回歸方法進行偽分類的識別：

(8)

順序類別自變量為有序類別，設置第一個類別為基準類，定義β1=0，這樣可以方便懲罰相鄰類別之間的系數差異。 將相鄰2類系數之差定義為Δβj=βj-βj-1。該方法實際上是使用l2懲罰的嶺回歸方法，但嶺回歸無法消除偏差，不能有效識別出偽分類，且沒有考慮連續型變量，因此模型不具有可解釋性。

之后，Huang等[22]提出TD-LASSO方法識別順序類別自變量中的偽分類，主要思想是運用TD對系數進行懲罰。由于TD-LASSO是基于最小二乘基礎上進行的，可知最小二乘對異常值敏感，且對模型誤差有正態假設條件。因此，根據LAD的不受異常值影響，不對模型誤差做假設還能得到穩健性估計結果的特點。用LAD替換最小二乘提出LAD-TD-LASSO。又由于LASSO采用無差別的懲罰系數，會導致部分自變量系數被過度壓縮降低估計效率，并影響選擇的一致性。因此，采用自適應LASSO對LAD-TD-LASSO方法加以改進提出LAD-TD-ADLASSO方法對偽分類進行識別，具體步驟如算法1所示。

算法1LAD-TD-ADLASSO算法

α1Z1,i+…+αkZk,i+εi

(9)

步驟2將步驟1中的參數進行合并整理：

ΔβmXm,i+α1Z1,i+…+αkZk,i+εi

(10)

α1Z1,i+…+αkZk,i+εi

(11)

步驟4使用分位數懲罰回歸識別偽分類，令φ=(α,Δβ2,…,Δβm,α1,…,αk)T，則

(12)

2 真實數據

在本節中，使用提出的LAD-TD-ADLASSO方法對2個實際數據例子進行分析研究，比較模型結果，觀察并分析該方法是否具有實用性和有效性。

2.1 實證分析1

分析一個來自醫學研究的真實數據集。數據來自新加坡眼科研究所眼科流行病學研究小組，見網站blog.uus.edu.sg/SEED。重點研究眼病、糖尿病、視網膜病變、青光眼、屈光不正和白內障等。選擇其中一個糖尿病數據集來展示所提出的LAD-TD-ADLASSO方法對具有順序類別自變量的模型的性能。一共3 280個數據，7個變量。數據存在缺失值，利用刪除法對數據進行處理，剩余3 242個數據。其中因變量為血壓，有4個順序類別自變量，分別是性別、受教育水平、家庭經濟狀況、糖尿病診斷。連續變量為年齡和BMI指標(身高與體重平方的比值)，根據以往的醫學研究可將年齡和BMI進行分類。年齡以10年為一分段將其劃分為一個5分類的順序類別自變量。BMI根據世界衛生組織(WHO)的劃分標準將其劃分成偏瘦、正常、肥胖前狀態、一級肥胖、二級肥胖、三級肥胖。整個糖尿病數據集有6個順序類別自變量。醫學研究中通常繪制table one向讀者展示研究人群的基本特征變量，見表1。特征變量中用(*)表示順序類別自變量。對順序類別自變量的分類情況做相應的統計描述，得到相應的頻數和頻率。連續型變量描述其均值、標準差、中位數和最大最小值。

表1 糖尿病數據集的統計描述

判斷順序類別自變量中是否存在偽分類，以受教育水平為例。繪制受教育水平影響血壓的階梯圖，如圖1所示。受教育水平中2、3分類對于血壓影響相似， 猜測受教育水平中2、3分類之間沒有顯著差異，可能存在偽分類。同樣，受教育水平中4、5分類對于血壓影響非常接近，猜測受教育水平中4、5分類可能存在偽分類。但不能單憑圖1就猜測哪些類別存在偽分類，需要通過提出的LAD-TD-ADLASSO方法進行偽分類的識別。

圖1 受教育水平影響血壓的階梯圖

首先，構造線性模型(模型1)為：

然后，通過LAD-TD-ADLASSO方法識別出順序類別自變量中隱藏的偽分類，檢測出受教育水平中3、4分類為偽分類，5、6分類為偽分類，將6分類的受教育水平進行偽分類識別與融合得到一個4分類的順序類別自變量。 家庭經濟狀況中2、3、4分類為偽分類，將其合并為一個2分類順序類別的自變量。BMI分類中4、5、6分類為偽分類，融合成一個4分類變量。建立模型2：

為了評價未做TD-ADLASSO偽分類識別的模型1和做過TD-ADLASSO偽分類識別的模型2的差異性。 同時，檢驗做LAD與做最小二乘下的模型的性能。 采用交叉驗證進行評價，交叉驗證可以用于評估模型的預測性能。比較平均絕對誤差(mean absolute error，MAE)和均方誤差(mean square error，MSE)指標，MAE和MSE越小，則表示模型效果越好：

表2 添加干擾后交叉驗證的MAE和MSE結果

將未做TD-ADLASSO偽分類識別的模型1基于最小二乘得到的結果記做OLS，模型1基于LAD得到的結果記做LAD。 將做過TD-ADLASSO偽分類識別的模型2基于最小二乘得到的結果記做TD+OLS，模型2基于LAD得到的結果記做TD+LAD。

根據表2可知，TD+OLS得到的MAE和MSE均小于OLS的MAE和MSE，TD+LAD下的MAE和MSE結果也均小于LAD的結果。說明經過TD-ADLASSO偽分類識別下的模型較未經過TD-ADLASSO偽分類識別的模型更具穩健性。同時，比較模型1和模型2內部的OLS和LAD結果，經過偽分類識別的模型2中，基于LAD得到的MAE和MSE均小于基于OLS得到的MAE和MSE。同樣，未經過TD-ADLASSO偽分類識別的模型1中，基于LAD得到的MAE和MSE均小于基于OLS得到的MAE和MSE。表明LAD較最小二乘更具穩健性。通過TD+LAD得到的MAE和MSE結果發現，使用了TD的模型MAE、MSE不僅更小，而且啞變量個數變少，模型更精簡，這符合解釋型建模的要求。而且使用了LAD的模型，在隨機干擾下表現更穩健。因此，提出的LAD-TD-ADLASSO方法在實際數據處理上具有可行性。

2.2 實證分析2

數據來自于成都市某中學不同年級學生的若干項目指標，該組數據包含了2 550位同學的大樣本數據，主要統計了2類指標：一類是身體素質，包括性別、身高、體重；另一類是體測成績，此類指標包括肺活量、50米跑、坐位體前屈、一分鐘跳繩。將50米跑作為因變量，考慮其他因素的影響，其中有6個順序類別的自變量，分別是年級、性別、體重、肺活量、坐位體前屈、一分鐘跳繩。 性別變量是二分類變量，1代表男性，2代表女性; 體重中的類別分為4個等級，等級1～4分別代表“正常”，“低體重”，“超重”，“肥胖”; 肺活量中的類別分為4個等級，等級1～4分別代表“不及格”，“及格”，“良好”，“優秀”; 坐位體前屈和一分鐘跳繩中的類別也分為4個等級，等級1～4代表的分類情況同上。變量及其變量解釋具體見表3。特征變量中用(*)表示順序類別自變量，其余的為連續型變量。對順序類別自變量的分類情況做相應的統計描述，連續型變量描述其均值、標準差、中位數和最大最小值。

表3 成都市某中學不同年級學生的若干項目指標的統計描述

判別順序類別自變量中是否存在偽分類，以體重水平為例，繪制體重的核密度估計圖，如圖2所示。由圖2可知，體重的第2、3、4分類非常近似，猜測體重的2、3、4分類存在偽分類。為了判斷是否確實存在偽分類，通過提出的LAD-TD-ADLASSO判別偽分類。

圖2 體重的核密度估計圖

首先，構造線性模型(模型3)為：

然后，通過LAD-TD-ADLASSO方法識別出順序類別自變量中隱藏的偽分類，如體重水平中的2、3、4分類為偽分類，將其合并為同一類，建立模型4：

為了評價未做TD-ADLASSO偽分類識別的模型3和做過TD-ADLASSO偽分類識別的模型4的差異性，以及檢驗做最小二乘和做LAD的模型的性能。同樣，在交叉驗證過程中，選擇訓練集殘差最大的前5%的數據，將其因變量Y添加總體因變量的2倍標準誤差得到其MAE和MSE結果，見表4。

表4 添加干擾后交叉驗證的MAE和MSE結果

根據表4可知，TD+LAD的MAE和MSE小于LAD的結果，且TD+OLS的MAE和MSE小于OLS的結果，表明通過TD-ADLASSO偽分類識別后的模型較未經過TD-ADLASSO偽分類識別模型更具穩健性，模型效果更好。比較模型3和模型4內部的OLS和LAD結果，經過TD-ADLASSO偽分類識別的模型4中，基于LAD得到的MAE和MSE均小于基于OLS得到的MAE和MSE。同樣，未經過TD-ADLASSO偽分類識別的模型3中，基于LAD得到的MAE和MSE均小于基于OLS得到的MAE和MSE，表明LAD較最小二乘更具穩健性。通過TD+LAD得到的MAE和MSE結果發現，使用了TD的模型MAE、MSE不僅更小，而且啞變量個數變少，模型更精簡，這符合解釋型建模的要求。而且使用了LAD的模型，在隨機干擾下表現更穩健。因此，提出的LAD-TD-ADLASSO方法在實際數據處理上具有可行性。

3 結論

針對含順序類別自變量的線性回歸模型，研究LAD和一種新型的順序類別自變量的偽分類識別及融合技術。對提出的偽分類識別技術TD-LASSO加以改進，利用自適應LASSO代替LASSO對不同參數進行不同壓縮提出TD-ADLASSO方法。之后，介紹了LAD的性質并與最小二乘作比較得出LAD的優勢，用LAD替換最小二乘做參數估計。對TD-ADLASSO偽分類識別技術加以改進得到LAD-TD-ADLASSO方法，該方法能有效地進行變量選擇與參數估計，得到穩健性結果。最后通過糖尿病數據集和成都市某中學不同年級學生的若干項目指標2個真實數據集，分別展示所提出的LAD-TD-ADLASSO方法識別偽分類的實用性，通過模型比較、交叉驗證得到的結果證實提出的LAD-TD-ADLASSO方法具有可行性。 然而，本文中僅考慮τ=0.5的LAD，當τ取任意分位數時，也是值得研究的內容。