基于正則化回歸的組學數據變量篩選方法*

哈爾濱醫科大學衛生統計學教研室(150081)

謝宏宇　侯　艷　李　康△

?

基于正則化回歸的組學數據變量篩選方法*

哈爾濱醫科大學衛生統計學教研室(150081)

謝宏宇侯艷李康△

近年來，隨著各種生物檢測技術的發展，醫學研究中出現了各種高通量數據，如基因組、蛋白質組和代謝組學數據等，變量選擇是生物標志物識別和建立分類模型的重要環節，由于高維組學數據中的絕大多數變量對分類并不起作用，并且存在多重共線性、模型過擬合等問題，傳統的基于最小二乘方法估計的線性回歸并不適用于高維數據[1]。在高維組學數據特征變量篩選過程中，由于變量數目很多，子集篩選方法計算量巨大，并可能由于選擇不同的篩選準則導致篩選的結果有很大差異[2]；維數縮減方法雖然能夠避免計算量大的問題，但是由于模型中的變量不再是原始變量，模型的可解釋性差；而正則化回歸方法由于在解回歸方程的過程中可以同時實現參數估計和變量篩選，且計算速度快，對變量數目沒有限制，因此受到研究者的關注[3]。這類方法不僅能夠用于單一組學數據的變量篩選，同時也能拓展到多組學數據融合的情況，因此在實際中具有很好的應用前景。本文將對正則化回歸方法及在高維組學數據中的應用做一綜述。

正則化的基本原理

正則化是指在原有的損失函數的基礎上增加懲罰回歸系數的正則項。記β=(β1,β2,…,βm)為回歸系數，m為總自變量的個數，則在線性回歸中，通過最小化損失函數可以得出對應模型的回歸系數估計值

(1)

其中‖·‖2是L2范數，表示向量各元素平方和的平方根，該式表示取右端函數最小值的模型參數，這實際是傳統的最小二乘估計。當變量個數較多時，利用該式估計得出的回歸模型存在過擬合的風險，正則化則可以在保留所有特征變量的情況下，避免過擬合的發生，其基本原理是通過增加的正則項，減少所有特征變量回歸系數估計值的數量級，具體的表現形式如下：

(2)

其中，λP(β)表示正則化項，λ為正則化參數，P(β)為回歸系數的懲罰函數，主要目的是用來平衡模型對樣本數據的擬合程度以及回歸模型的預測能力。在正則化項中，如果正則化參數設定較大會使得每個回歸系數估計值偏小；如果回歸系數估計值小到一定程度時，相當于因變量只等于常數項，類似于擬合了一條水平直線，導致欠擬合，產生過高的偏差。如果模型中涉及到高階項，則回歸系數的估計值越小，對應的曲線越光滑，從而使函數得到簡化，實際中需要選擇合適的正則化參數值。目前正則化參數的選擇可以通過偏差原理、Engl誤差極小原理、Hansen 的L曲線準則、擬最優準則和交叉驗證等方法進行確定[4]。回歸模型的系數估計值可以通過梯度下降等方法進行求解。

基于正則化回歸的單一組學數據變量篩選方法

1.嶺回歸

嶺回歸(ridge regression)方法由Hoerl 和Kennard提出[5]，其基本思想是在傳統最小化殘差平方和基礎上加入回歸系數的L2范數懲罰項從而收縮回歸系數。最小化回歸系數的L2范數，會使稀疏矩陣中每個元素的值都很小，但并不一定為0。回歸系數估計值的表達式如下：

(3)

其中λ是正則化參數。由于L2范數可以收縮回歸系數估計值，因此能夠在一定程度上避免模型的過擬合。嶺回歸的主要特點是通過L2范數對回歸系數的連續收縮，能夠使每個變量的系數變小，從而通過損失無偏性提高了模型的預測能力。主要缺點是，嶺回歸將所有的預測變量均保留在模型中，因此在分析高維組學數據時會導致模型的可解釋性較差。

2.lasso回歸

Tibshirani于1996年提出了基于線性回歸的最小化的絕對收縮和選擇算子(least absolute shrinkage and selection operator,lasso)來收縮回歸系數，這種方法在損失函數中增加了回歸系數的L1范數懲罰項，表示為‖·‖1，代表向量中各個元素絕對值之和。在回歸系數的絕對值之和小于一個常數的約束條件下，使殘差平方和最小化，能夠使部分回歸系數等于0，同時實現回歸系數收縮和變量篩選，從而提高了模型的可解釋性[6]。回歸系數估計值的表達式如下：

(4)

