LASSO方法在Cox回歸模型中的應(yīng)用*

閆麗娜 覃 婷 王 彤△

通過揭露癌癥死亡等終點(diǎn)事件發(fā)生的時(shí)間和基因表達(dá)譜數(shù)據(jù)之間的關(guān)系來研究不同患者的預(yù)后從而改進(jìn)治療策略，故而基于基因資料的生存分析越來越得到重視。生存資料的經(jīng)典方法是Cox比例風(fēng)險(xiǎn)回歸模型，該模型要求自變量之間相互獨(dú)立，且樣本量大于預(yù)測變量，但在基因表達(dá)譜資料中，預(yù)測變量遠(yuǎn)遠(yuǎn)大于樣本含量且各變量之間常具有強(qiáng)相關(guān)，呈現(xiàn)高維度和共線性，此時(shí)傳統(tǒng)Cox模型就不再適用。本文介紹的LASSO就是在系數(shù)的絕對值之和上增加一個(gè)約束條件來對高維資料進(jìn)行降維〔1〕，可得到更好的擬合效果。

基本思想與原理介紹

LASSO(the least absolute shrinkage and selection operator)由 Tibshirani〔2〕提出，由于它是對系數(shù)的絕對值而非系數(shù)的平方項(xiàng)進(jìn)行懲罰，也叫L1懲罰，它是在回歸系數(shù)的絕對值之和小于等于一個(gè)常數(shù)λ的約束條件下，使logL(β)達(dá)到最大來產(chǎn)生某些嚴(yán)格等于0的回歸系數(shù)，從而得到參數(shù)估計(jì)值。即:

調(diào)整參數(shù)λ的估計(jì)

該模型的復(fù)雜性主要在于確定最優(yōu)調(diào)整參數(shù)λ，調(diào)整參數(shù)λ的估計(jì)一般有三種方法:交叉驗(yàn)證法，廣義交叉驗(yàn)證法和無偏估計(jì)的風(fēng)險(xiǎn)分析。其中最常用的交叉驗(yàn)證法〔4〕由 Verweij和 van Houwelingen提出，之后他們又提出的留一法是交叉驗(yàn)證法的特例，要求K=n，方法是輪流將其中的一個(gè)研究對象作為驗(yàn)證集，剩下的n－1個(gè)研究對象作為訓(xùn)練集，用訓(xùn)練集擬合得到預(yù)測模型，把預(yù)測模型用到驗(yàn)證集中對該研究對象進(jìn)行預(yù)報(bào)評價(jià)，重復(fù)K(K=n)次。Cox模型中應(yīng)用交叉驗(yàn)證法決定最優(yōu)調(diào)整參數(shù)λ是以log L(β)為基礎(chǔ)的〔5〕，每一個(gè)研究對象對 log L(β)的貢獻(xiàn)為 l指所有研究對象在內(nèi)估計(jì)得到的log L(β)，l(－i)(f)是當(dāng)?shù)?i個(gè)對象作為驗(yàn)證集后估計(jì)的 log L(β)，i=1，2，3，…，n，^f(－i)(λ)是第i個(gè)研究對象作為驗(yàn)證集取出后剩下的數(shù)據(jù)資料中調(diào)整參數(shù)取λ時(shí)LASSO程序運(yùn)行得到的得分函數(shù)估計(jì)值。每一個(gè)對象i對log L(β)的貢獻(xiàn)和就是交叉驗(yàn)證 偏 似 然 值 CVL， CVL(λ) =，CVL取最大值時(shí)所對應(yīng)的值就是該模型的最優(yōu)調(diào)整參數(shù)〔5〕。

算 法

