孫紅衛(wèi)　楊文越　王　慧　羅文?！『藢殹⊥酢⊥?/p>

?

懲罰logistic 回歸用于高維變量選擇的模擬評價*

孫紅衛(wèi)1,2楊文越2王慧1羅文海2胡乃寶2王彤1△

【提要】目的logistic回歸是生物醫(yī)學研究中常用的方法，可以進行影響因素篩選、概率預測、分類等。高通量測序技術得到的數據給高維變量選擇問題帶來挑戰(zhàn)。懲罰logistic回歸可以對高維數據進行變量選擇和系數估計，且其有效的算法保證了計算的可行性。方法本文介紹了常用的懲罰logistic算法如LASSO(least absolutes shrinkage and selection operator)、EN(elastic net)、SCAD(smoothly clipped absolute deviation)、MCP(minimax concave penalty)以及SIS(sure independence screening)等，并用模擬數據對各方法進行評價。結果(1)各方法的結果與自變量間的相關程度有關，不同懲罰logistic回歸的精確性與自變量間的相關程度有關，如果相關較高，LASSO或EN的結果較好，而在相關較低時，MCP或SCAD結果較好；(2)結合SIS的方法傾向于少選變量，誤選率低，但敏感度也低，而LASSO、MCP、SCAD選擇變量較多，誤選率高，但敏感度較高;(3)當自變量間低度相關時，SIS的三種方法結果非常接近，但相關較高時，SIS+LASSO的結果表現較好。結論采用非小細胞型肺癌的基因數據集進行實例分析，并表明如何根據模擬實驗的結論，在多種方法的不同結果間進行選擇。

高維變量選擇懲罰logistic回歸LASSOMCPSCADSIS

logistic回歸模型已經被廣泛應用在生物醫(yī)學領域，它適用于響應變量為分類資料的情況。它通過對概率進行l(wèi)ogit變換，并對其與協變量的線性組合建立模型，用來探索影響因素或者預測疾病的發(fā)生概率。

隨著高通量技術的快速發(fā)展，現在的遺傳學研究已經提供了豐富的數據集，用來識別與疾病(如癌癥、自身免疫性疾病、心臟病和精神疾病等)有關的遺傳變異[1-3]。這些數據共同的特點是變量維數遠遠大于樣本量，所以傳統統計方法難以進行模型選擇和參數估計；同時存在著計算成本大、最優(yōu)化難以實現等問題[4]。

對高維數據的研究，通常預計只有部分基因位點或SNP與疾病有關，即滿足稀疏性。懲罰回歸方法用來解決高維變量選擇問題，比如LASSO(least absolutes shrinkage and selection operator)[5]，SCAD(smoothly clipped absolute deviation)[6]，MCP(minimax concave penalty)[7]。很多學者將其擴展到logistic模型中。Shevade 和Keerthi[8]為L1規(guī)則化logistic回歸提出了一種有效的算法，并將其用于基因篩選研究。Jiang 等[9]用Coordinate Descent算法求解了MCP 懲罰的logistic 回歸。Jiang和Zhang[10]給出了GPLUS算法來計算SCAD懲罰和MCP懲罰的logistic估計。

本文對目前的高維懲罰logistic方法進行介紹，并用數值模擬的方法來評價各種高維logistic模型的優(yōu)劣。運用一個非小細胞型肺癌的實際案例，篩選出可能影響其五年生存率的基因位點，為下一步的研究作參考。

logistic回歸模型的介紹

設響應變量Y為二分類資料，yi～B(1,πi)，i=1,2,…,n，即共有n個觀測，影響πi的有p個自變量x1,x2,…,xP。

懲罰logistic回歸方法

1.橋回歸、嶺回歸、LASSO懲罰模型、EN模型

Frake和Friedman[11](1993)提出了橋回歸?？芍苯訉蚧貧w推廣到logistic回歸模型。橋回歸logistic模型的目標函數：-(g(X,β),Y)+λ‖β‖q。