隨著正則化參數λ的增大，lasso方法能夠不斷地縮小回歸系數的估計值，使其趨近于0，實現回歸系數的稀疏化。在高維組學數據中，最常用于估計lasso回歸系數的方法為最小角算法(least angle regression，LARS)[7]，這種算法相對于最小二乘回歸能夠很好地解決lasso回歸的計算問題。lasso回歸存在一定的局限性，即在自變量個數m遠大于樣本量n時，只能保證lasso回歸中最多選擇n個變量；同時，如果一組變量高度相關時，這種算法只傾向于選擇其中之一，而不關心選擇的究竟是哪個變量[5]。

3.自適應lasso回歸

Zou(2006)發現lasso回歸中L1范數懲罰項對所有回歸系數懲罰強度相同，從而導致了回歸系數估計值不具有漸進正態性[8]。另外，lasso回歸只有在兩種特定條件下的變量篩選才具有相合性[5]。為了解決這個問題，Zou(2006)提出了一種新的lasso回歸方法，即自適應lasso回歸，回歸系數估計值的表達式為

(5)

其中wj代表回歸系數的權重。該方法的特點是可以對不同的系數設置不同的懲罰權重，從而改進其估計值的準確性。如果依賴于數據本身對權重做出恰當的選擇，則自適應lasso回歸具有相合性和漸進正態性，并且能夠避免局部最優化的問題[9]。Breheny在2013年將自適應lasso回歸方法用于微陣列數據進行變量篩選[10]，并采用預測誤差均方作為評價指標與lasso回歸進行比較，結果表明自適應lasso回歸方法的預測誤差均方小于lasso回歸，并且篩選出的差異變量更少[11]。

4.樸素彈性網和彈性網算法

高維組學數據往往具有高度相關性和分組特征(例如來自于同一通路的基因)，如前所述lasso回歸方法針對以上兩種情況進行變量篩選效果不理想。因此Zou于2003年提出了彈性網算法(elastic net)[12]，該方法既能夠同時實現變量的自動篩選和回歸系數的連續收縮，又能夠保證選擇出同一分組內與因變量相關性大的變量。樸素彈性網(naive elastic net)算法是最基礎的彈性網算法，主要是將lasso回歸的懲罰項和嶺懲罰項相結合，其表達式為

(6)

其中，λ1和λ2均為非負的正則化參數。若記α=λ2/(λ1+λ2)，上式等價于

(7)

并且(1-α)‖β‖1+α‖β‖2≤t

其中，t為一個常數界值，(1-α)‖β‖1+α‖β‖2稱為彈性網懲罰，α∈[0,1]。當α=0時，該式為lasso回歸懲罰項；當α=1時，為嶺回歸懲罰項。

樸素彈性網算法的參數估計分為兩個階段：首先固定λ2找到嶺回歸系數，然后通過λ1進行系數壓縮。雖然樸素彈性網算法能夠克服傳統lasso回歸的部分不足，但模擬實驗表明只有當它接近嶺回歸或者lasso回歸時，才能獲得較理想的變量篩選結果。因此，Zou于2005年又提出了對樸素彈性網系數進行重縮放，這種方法即為目前的彈性網算法。當對預測變量進行標準化后，彈性網方法的回歸系數與樸素彈性網的回歸系數之間具有如下關系：

(8)

其中1+λ2為收縮因子。

研究結果表明，彈性網算法與lasso回歸和嶺回歸相比具有較好的篩選變量的性能，L1范數可以實現自動變量篩選，L2范數可以實現連續收縮，尤其在自變量之間存在較強的相關性時，彈性網算法能夠明顯的提高預測的準確性[3]。Zou將幾種變量篩選方法應用于白血病患者的基因表達數據，目的是篩選用于診斷和預測白血病分型的基因。結果表明，彈性網算法構建的模型分類效果優于支持向量機和懲罰logistic回歸等方法，并且能對組內基因進行篩選。由于彈性網算法估計出的系數不具有相合性和漸進正態性[5,8]，因此Zou于2009年提出將自適應加權L1懲罰納入到彈性網算法正則項中提高估計準確性，即自適應彈性網算法(adaptive elastic-net)，該方法可以將其視為彈性網算法和自適應lasso的結合，具有相合性和漸進正態性[13]。

正則化參數λ決定了模型中回歸系數估計值的大小和稀疏化的程度。確定正則化參數的基本方法有交叉驗證[14]、貝葉斯信息準則(BIC)[15]、Cp統計量[7]和赤池信息量準則(AIC)[16]。Zou于2007年從變量篩選的角度證明了BIC相對于其他方法更適用于參數值的選擇，該方法能夠產生一個更加稀疏的模型[12]。Chen等認為使用BIC準則在高維數據中篩選變量的標準具有一定的任意性，因此提出了擴展的BIC方法(EBIC)，這種方法既考慮了未知參數的個數，也考慮了模型空間的復雜性，并且能夠更加嚴格地控制差異變量錯誤發現率[17]。