Tibshirani提出了LASSO程序的兩種計(jì)算方法，均以二次規(guī)劃為基礎(chǔ)，這兩種計(jì)算都是迭代的過程且涉及到重復(fù)的最小二乘求解問題，需要經(jīng)過p到2p的迭代，p是自變量的個(gè)數(shù)。為求解上式^β=arg max，就是把牛頓迭代(Newton-Raphson update)表達(dá)為迭代再加權(quán)最小二乘(IRLS)，然后用帶有限制條件的加權(quán)最小二乘程序取代了加權(quán)最小二乘程序。X表示變量的設(shè)計(jì)矩陣，且η =Xβ'，令 μ = ?l/?η，A= － ?2l/?ηηT，z= η +A－1μ。l(β)=log L(β)形式為:(z－ η)TA(z－ η)，由于 A 里面每一行或每一列的所有元素加起來為0，A很顯然是一個(gè)奇異矩陣，可以利用它的廣義逆，Tibshirani提出了用一個(gè)對角矩陣D代替設(shè)計(jì)矩陣A，兩矩陣含有相同的對角元素。在許多應(yīng)用中，n一般較小，且廣義逆的計(jì)算可以實(shí)現(xiàn)。LASSO實(shí)現(xiàn)程序包括以下步驟〔2〕:

(3)在限制條件∑|βj|≤λ下最大化(z－β'X)TA(z－ β'X)，估計(jì)出 β'。

(4)重復(fù)步驟(2)和步驟(3)，直到β'不再變化。

模擬研究

模擬生成高維、高度相關(guān)的微陣列基因數(shù)據(jù)，分別采用Cox比例風(fēng)險(xiǎn)回歸模型(逐步法)與L1懲罰Cox回歸模型對模擬數(shù)據(jù)進(jìn)行變量篩選與模型擬合。

1．協(xié)變量矩陣X的生成

真實(shí)的微陣列數(shù)據(jù)在一次實(shí)驗(yàn)時(shí)可以得到數(shù)千至上萬個(gè)基因〔6〕，為了計(jì)算的方便，我們僅模擬設(shè)置了1000個(gè)基因，并設(shè)樣本含量為100，保證協(xié)變量數(shù)量遠(yuǎn)遠(yuǎn)大于樣本含量，數(shù)據(jù)呈現(xiàn)高維性。100×1000的基因協(xié)變量矩陣中，每一行表示一條記錄，每一列表示一個(gè)基因，協(xié)變量矩陣服從均數(shù)為零的多元正態(tài)分布。將數(shù)據(jù)分成十塊等大小的基因塊，每塊包括100個(gè)基因變量，令它們的方差協(xié)方差矩陣的對角元素為1，非對角元素為0.8。10個(gè)基因塊對應(yīng)于基因表達(dá)的10個(gè)不同類別，不同類別的基因表達(dá)是獨(dú)立的，但是在同一個(gè)類別中的基因表達(dá)是兩兩相關(guān)的。

2．參數(shù)的設(shè)定

每個(gè)回歸系數(shù)對應(yīng)于它對應(yīng)變量的影響。在本文中，回歸系數(shù)參數(shù)的設(shè)定如下:當(dāng) 1≤j≤100，βj=0.01;當(dāng)101≤j≤200，βj從 －0.50 到 0.05;當(dāng) 201≤j≤1000，βj=0。表示在模擬基因矩陣數(shù)據(jù)中，只有少數(shù)協(xié)變量與應(yīng)變量有關(guān)，大部分協(xié)變量都是無關(guān)的。

3.生存時(shí)間變量的設(shè)定

首先生成服從(0，1)均勻分布的隨機(jī)數(shù) S，令S(t)=S，利用產(chǎn)生相應(yīng)的生存時(shí)間t，因此每個(gè)個(gè)體所對應(yīng)的生存時(shí)間為:，其中l(wèi)(·)表示對數(shù)似然函數(shù)。R2越大，則表示協(xié)變量所能解釋的那部分變異所占的百分比越大，模型擬合也越好。