橋回歸又稱Lq懲罰模型，參數估計式中‖β‖q為懲罰項，λ為調整參數，常常用AIC、BIC或者交叉驗證等準則確定最合適的λ值。當λ=0時，即不對原有的模型做懲罰，隨著λ的增大，原模型的懲罰力度會越來越大，則被選入的變量越來越少。當0

其中除嶺回歸在零點不存在奇異性，LASSO、EN和末稀疏模型的懲罰項在零點均存在奇異性，正是由于這三種方法的懲罰項在零點的導數不存在，所以將接近于零的系數壓縮為零(壓縮范圍由調整參數λ決定)，從而實現模型的稀疏性假定。

圖1　橋回歸不同懲罰下的估計系數

Fan和Li(2001)[6]提出評判模型選擇優(yōu)劣的標準，即Oracle的三個性質:(1)稀疏性：模型選擇中對參數的估計應自動實現系數的稀疏性，將一些不重要的變量的系數壓縮為零。(2)無偏性：估計的參數值應該是無偏的或者近似無偏的。(3)連續(xù)性：參數估計與對應的數據應該是連續(xù)的，從而避免模型擬合的偏差與預測的不穩(wěn)定性。

從圖1可知，當q=0.1時即末稀疏模型，類似于最優(yōu)子集選擇(q=0)，系數較小時壓縮為0，滿足稀疏性，系數較大時不進行壓縮，滿足無偏性，但不滿足oracle性質中的連續(xù)性，即在不連續(xù)點的估計系數會不穩(wěn)定；當q=1時即LASSO模型不滿足Oracle性質中的無偏性，但滿足稀疏性和連續(xù)性；當q=2即嶺回歸，不滿足Oracle中的無偏性和稀疏性，滿足連續(xù)性；EN模型介于LASSO和嶺回歸之間，滿足稀疏性，但是同樣的λ下壓縮為0的范圍小于LASSO，同樣不滿足Oracle性質中的無偏性。

2.SCAD懲罰模型和MCP懲罰模型

Fan 和Li[6](2001)提出了非凹懲罰似然的方法SCAD(smoothly clipped absolute deviation)來選擇變量，Zhang[7](2007)提出了MCP(minimax concave penalty)方法，MCP方法和SCAD方法類似，具備SCAD的優(yōu)點，滿足Oracle的三個性質。

SCAD懲罰函數：

MCP懲罰函數：

MCP懲罰用于logistic回歸的目標函數：

如圖2所示，SCAD和MCP這兩種方法，在自變量系數較小時，都被壓縮為零；自變量系數很大時，不進行壓縮，這些自變量系數是無偏的；介于二者之間時，都進行部分壓縮，以保證連續(xù)性。SCAD和MCP滿足Oracle的三個性質。MCP是在所有滿足無偏性條件的懲罰函數中，擁有最小最大凸性，有很好的理論性質，比如不需要很強的不可代表條件(LASSO所需要的)，就可證明以很高的概率能夠正確選擇變量[7]。

圖2　SCAD 和MCP的估計系數

3.SIS方法

Fan and Lv(2007)[13]提出了SIS(sure independence screening)方法，通常用于超高維的情況下，它可以快速地從超高維降到高維。該方法基本原理如下：

令W={1≤i≤p∶βi≠0}為真實的稀疏模型，令d=(d1,…,dp)T；通過如下方式得到p維回歸變量：d=XTy。對任意給定的0<γ<1建立子模型：

Wγ={1≤i≤p∶di為前|γn|相關性最大的子集}

在子模型中我們選取與響應變量有關聯的|γn|個自變量，由此一來就可以將模型從p維降到|γn|維。在將超高維數據利用SIS降到高維后，可以利用上述方法如LASSO，MCP，SCAD等方法對模型進行懲罰回歸。

模擬實驗設計

1.模擬設計

在利用R軟件實現各種情況的模擬實驗設計后，采用glmnet軟件包實現LASSO，EN，嶺回歸；利用ncvreg軟件包實現MCP，SCAD。利用SIS軟件包實現SIS+LASSO，SIS+MCP，SIS+SCAD。懲罰參數λ的選取采用交叉驗證。