5.分組lasso回歸

(9)

其中L為變量的組別數，l=1,2,…,L，X(l)代表組l中與因變量有關的X列的子矩陣，β(l)是組的系數向量，pl是第l組中包含的變量個數。這種方法利用了‖β(l)‖2在β(l)=0處不可微的性質，將該組從模型中剔除。其主要思想是篩選出對因變量有影響的特征組，同時通過選擇合適的參數λ調整組別個數，λ值越大，對各分組作用的懲罰越大，則模型中保留的組數越少。雖然分組lasso回歸能夠實現對組別的篩選，但是只能篩選出模型中整個組內的變量回歸系數β(l)=0的特征組，這一缺點限制了其應用。當每個特征組內只包含一個自變量時，則該方法即為傳統的lasso方法。

6.稀疏組lasso回歸

實際研究中不僅僅需要實現組別的稀疏化，同時還需要實現組內變量的稀疏化，例如，研究者識別感興趣基因通路，同時對該條基因通路中的關鍵基因進行篩選。因此Simon(2010)提出將lasso回歸和分組lasso回歸相結合，引進了稀疏組lasso(sparse-group lasso)回歸[19]，其表達式為

(10)

稀疏組lasso回歸的方法與彈性網方法相似，不同點在于該種方法是利用在懲罰回歸系數為0時不可微的性質，將稀疏為0的組別從模型中去除，實現組間稀疏化，而彈性網方法則保留了所有的分組。Simon(2013)將lasso回歸、分組lasso回歸和稀疏組lasso回歸的方法應用于乳腺癌患者的基因表達數據中，并比較了三種方法的篩選效果。結果表明稀疏組lasso回歸的變量篩選性能優于lasso回歸和分組lasso回歸：稀疏組lasso回歸的分類正確率達到70%，而分組lasso回歸和lasso回歸的分類正確率分別為60%和53%。由于在癌癥數據中添加分組的信息對于分類非常有意義，同時分組信息可以幫助更加深入的了解生物學機制，因此對于癌癥數據的分析，稀疏組lasso回歸的方法有很大的優勢[19]。

多組學融合數據變量篩選方法

傳統的變量篩選方法一般均可應用于單一組學數據的變量篩選，如基于回歸、基于機器學習和基于網絡的方法等，目前應用于多組學融合數據變量篩選方法相對較少，而基于正則化的變量篩選方法可以實現多組學數據的融合和變量篩選。

1.稀疏廣義典型相關分析

典型相關分析(canonical correlation analysis，CCA)是用于研究兩組變量之間關系的常用方法。Tenenhaus于2011年提出了正則化的廣義典型相關分析(regularized generalized canoncial analysis，RGCCA)方法，該方法可用于分析三個或者更多的變量集合間的關系[20]。RGCCA是一種基于主成分分析的方法，用于研究多個數據集中變量之間的關系。RGCCA成分的性質和解釋性受每組變量之間有用性和相關性的影響。RGCCA主要基于使多個數據樣本中生成新的綜合變量的相關程度最大化的思想進行求解。

實際中，在每組變量中識別出在組間關系中起顯著作用的變量子集非常重要，因此Tenenhaus在2014年提出了稀疏的廣義典型相關分析(sparse generalized canonical correlation analysis，SGCCA)，這種方法通過對外部權重向量加上L1懲罰，在同一方法中結合了RGCCA和L1懲罰項[20]。Tenenhaus將SGCCA方法應用到兒科神經膠質瘤數據，結果表明：與RGCCA比較，SGCCA方法能夠篩選出在組間相關作用中具有更小差異的變量組合[21]。

2.稀疏偏最小二乘回歸

偏最小二乘(PLS)的變量篩選方法已經成功地應用于代謝組學數據。其主要原理是分別在自變量和因變量中提取出成分，使各自提取出的成分盡可能多的解釋各數據的變異信息，同時使提取成分的相關程度達到最大。

Cao(2008)在此基礎上提出了稀疏偏最小二乘回歸(sparse least squares regression,SPLS)的變量篩選方法，該方法能夠同時實現數據的整合和變量篩選，其主要思想是在PLS的基礎上，通過Q2值作為評價指標對構建模型的成分數量進行選擇，同時對每個成分加上lasso懲罰，實現變量篩選。研究表明：該方法應用于高維數據集分析時，與PLS相比具有更高的穩定性，能夠更好地進行變量篩選[22]。

3.結構組稀疏算法