模擬數(shù)據(jù)分析結(jié)果顯示，對于高緯度、強(qiáng)相關(guān)的基因模擬數(shù)據(jù)，采用逐步法進(jìn)行變量篩選，篩選出的自變量個(gè)數(shù)42大于LASSO選出的11，而模型評價(jià)指標(biāo)顯示逐步法Cox模型決定系數(shù)僅為0.3078，低于LASSO的0.6456，說明LASSO方法在將許多沒有意義的解釋變量壓縮為0之后，模型反而更優(yōu)，在Cox模型中進(jìn)行變量篩選用LASSO方法要比逐步篩選更具有競爭力。=1，2，…，p;xi=xi1，xi2，…，xip;i=1，2，…，n

4．刪失指示變量的設(shè)定

產(chǎn)生一列服從二項(xiàng)分布的隨機(jī)變量，發(fā)生1的概率為80%，即截尾比例為20%。

5．模擬數(shù)據(jù)分析及結(jié)果

按照以上步驟產(chǎn)生模擬微陣列數(shù)據(jù)，重復(fù)模擬50次，用逐步法擬合Cox比例風(fēng)險(xiǎn)回歸模型，變量入選標(biāo)準(zhǔn)為α=0.05，剔除標(biāo)準(zhǔn)為α=0.10。同時(shí)對每一數(shù)據(jù)集擬合基于LASSO的Cox回歸，調(diào)整參數(shù)的選擇采用交叉驗(yàn)證法，CV(λ)值最大時(shí)，即對應(yīng)最優(yōu)的調(diào)整參數(shù)λ。

模型評價(jià)比較采用Nagelkerke給出的一個(gè)可以用在刪失生存數(shù)據(jù)條件下的R2統(tǒng)計(jì)量，計(jì)算如下:R2=1

表1 模擬數(shù)據(jù)Cox逐步回歸與LASSO變量篩選個(gè)數(shù)表

表2 模擬數(shù)據(jù)Cox回歸與LASSO方法模型評價(jià)(R2)

實(shí)例分析

本實(shí)例來自于Van't Veer(2002)〔7〕等學(xué)者乳腺癌研究數(shù)據(jù)集，該數(shù)據(jù)集共包括259例乳腺癌患者，25000個(gè)微陣列基因數(shù)據(jù)。我們從中選擇沒有發(fā)生淋巴結(jié)轉(zhuǎn)移的乳腺癌患者78例，基因4751個(gè)。觀察事件的結(jié)局為乳腺癌是否發(fā)生遠(yuǎn)端轉(zhuǎn)移，其中44例沒有發(fā)生遠(yuǎn)端轉(zhuǎn)移，平均隨訪期為8.7年;34例在5年內(nèi)發(fā)生遠(yuǎn)端轉(zhuǎn)移，平均隨訪期為2.5年，截尾比例為56.4%。

上述實(shí)例資料顯示所研究變量個(gè)數(shù)4751遠(yuǎn)遠(yuǎn)大于樣本量78，存在高維度現(xiàn)象，提示不符合經(jīng)典Cox比例風(fēng)險(xiǎn)回歸模型的條件。

(1)首先采用SAS 9.2中PHREG語句，對該數(shù)據(jù)擬合Cox比例風(fēng)險(xiǎn)回歸模型(逐步法)，變量入選標(biāo)準(zhǔn)為α=0.05，剔除標(biāo)準(zhǔn)為α=0.10，擬合結(jié)果見表3。

表3 乳腺癌數(shù)據(jù)Cox逐步回歸模型變量篩選結(jié)果