2.評價指標FDR、PSR、RMSPE、RMSE

FDR(false discovery rate)和PSR(positive select rate)自1995年被Benjamini和Hochberg提出以來，被廣泛研究，特別是在高維數據回歸建模和復雜數據的多重比較領域有很好的應用[14]。

其中FP(falsepositive)表示假陽性的個數，即在真實模型中的系數是零，但是被估計成非零。TP(truepositive)表示真陽性的個數，即在真實模型中的系數是非零，估計的結果也是非零的個數。m為真實模型中非零系數的個數。FDR的意義是估計為非零的系數中假陽性占的比例。PSR的意義是真實模型中非零系數中真陽性占的比例。一般而言，FDR越接近于0，PSR越接近于1，說明該選擇方法越好。

模擬實驗結果

1.模擬實驗1p=150，n=100，自變量間最大相關系數r=0.2，見表1。

表1　p=150且自變量間低度相關時五種方法的比較結果

當自變量最大、相關系數較小時，MCP在FDR，PSR，RMPSE和RMSE四個方面均表現較好，即誤選率最低，預測準確性最好，估計系數準確性最高，選取了67.41%的真實變量。嶺回歸由于沒有進行變量選擇，所以FDR接近于1。EN和LASSO的表現相似，FDR要大于MCP，而PSR卻低于MCP。SCAD預測誤差和參數估計誤差都較低，但是FDR較高。

2.模擬實驗2p=150，n=100，自變量間最大相關系數r=0.8，見表2。

表2　p=150且自變量間高度相關時五種方法的比較結果

當自變量最大相關系數較大時，EN和LASSO在FDR、PSR和RMSE均表現均較好，誤選率較低，且能選出60%以上的真實變量，估計系數的準確性也較好。MCP和SCAD在RMSPE方面表現優(yōu)于其他三種方法，即預測準確性方面較好。綜合四個指標來看，EN和LASSO方法表現最好。

3.模擬實驗3n=200，p=1000，自變量間最大相關系數r=0.2，見表3。

表3　p=1000且自變量間低度相關時六種方法的比較結果

當自變量間最大相關系數較小時，SIS+MCP、SIS+SCAD和SIS+LASSO 三種方法與LASSO、MCP和SCAD三種方法相比較，其FDR比較小，即誤選率比較低，PSR比較低，即選出的真實變量比較少。SIS+MCP、SIS+SCAD和SIS+LASSO傾向于少選變量，而LASSO、MCP、SCAD則傾向于多選變量。其中MCP的FDR較低，SCAD其次，LASSO的FDR較高，RMSPE指標和RMSE指標上,MCP和SCAD的表現好于LASSO。所以在自變量間相關較小時，如果傾向于誤選率低，那么采用SIS+MCP和SIS+SCAD，如果傾向于盡可能多地選出真實變量，且預測的準確性較高，估計系數的準確性較好，則應該選擇SCAD和MCP。

4.模擬實驗4n=200，p=1000，自變量間最大相關系數r=0.8，見表4。

當自變量間最大相關系數較大時，SIS的三個方法FDR較低，其中SIS+LASSO的FDR最低，且PSR最高，表現較好，LASSO的表現要好于MCP和SCAD。SIS+LASSO與LASSO相比，SIS+LASSO的FDR較低，而LASSO的PSR較高，這是由于SIS的方法傾向于少選變量，而LASSO的方法傾向與多選變量。所以當自變量間最大相關系數較大時，如果傾向于少選變量，誤選率低，可以用SIS+LASSO，如果傾向于盡可能多地選出真實變量則應選擇LASSO。

表4　p=1000且自變量間高度相關時六種方法的比較結果

綜上，自變量最大相關系數較小時，可結合SIS+MCP或者SIS+SCAD與MCP或者SCAD相結合分析；自變量最大相關系數較大時可結合LASSO和SIS+LASSO分析。

實例分析