多尺度數據分析的關鍵性問題是，數據結構異質性整合和特征變量篩選的穩定性。基于結構組稀疏算法(structure grouping sparsity，SGS)的多尺度數據變量篩選方法的目的是根據實際數據給出一個可以解釋和預測的模型。其主要思想是，根據實際數據建立應變量Y與自變量X=(X1,X2,…,Xm)關系的廣義線性模型(可拓展至非線性)，實現對不同來源的異質性數據在不同水平上進行組間和組內特征篩選。其表達式為：

(11)

展　　望

基于正則化回歸的變量篩選方法，克服了傳統變量篩選方法的不足，且隨著研究的深入，需要不斷的更新和發展。該方法的發展一直圍繞著擬合較好的模型應該具有預測準確度高、模型的可解釋性強的特點；本著模型本身具有優良的參數估計性質，即無偏性、有效性、相合性和漸進正態性。基于正則化的變量篩選方法不僅能夠應用于單一組學數據的變量篩選，也能夠應用于多組學數據的融合和變量篩選。然而，這種方法的懲罰項選擇及其統計性質，以及參數求解等問題都有待進一步研究。展望未來，高維組學的數據研究將實現跨組學的超高維變量篩選，從而更全面的研究疾病的發生機制，因此這類方法將會具有較好的前景。

[1]趙奕林，朱真峰，周清雷.適用于大規模高維多類別數據分類的并行非線性最小二乘分類器.小型微型計算機系統，2014,3:579-583.

[2]Daniel PB,Pierluigi C.Introduction to the theory of complexity.Prentice Hall.ISBN 0-13-915380-2,1994.

[3]Zou H,Hastie T.Regularization and variable selection via the elastic net.J.R.Statist.Soc.B,2005,67:301-320.

[4]閔濤，葛寧國，黃娟，等.正則參數求解的微分進化算法.應用數學與計算數學學報,2010,24(2):23-27.

[5]Hoerl A,Kennar R.Ridge regression.In Encyclopedia of Statistical Sciences,1998,8:129-136.

[6]Tibshirani R.Regression shrinkage and selection via the lasso.J.R.Statist.Soc.B,1996,58:267-288.

[7]Efron B,Hastie T,Johnstone I,et al.Least angle regression.Ann.Statist.,2004,32:407-499.

[8]Zou H.The adaptive lasso and its oracle properties.Journal of the American Statistical Association,2006,101:1418-1429.

[9]Li ZT,Mikko J,Sillanpaa.Overview of lasso-related penalized regression methods for quantitative trait mapping and genomic selection.Theor Appl Genet,2012,125:419-435.

[10]Scheetz T,Kim K,et al.Regulation of gene expression in the mammalian eye and its revevance to eye disease.Proc.Natl.Acad.Sci,2006,103:14429-14434.

[11]Patrick B,Jian H.Group descent algorithms for nonconvex penalized linear and logistic regression models with grouped predictors.Stat Comput,2015,25:173-187.

[12]Zou H,Hastie T.Regression shrinkage and selection via the elastic net,with application to microarrays,2003，1-26.

[13]Zou H,Zhang H.On the adaptive elastic-net with a diverging number of parameters.Ann.Statist.,2009,37:1733-1751.

[14]Hastie,Tibshirani R,Friedman JH.The elements of statistical learning.Springer,New York,2009.

[15]Zou H,Hastie T,Tibshirani R.On the “degrees of freedom” of the lasso.Ann Stat,2007,35:2173-2192.

[16]Akaike H.New look at the statistical model identification.IEEE T Autom Contr,1974,19:716-723.

[17]Chen J,Chen Z.Extended Bayesian information criteria for model selection with large model spaces.Biometrika,2008,95:759-771.

[18]Yuan M,Lin Y.Model selection and estimation in regression with grouped variables.Journal of the Royal Statistical Society,Series B,2007,68(1):49-67.

[19]Simon N,Friedman J,et al.A Sparse-Group lasso.Journal of computational and Graphical Statistics,2013,22:231-245.

[20]Tenenhaus A,Tenenhuas M.Regularized generalized canonical analysis.Psychometrika,2011,76:257-284.

[21]Tenenhaus.Variable selection for generalized canonical correlation analysis.Biostatistics,2014:1-15.

[22]Le Cao KA,Rossouw D,et al.A sparse PLS for variable selection when integrating omics data.Stat Appl Genet Mol Biol,2008,7(1)：1-32.

(責任編輯：郭海強)

國家自然科學基金資助(81573256，81473072)；黑龍江省博士后資助經費(LBH-Z14174)

李康,E-mail:likang@ems.hrbmu.edu.cn