表3結(jié)果顯示應(yīng)用逐步法進(jìn)行變量篩選，4751個(gè)基因中與乳腺癌發(fā)生遠(yuǎn)端轉(zhuǎn)移有關(guān)的基因有17個(gè)，且根據(jù)擬合模型的評價(jià)標(biāo)準(zhǔn)R2統(tǒng)計(jì)量R2=1－exp{－得出R2為0.1947。

(2)進(jìn)行基于LASSO的生存分析，調(diào)整參數(shù)λ的選擇采用交叉驗(yàn)證法得到圖1，2。

圖1和圖2分別為LASSO方法決定最優(yōu)調(diào)整參數(shù)λ和和篩選變量過程，結(jié)果顯示CVL(λ)取得最大值時(shí)為－167.8447，對應(yīng)的λ為5.95。在最優(yōu)調(diào)整參數(shù)λ為5.95時(shí)，LASSO篩選變量為13個(gè)，具體見表4。

圖1 調(diào)整參數(shù)λ與交叉驗(yàn)證CVL(λ)值變化圖

圖2 調(diào)整參數(shù)λ與LASSO篩選自變量變化圖

表4 乳腺癌數(shù)據(jù)LASSO變量篩選結(jié)果

表4結(jié)果顯示應(yīng)用LASSO進(jìn)行變量篩選，4751個(gè)基因中與乳腺癌發(fā)生遠(yuǎn)端轉(zhuǎn)移有意義的基因有13個(gè)，且根據(jù)擬合模型的評價(jià)標(biāo)準(zhǔn)R2統(tǒng)計(jì)量R2=1－得出R2為0.3923。

討 論

在腫瘤和其他疾病研究中，微陣列數(shù)據(jù)和其他的高通量檢測技術(shù)得到的數(shù)據(jù)正逐漸地用于診斷疾病的結(jié)果〔6〕。知道病人的病變轉(zhuǎn)移(或死亡)的風(fēng)險(xiǎn)信息對于成功地處理癌癥是很有必要的。因此如果能夠揭露死亡時(shí)間(或者其他終點(diǎn)事件的時(shí)間)和基因表達(dá)譜之間的關(guān)系就有可能得到更精確的診斷和改進(jìn)治療策略。本文介紹的LASSO方法是處理基因表達(dá)譜等高維數(shù)據(jù)生存分析的眾多方法中的一種〔8〕。通過Van't Veer等的乳腺癌數(shù)據(jù)，研究乳腺癌是否發(fā)生遠(yuǎn)端轉(zhuǎn)移與檢測到的4751個(gè)基因的關(guān)系。采用逐步法進(jìn)行變量篩選篩出有意義的自變量個(gè)數(shù)17大于LASSO篩出的13，而模型評價(jià)指標(biāo)Cox模型的決定系數(shù)R2僅為0.1947，低于LASSO的0.3923，LASSO模型優(yōu)于Cox模型，說明LASSO方法在將一些沒有意義或者意義很小的解釋變量系數(shù)壓縮為0之后，模型反而更優(yōu)。無疑證明LASSO模型是通過將一些無意義或者意義很小的自變量的系數(shù)壓縮為0而對高維數(shù)據(jù)進(jìn)行降維，而得到的一個(gè)更為穩(wěn)定科學(xué)且容易解釋的模型，適合于基因數(shù)據(jù)的生存資料分析。

1．Tibshirani RJ．Regression shrinkage and selection via the lasso．Journal of the Royal Statistical Society，1996，58:267-288．

2．Tibshirani RJ．The Lasso method for variable selection in the Cox model．Statistics in Medicine，1997:385-395．

3．Gui J，Li H．Penalized Cox regression analysis in the high dimensional and low-sample size settings with applications to microarray gene expression data．Bioinformatics，2005:3001-3008．

4．Verweij PJ．Cross-validation in survival analysis．Statistics in Medicine，1993，12:2305-2314．

5．Van HC，Bruinsma T，Van't Veer LJ，et al．Cross-validated Cox regression on microarray gene expression data．Statistics in Medicine，2006，25:3201-3216．

6．Segal MR，Dahlquist KD，Conklin BR．Regression approaches for microarray data analysis．Journal of Computational Biology，2003，10:961-980．

7．van de Vijver MJ，He YD，van't Veer LJ，et al．A gene-expression signature as a predictor of survival in breast cancer．N Engl J Med，2002，347:1999-2009．

8．Tim H，Nam HC，Lukas M，et al．Least angle and ?1penalized regression．Statistics Surveys，2008:61-93．