本文采用Sandy D.Der et al[16]關于非小細胞型肺癌(non-small cell lung cancer，NSCLC)的數據，其中包含181例I期和II期病人，共50248個基因位點的表達水平數據，以病人五年生存率作為響應變量。

首先求得該數據集前5000個基因之間的相關矩陣，并求出每個變量與其他變量的最大相關系數，相關系數的描述見表5和圖5。

表5　肺癌數據5000個基因表達之間最大相關系數的統計描述

圖3　肺癌基因數據5000基因表達之間最大相關系數的分布

從最大相關系數的分布來看，基因表達之間的最大相關系數較高，以0.57為中心對稱分布，最大相關系數大于0.5的比例達到68.26%，而最大相關系數大于0.6的比例達到38.26%。由于維數p=50248，LASSO、MCP和SCAD等無法處理超高維數據，而SIS的方法是先降維后，再用LASSO、MCP和SCAD等方法進行變量選擇，所以本實例分析中采用SIS+LASSO、SIS+MCP和SIS+SCAD三種方法分別對肺癌基因數據進行分析，計算其RMSPE以及AUC(area under curve)來衡量它們的預測結果。

表6　SIS的三種方法預測肺癌基因數據的RMSPE和AUC

由表6可以看到，SIS+LASSO的預測誤差稍高。SIS的三種方法的預測效果很好，AUC均達到了0.8以上，SIS+MCP和SIS+SCAD達到了0.83以上。

但由模擬實驗的結果看出，SIS+LASSO的預測誤差通常稍大，但是變量選擇的FDR和PSR卻在自變量相關較大時表現較好，所以不能僅通過預測性能的高低來判斷變量選擇的準確性。

表7　不同方法篩選的基因個數和基因列表

從SIS三種方法篩選出的基因列表(表7)可以看到，三種方法都篩選出了8個基因位點，且其中5個基因位點都是相同的，分別是探針集205433_at，217045_x_at，226476_s_at，234490_at，240078_at。這些探針集對應的基因很可能與Ⅰ、Ⅱ期非小細胞型肺癌的五年生存率有關。有2個探針集228748_at，236391_at是其中兩種方法共同篩選出的基因。

從對肺癌數據的相關性分析知，基因之間存在較大的相關性，所以通過模擬實驗的結果知，SIS+LASSO的FDR較低，PSR較高，這樣SIS+LASSO篩選出的基因位點很可能與I、II期非小細胞型肺癌的五年生存率有關。

這些方法篩選出的基因位點，可以作為候選的基因，但這些結果僅僅說明這些候選基因與五年生存率存在相關，不一定存在因果關系，所以需要進一步從專業(yè)意義上進行驗證。

討　　論

近年來，生物信息數據挖掘技術被廣泛應用，生物醫(yī)學中常見的響應變量是分類資料，比如通過病例對照研究或隊列研究，來探索與疾病相關的基因位點或SNP等。但是，高維變量選擇方法用于logistic回歸方面的研究較少。本文對用于高維數據變量選擇的懲罰logistic方法進行了比較，得出了以下結論：

不同懲罰logistic回歸的表現與自變量間的相關程度有關，如果相關較高，LASSO或EN的結果較好，而在相關較低時，MCP或SCAD結果較好，MCP的估計誤差較低，SCAD的預測誤差較低。結合SIS的方法傾向于少選變量，誤選率低，但敏感度也低，而LASSO、MCP、SCAD選擇變量較多，誤選率高，但敏感度較高。當自變量間低度相關時，SIS的三種方法結果非常接近，但相關較高時，SIS+LASSO的結果表現較好。

就像文獻[17]所揭示的，LASSO懲罰(q=1)通過對相關的變量比較平等對待，而不像最優(yōu)子集選擇(q=0)，一旦一個很強的變量進入了，它就會排除其他相關的變量。MCP和SCAD等懲罰介于LASSO懲罰和最優(yōu)子集選擇之間，所以相對于LASSO，會傾向于排除其他相關的變量。當自變量相關較高的時候，特別是相關較高的變量都是重要變量，LASSO會傾向于全選出來，而MCP等方法，當一個很強的變量被選入后，其他相關的重要變量就很難進入模型。所以當自變量相關較高時，且相關的變量都是重要變量時，LASSO的結果會較好。在實際數據分析時，要根據自變量間的相關性，選取性能較好的方法來進行變量選擇或概率預測。

本文的結論可以用來對高維logistic變量選擇方法的選用提供參考，但仍然存在一些不足。例如懲罰參數的選取只采用交叉驗證，沒有考慮采用AIC、BIC、EBIC等方法來調整參數，這有待后面的研究進一步解決。

[1]Cai T,Lin X,Carroll RJ.Identifying genetic marker sets associated with phenotypes via an efficient adaptive score test.Biostatistics,2012,13(4):776-790.

[2]Zhu J,Hastie T.Classification of gene microarrays by penalized logistic regression.Biostatistics,2004,5(3):427-443.

[3]Risch N,Merikangas K.The future of genetic studies of complex human diseases.Science,1996,273(5281):1516-1517.

[4]李根,鄒國華,張新雨.高維模型選擇方法綜述.數理統計與管理,2012,31(4):641-651.

[5]Tibshirani R.Regression shrinkage and selection via the lasso.Journal of the Royal Statistical Society.Series B(Methodological),1996,58(1):267-288.

[6]Fan J,Li R.Variable selection via non concave penalized likelihood and its oracle properties.Journal of the American Statistical Association,2001,96(4):1348-1360.

[7]Zhang CH.Penalized linear unbiased selection.Technical Report,2007,52(17):374-393.

[8]Shevade SK,Keerthi SS.A simple and efficient algorithm for gene selection using sparse logistic regression.Bioinformatics,19(17):2246-2253.

[9]Jiang D,Huang J,Zhang Y.The cross-validated AUC for MCP-logistic regression with high-dimensional data.Stat Methods Med Res,2013,22(5):505-518

[10]Jiang W,Zhange CH.Path following algorithm for penalized logistic regression using SCAD and MCP.Communications in Statistics-Simulation and Computation,2014,43(5):1064-1077.

[11]Frank I,Friedman J.A statistical view of some chemometrics regression tools.Technometrics,1993,35:109-148.

[12]Zou H,Hastie T.Regularization and variable selection via the elastic net.Journal of the Royal Statistical Society:Series B(Statistical Methodology),2005,67(2):301-320.

[13]Fan J,Lv J.Sure independence screening for ultra-high dimensional feature Space.Journal of the Royal Statistical Society Series B,2007,70(4):849-911.

[14]李瑞.SNP定位的一種降維及變量選擇方法.合肥:中國科技大學,2011.

[15]Alfons A,Croux C,Gelper S.Sparse least trimmed squares regression for analyzing high-dimensional large data sets.The Annals of Applied Statistics,2013,7(1):226-248.

[16]Der SD,Sykes J,Pintilie M,et al.Validation of a histology-independent prognostic gene signature for early-stage,non-small-cell lung cancer including stage IA patients.J Thorac Oncol,2014,9(1):59-64.

[17]Mazumder R，Friedman J,Hastie T.SparseNet:Coordinate Descent With Nonconvex Penalties.Journal of the American Statistical Association,2011:1125-1138.

[18]閆麗娜,覃婷,王彤.LASSO方法在Cox回歸模型中的應用.中國衛(wèi)生統計,2012,29(1):58-60,64.

[19]覃婷,閆麗娜,王彤.基于腫瘤患者高維生物信息的生存預測.中國衛(wèi)生統計,2011,28(1):101-103,105.

[20]張秀秀,王慧,田雙雙,等.高維數據回歸分析中基于LASSO的自變量選擇.中國衛(wèi)生統計,2013,30(6):922-926.

(責任編輯：鄧妍)

國家自然科學基金資助(81473073)，國家自然科學基金青年基金(81502891)

王彤，E-mail:wtstat@21cn.com

1.山西醫(yī)科大學衛(wèi)生統計學教研室(030001)

2.濱州醫(yī)學院衛(wèi)生統計學教研室