基于 LASSO 回歸篩選變量以構(gòu)建預(yù)測(cè)模型
在構(gòu)建疾病風(fēng)險(xiǎn)評(píng)估和預(yù)測(cè)模型時(shí),我們會(huì)想要納入盡可能多的指標(biāo)去構(gòu)建模型。這么做雖然會(huì)更準(zhǔn)確,但卻不實(shí)用。因?yàn)獒t(yī)生不可能讓所有剛?cè)朐旱幕颊甙阉械臋z查都做一遍。我們更希望在盡可能知道較少指標(biāo)的情況下而不損失風(fēng)險(xiǎn)評(píng)估準(zhǔn)確性,這樣才可以實(shí)現(xiàn)效益和效率最大化。又或者,在大量基因中選取可能對(duì)疾病有重要影響的基因(比如我們前面在看完還不會(huì)來揍我 | 差異分析三巨頭 —— DESeq2、edgeR 和 limma 包 | 附完整代碼 + 注釋中得到的差異基因),用于構(gòu)建風(fēng)險(xiǎn)預(yù)測(cè)模型。這個(gè)時(shí)候,我們就可以使用 LASSO 回歸在眾多變量中進(jìn)行篩選,獲取重要變量,也就是特征,用于構(gòu)建模型啦!今天我們就來分享一下具體要怎么做嘍!
我們先介紹一下Lasso回歸的原理,大家如果不需要的話,也可以跳過這部分直接去看用法和代碼部分喲!
Lasso回歸(Least Absolute Shrinkage and Selection Operator)是一種用于線性回歸問題的正則化方法,它可以用來降低模型的復(fù)雜性,防止過擬合,并選擇重要的特征變量。
LASSO 回歸的原理
Lasso回歸的原理基于在線性回歸(想深入了解線性回歸的小伙伴們,詳見跟我一起啃西瓜書系列中的第三章線性模型3.2小節(jié))的基礎(chǔ)上添加一個(gè)L1正則化項(xiàng)(不要慌!后面會(huì)講!),以約束模型的參數(shù)估計(jì)。L1正則化項(xiàng)是參數(shù)向量中絕對(duì)值之和,它的存在迫使許多參數(shù)變?yōu)榱?/strong>,從而實(shí)現(xiàn)了變量選擇(特征篩選)的功能。Lasso回歸的優(yōu)化目標(biāo)是最小化以下?lián)p失函數(shù):
其中:
- 是第 個(gè)觀測(cè)值的目標(biāo)值。
- 是第 個(gè)觀測(cè)值的第 個(gè)特征變量的值。
- 是第 個(gè)特征變量的系數(shù),當(dāng)其為 0 時(shí), 就被踢出方程咯,這樣就達(dá)到篩選變量的效果啦!
- 是正則化參數(shù),用于控制正則化的強(qiáng)度。較大的 值將導(dǎo)致更多的參數(shù)估計(jì)為零。
簡(jiǎn)單介紹一下正則化。
正則化(Regularization)是一種在機(jī)器學(xué)習(xí)和統(tǒng)計(jì)建模中常用的技術(shù),旨在幫助模型更好地泛化到未見過的數(shù)據(jù),避免過擬合問題。過擬合(overfitting)是指模型在訓(xùn)練數(shù)據(jù)上表現(xiàn)得非常好,但在未見過的測(cè)試數(shù)據(jù)上表現(xiàn)較差的情況。簡(jiǎn)單來說,模型過度“記憶”了訓(xùn)練數(shù)據(jù)的細(xì)節(jié)和噪聲,導(dǎo)致在新數(shù)據(jù)上的預(yù)測(cè)表現(xiàn)不佳。過擬合通常出現(xiàn)在模型過于復(fù)雜或者訓(xùn)練數(shù)據(jù)較少的情況下。就好像是一個(gè)學(xué)生過度死記硬背了大量的題目答案,但實(shí)際上并沒有理解問題的本質(zhì)。因此,當(dāng)遇到新的類似問題時(shí),雖然能夠準(zhǔn)確回答原有題目,但對(duì)于新問題的推廣能力較差。
正則化就好像是在機(jī)器學(xué)習(xí)中的模型訓(xùn)練過程中施加的一種規(guī)則,它會(huì)對(duì)模型進(jìn)行“限制”或“懲罰”,以確保它不會(huì)變得過于復(fù)雜。
想象一下,你正在嘗試通過擬合一條曲線來擬合一組散點(diǎn)數(shù)據(jù),但你擔(dān)心曲線會(huì)過于彎曲,導(dǎo)致在散點(diǎn)之間出現(xiàn)怪異的起伏。正則化就像是在你的曲線擬合問題中引入一個(gè)規(guī)則,告訴你不要讓曲線變得太彎曲。嶺回歸(L2正則化)和Lasso回歸(L1正則化)就是為了解決線性回歸中的過擬合問題而出現(xiàn)的,彈性網(wǎng)絡(luò)是L1和L2正則化的結(jié)合,可以平衡兩者的影響。
- L1正則化(也稱為L(zhǎng)asso正則化):它告訴你,如果曲線中有一些不太重要的部分,就將它們收縮到零,從而削減模型的復(fù)雜度。這相當(dāng)于在模型中刪除不必要的特征或參數(shù),使模型更簡(jiǎn)單。
- L2正則化(也稱為嶺正則化):它告訴你,如果曲線的某些部分變得過于陡峭,就要降低它們的陡峭程度。這使得曲線的參數(shù)不能太大,從而控制了模型的復(fù)雜度(在本文的后半部分,L2正則化馬上就要登場(chǎng)啦)。
正則化的目標(biāo)是在模型的損失函數(shù)中添加一個(gè)額外的項(xiàng),這個(gè)項(xiàng)與模型參數(shù)相關(guān),通過這個(gè)項(xiàng)來限制參數(shù)的大小或稀疏性,從而防止模型在訓(xùn)練數(shù)據(jù)上過度擬合。這種“限制”或“懲罰”有助于模型更好地泛化到新數(shù)據(jù),降低了過擬合的風(fēng)險(xiǎn)。
LASSO 回歸與嶺回歸、彈性網(wǎng)絡(luò)的區(qū)別
Lasso回歸(Least Absolute Shrinkage and Selection Operator):Lasso回歸引入了L1正則化項(xiàng),通過約束參數(shù)向量的絕對(duì)值之和來實(shí)現(xiàn)稀疏性。這使得Lasso能夠執(zhí)行特征選擇,將許多參數(shù)估計(jì)為零,從而簡(jiǎn)化模型。Lasso回歸傾向于在高維數(shù)據(jù)中選擇少數(shù)關(guān)鍵特征,提高了模型的可解釋性。
Lasso回歸的優(yōu)化目標(biāo)如下:
嶺回歸(Ridge Regression):與Lasso不同,嶺回歸使用L2正則化項(xiàng),其目標(biāo)是最小化參數(shù)的平方和,而不是絕對(duì)值和。這意味著嶺回歸不會(huì)將參數(shù)估計(jì)為零,但會(huì)減小參數(shù)的絕對(duì)值,從而減小參數(shù)的大小。因此嶺回歸適合在存在多個(gè)相關(guān)特征的情況下穩(wěn)定模型,但不執(zhí)行特征選擇。
嶺回歸的優(yōu)化目標(biāo)如下:
彈性網(wǎng)絡(luò)(Elastic Net):彈性網(wǎng)絡(luò)是L1和L2正則化的組合,綜合了嶺回歸和Lasso的優(yōu)點(diǎn)。它的損失函數(shù)包括L1和L2正則化項(xiàng),可以同時(shí)控制參數(shù)的大小和稀疏性,在特征選擇和參數(shù)縮減之間進(jìn)行權(quán)衡。彈性網(wǎng)絡(luò)適用于數(shù)據(jù)集中既存在高度相關(guān)特征又需要進(jìn)行特征選擇的情況。
彈性網(wǎng)絡(luò)的優(yōu)化目標(biāo)如下:
為何選擇 LASSO 回歸進(jìn)行特征篩選
選擇在什么時(shí)候使用嶺回歸、Lasso回歸或彈性網(wǎng)絡(luò)來篩選變量,其實(shí)取決于具體情況和數(shù)據(jù)性質(zhì),有一些指導(dǎo)原則大家可以了解一下。
- 當(dāng)你懷疑有多個(gè)特征對(duì)目標(biāo)變量有影響,而且這些特征可能存在高度相關(guān)性時(shí),可以考慮使用嶺回歸:嶺回歸通過L2正則化可以減小參數(shù)的絕對(duì)值,但不會(huì)將參數(shù)估計(jì)為零,因此適合在存在多個(gè)相關(guān)特征的情況下穩(wěn)定模型。它不會(huì)執(zhí)行嚴(yán)格的特征選擇,而是縮小了參數(shù)的幅度。
- 當(dāng)你懷疑只有少數(shù)幾個(gè)特征對(duì)目標(biāo)變量有顯著影響,并且希望執(zhí)行特征選擇時(shí),可以考慮使用Lasso回歸:Lasso回歸通過L1正則化傾向于將許多參數(shù)估計(jì)為零,從而實(shí)現(xiàn)了特征選擇的功能。如果你想要簡(jiǎn)化模型并提高可解釋性,Lasso通常是一個(gè)不錯(cuò)的選擇。
- 當(dāng)你既想要特征選擇,又想要保留某些相關(guān)特征時(shí),可以考慮使用彈性網(wǎng)絡(luò):彈性網(wǎng)絡(luò)結(jié)合了L1和L2正則化,允許你在特征選擇和參數(shù)縮減之間進(jìn)行權(quán)衡。這在存在多個(gè)相關(guān)特征,但你仍然希望進(jìn)行一定程度的特征選擇時(shí)非常有用。通過調(diào)整兩個(gè)正則化參數(shù)的值,你可以控制模型的行為。
總的來說,選擇哪種正則化方法應(yīng)該基于以下考慮:
- 數(shù)據(jù)性質(zhì):了解數(shù)據(jù)的特征和相關(guān)性,以確定哪種正則化方法更適合。
- 目標(biāo):確定你是否更關(guān)注特征選擇、參數(shù)縮減或平衡兩者。
- 交叉驗(yàn)證:通常需要使用交叉驗(yàn)證來選擇適當(dāng)?shù)恼齽t化參數(shù),以確保模型的性能。
通過上面的介紹,大家應(yīng)該對(duì)我們選擇Lasso進(jìn)行特征篩選有了一定了解。一般我們?cè)谂R床預(yù)測(cè)模型構(gòu)建過程中,都會(huì)期望找到相對(duì)較少的特征數(shù)量,希望用更少的特征就可以對(duì)目標(biāo)變量產(chǎn)生較大影響,且希望模型更具可解釋性。在臨床預(yù)測(cè)中,可能存在大量的生物醫(yī)學(xué)特征,其中一些可能不是很重要或相關(guān)。Lasso可以自動(dòng)將某些特征的系數(shù)估計(jì)為零,從而實(shí)現(xiàn)特征選擇,幫助減少模型中的不必要的特征,這樣生成的模型更容易解釋,只包含最重要的特征,更便于醫(yī)生或臨床決策者理解模型的預(yù)測(cè)結(jié)果。而且Lasso在樣本數(shù)量相對(duì)較少的情況下也可以工作良好,尤其是在高維數(shù)據(jù)中,它有助于防止過擬合,即使數(shù)據(jù)點(diǎn)有限(試問大家能拿到多少有效樣本哈哈哈哈哈哈哈嗚嗚嗚嗚嗚嗚嗚嗚嗚嗚,幫你們哭一下)。所以,我們見到的臨床模型多數(shù)是基于Lasso回歸構(gòu)建的!
當(dāng)然!咱也不能無腦Lasso,還是要多嘗試滴!萬一其他模型更合適呢!
這個(gè)時(shí)候,可能會(huì)有小伙伴有疑問,那既然彈性網(wǎng)絡(luò)是L1和L2正則化的組合,看起來似乎它考慮更全面吶,我們?yōu)槭裁床恢苯佑盟鼇磉M(jìn)行特征篩選呢?
其實(shí),當(dāng)數(shù)據(jù)集中存在高度相關(guān)的特征時(shí),Lasso有可能隨機(jī)選擇其中的一個(gè)特征并將其他特征的系數(shù)設(shè)為零,這可能導(dǎo)致結(jié)果不穩(wěn)定,因?yàn)榫唧w選擇哪個(gè)特征取決于模型的隨機(jī)性。這在高維數(shù)據(jù)中尤其成問題。而彈性網(wǎng)絡(luò)可以平衡這兩種正則化的影響,既能夠進(jìn)行特征選擇,又能夠處理多重共線性問題。它不會(huì)像Lasso那樣過于隨機(jī)地選擇特征,所以更加穩(wěn)定。
但是吧,彈性網(wǎng)絡(luò)也不一定在所有數(shù)據(jù)中都表現(xiàn)最好。當(dāng)數(shù)據(jù)量較大、質(zhì)量較好,且問題適用于線性模型時(shí),彈性網(wǎng)絡(luò)可能表現(xiàn)良好。當(dāng)數(shù)據(jù)量較小、數(shù)據(jù)質(zhì)量較差,或問題具有高度非線性性質(zhì)時(shí),彈性網(wǎng)絡(luò)的效果可能就會(huì)表現(xiàn)不如其他模型。
在實(shí)際應(yīng)用中,我們還是要嘗試不同的模型(包括彈性網(wǎng)絡(luò)、嶺回歸、Lasso等)并使用交叉驗(yàn)證來確定哪種模型在給定問題上表現(xiàn)最佳。
LASSO 回歸在預(yù)測(cè)模型中的用法
用法一
這也是最常見的用法,采用Lasso進(jìn)行回歸,根據(jù)十折交叉驗(yàn)證,篩選得到候選預(yù)測(cè)因子。然后用候選預(yù)測(cè)因子繼續(xù)做后續(xù)的多因素Logsitic回歸或者多因素Cox回歸,得到最終的臨床預(yù)測(cè)模型,用于后續(xù)的Nomogram等分析。
比如:Mo R, Shi R, Hu Y, Hu F. Nomogram-Based Prediction of the Risk of Diabetic Retinopathy: A Retrospective Study. J Diabetes Res. 2020 Jun 7;2020:7261047.
這篇文章采用Lasso回歸,對(duì)19個(gè)預(yù)測(cè)因子進(jìn)行篩選,最終選得7個(gè)預(yù)測(cè)因子。然后對(duì)7個(gè)預(yù)測(cè)因子構(gòu)建多因素Logistic回歸模型。
然后作者就基于這7個(gè)因素構(gòu)建列線圖(Nomogram),并進(jìn)行了后續(xù)3個(gè)度的評(píng)價(jià)。
這些圖的繪制,我們后續(xù)都會(huì)陸續(xù)更新喲!
用法二
當(dāng)我們研究的因素較多,可以先進(jìn)行單因素Logistic或Cox回歸,先篩選一批可能的預(yù)測(cè)因子;然后再采用Lasso回歸進(jìn)行篩選,將篩選得到的預(yù)測(cè)因子,再次進(jìn)行多因素Logistic或Cox回歸,確定最終模型。
比如:Liu J, Huang X, Yang W, Li C, Li Z, Zhang C, Chen S, Wu G, Xie W, Wei C, Tian C, Huang L, Jeen F, Mo X, Tang W. Nomogram for predicting overall survival in stage II-III colorectal cancer. Cancer Med. 2020 Apr;9(7):2363-2371.
這篇文章通過單變量分析發(fā)現(xiàn)59個(gè)因素,基于專業(yè)背景知識(shí)又加了5個(gè) P > 0.05 的因素(有的時(shí)候我們建模,發(fā)現(xiàn)某個(gè)專業(yè)背景上有意義的變量,卻沒有進(jìn)入多因素分析階段,就可以通過這樣的表述把它加進(jìn)去),共64個(gè)因素。然后對(duì)這個(gè)64因素進(jìn)行Lasso篩選,發(fā)現(xiàn)了6個(gè)系數(shù)不為零,也就是有意義的預(yù)測(cè)因子。
作者對(duì)這6個(gè)預(yù)測(cè)因子,做了多因素Cox回歸,最終發(fā)現(xiàn)4個(gè)因子,并使用這4個(gè)因子構(gòu)建最終的最優(yōu)模型。
然后構(gòu)建了4個(gè)因子的Nomo圖以及后續(xù)的3個(gè)度的驗(yàn)證等等!
用法三
這種用法相對(duì)少見,就是直接用Lasso-Logit或者Lasso-Cox進(jìn)行變量篩選,篩選后,用最小誤差解構(gòu)建模型,無需再對(duì)Lasso篩選得到的多個(gè)因子,再進(jìn)行多因素的Logistic或Cox進(jìn)行分析。
LASSO 回歸存在的問題
盡管Lasso回歸應(yīng)用廣泛且優(yōu)點(diǎn)多多,但不可避免的,它也存在一些問題和局限性:
- 參數(shù)選擇問題:選擇合適的正則化參數(shù) 是一個(gè)挑戰(zhàn),需要進(jìn)行交叉驗(yàn)證或其他模型選擇方法。
- 估計(jì)偏差:Lasso傾向于對(duì)參數(shù)進(jìn)行嚴(yán)格的稀疏化,可能導(dǎo)致估計(jì)偏差,尤其是在特征之間存在高度相關(guān)性時(shí)。
- 難以處理大規(guī)模數(shù)據(jù)集:對(duì)于大規(guī)模數(shù)據(jù)集,Lasso的計(jì)算成本可能很高,因?yàn)橐紤]所有特征的組合。
LASSO 回歸代碼實(shí)現(xiàn)及結(jié)果解讀
數(shù)據(jù)用到的是我們之前看完還不會(huì)來揍我 | 差異分析三巨頭 —— DESeq2、edgeR 和 limma 包 | 附完整代碼 + 注釋中得到的差異基因表達(dá)矩陣,進(jìn)行Lasso回歸需要樣本分類信息,我們使用生存數(shù)據(jù),大家可以通過 UCSC Xena(http://xena./)進(jìn)行下載。從 UCSC Xena 下載數(shù)據(jù)的方法詳見:看完還不會(huì)來揍/找我 | TCGA 與 GTEx 數(shù)據(jù)庫(kù)聯(lián)合分析 | 附完整代碼 + 注釋的數(shù)據(jù)準(zhǔn)備部分。
如果大家想快快實(shí)戰(zhàn),不想耗費(fèi)時(shí)間去運(yùn)行前面的差異分析步驟得到差異分析表達(dá)矩陣或者下載生存數(shù)據(jù),完全沒問題!我已經(jīng)把用到的數(shù)據(jù)上傳到了GitHub,大家可以在公眾號(hào)后臺(tái)回復(fù)lasso回歸
,即可得到存放這些數(shù)據(jù)的鏈接。
咱們正式開始!
代碼看起來有點(diǎn)長(zhǎng),但其實(shí)正式的Lasso回歸部分一點(diǎn)都不長(zhǎng)!前面大多是數(shù)據(jù)處理部分,大家如果有現(xiàn)成的數(shù)據(jù),直接從”現(xiàn)在正式開始lasso回歸!”這里開始看就好啦!
################################# Lasso 回歸 ###################################
# 加載包
# 如果大家在運(yùn)行的時(shí)候發(fā)現(xiàn)某個(gè)函數(shù)不存在或者不能用,可能是忘記加載包啦
# 我這里如果沒有寫全的話,先給大家說聲抱歉,因?yàn)橛袝r(shí)候我的環(huán)境里可能已經(jīng)加載了就直接用了
# 所以可能會(huì)忘記寫某個(gè)包,大家不要慌,如果遇到了且不知道函數(shù)來自哪個(gè)包的時(shí)候
# 咱們?nèi)グ俣裙雀栲枥锱纠膊橐幌戮秃茫〔慌?!你是最棒的?br>library(glmnet)
library(survival)
library(dplyr)
### 加載數(shù)據(jù)
# 咱就今天使用的數(shù)據(jù)是BRCA數(shù)據(jù)集中的癌癥樣本,如果大家想了解這個(gè)數(shù)據(jù)是如何得出的,
# 可以查看:https://mp.weixin.qq.com/s/0DVzPXrNCQgG7kTtMw74jQ
brca_tumor <- readRDS('./lasso_data/brca_tumor.rds')
head(brca_tumor)[1:5, 1:4]
# TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A TCGA-BH-A201-01A TCGA-E2-A14T-01A
# TSPAN6 8.787903 12.064743 11.801304 10.723661
# TNMD 0.000000 2.807355 4.954196 6.658211
# DPM1 11.054604 11.292897 11.314017 11.214926
# SCYL3 10.246741 9.905387 11.117643 12.093748
# C1orf112 8.965784 10.053926 9.957102 9.503826
# deseq2得到的差異基因
load('./lasso_data/DEG_deseq2.Rdata')
logFC = 2.5
P.Value = 0.01
k1 <- (DEG_deseq2$pvalue < P.Value) & (DEG_deseq2$log2FoldChange < -logFC)
k2 <- (DEG_deseq2$pvalue < P.Value) & (DEG_deseq2$log2FoldChange > logFC)
DEG_deseq2 <- mutate(DEG_deseq2, change = ifelse(k1, 'down', ifelse(k2, 'up', 'stable')))
table(DEG_deseq2$change)
# down stable up
# 693 43491 1339
# edgeR得到的差異基因
load('./lasso_data/DEG_edgeR.Rdata')
logFC = 2.5
P.Value = 0.01
k1 <- (DEG_edgeR$PValue < P.Value) & (DEG_edgeR$logFC < -logFC)
k2 <- (DEG_edgeR$PValue < P.Value) & (DEG_edgeR$logFC > logFC)
DEG_edgeR <- mutate(DEG_edgeR, change = ifelse(k1, 'down', ifelse(k2, 'up', 'stable')))
table(DEG_edgeR$change)
# down stable up
# 634 30531 1904
# limma得到的差異基因
load('./lasso_data/DEG_limma_voom.Rdata')
logFC = 2.5
P.Value = 0.01
k1 <- (DEG_limma_voom$P.Value < P.Value) & (DEG_limma_voom$logFC < -logFC)
k2 <- (DEG_limma_voom$P.Value < P.Value) & (DEG_limma_voom$logFC > logFC)
DEG_limma_voom <- mutate(DEG_limma_voom, change = ifelse(k1, 'down', ifelse(k2, 'up', 'stable')))
table(DEG_limma_voom$change)
# down stable up
# 404 17753 248
# 定義函數(shù)挑選差異基因
deg_filter <- function(df){
rownames(df)[df$change != 'stable']
}
# 取交集
all_degs <- intersect(intersect(deg_filter(DEG_deseq2), deg_filter(DEG_edgeR)), deg_filter(DEG_limma_voom))
length(all_degs)
# [1] 442
# 三個(gè)包的結(jié)果取交集得到442個(gè)差異基因
# 獲取取交集得到的差異基因的表達(dá)矩陣
exp_brca_final <- exp_brca %>% filter(rownames(exp_brca) %in% all_degs)
dim(exp_brca_final)
# [1] 442 1217
head(exp_brca_final)[1:5, 1:4]
# TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A TCGA-BH-A201-01A TCGA-E2-A14T-01A
# KLHL13 6.044394 7.734710 10.702173 7.321928
# HSPB6 6.189825 8.312883 9.709084 10.121534
# PDK4 9.306062 11.454299 11.350387 11.959640
# MEOX1 3.700440 7.599913 8.455327 8.199672
# SCN4A 3.169925 4.459432 7.515700 6.643856
# 加載生存數(shù)據(jù)
sur_brac <- read.table('./lasso_data/TCGA-BRCA.survival.tsv', header = T)
head(sur_brac)
# sample OS _PATIENT OS.time
# 1: TCGA-C8-A275-01A 0 TCGA-C8-A275 1
# 2: TCGA-BH-A1F8-11B 1 TCGA-BH-A1F8 1
# 3: TCGA-BH-A1F8-01A 1 TCGA-BH-A1F8 1
# 4: TCGA-AC-A7VC-01A 0 TCGA-AC-A7VC 1
# 5: TCGA-AN-A0AM-01A 0 TCGA-AN-A0AM 5
# 6: TCGA-C8-A1HJ-01A 0 TCGA-C8-A1HJ 5
# 我們?cè)倏匆幌卤磉_(dá)數(shù)據(jù)
head(exp_brca_final)[1:5, 1:4]
# TCGA-E9-A1NI-01A TCGA-A1-A0SP-01A TCGA-BH-A201-01A TCGA-E2-A14T-01A
# KLHL13 6.044394 7.734710 10.702173 7.321928
# HSPB6 6.189825 8.312883 9.709084 10.121534
# PDK4 9.306062 11.454299 11.350387 11.959640
# MEOX1 3.700440 7.599913 8.455327 8.199672
# SCN4A 3.169925 4.459432 7.515700 6.643856
# 咱們需要轉(zhuǎn)置一下,然后合并!
# 整合表達(dá)數(shù)據(jù)和生存數(shù)據(jù)
exp_brca_final <- as.data.frame(t(exp_brca_final))
exp_brca_final$sample <- rownames(exp_brca_final)
lasso_data <- left_join(sur_brac, exp_brca_final, by = 'sample')
head(lasso_data)[1:5, 1:8]
# sample OS _PATIENT OS.time KLHL13 HSPB6 PDK4 MEOX1
# 1: TCGA-C8-A275-01A 0 TCGA-C8-A275 1 6.714246 6.409391 10.108524 8.339850
# 2: TCGA-BH-A1F8-11B 1 TCGA-BH-A1F8 1 11.071462 12.134105 13.323336 10.804131
# 3: TCGA-BH-A1F8-01A 1 TCGA-BH-A1F8 1 7.483816 7.228819 9.620220 5.459432
# 4: TCGA-AC-A7VC-01A 0 TCGA-AC-A7VC 1 7.098032 7.523562 9.552669 5.930737
# 5: TCGA-AN-A0AM-01A 0 TCGA-AN-A0AM 5 7.330917 6.491853 7.209453 4.392317
# 咱們整理一下下,樣本名設(shè)為行名,剔除沒用的列
rownames(lasso_data) <- lasso_data$sample
lasso_data <- lasso_data[ , -c(1, 3)]
lasso_data <- na.omit(lasso_data)
head(lasso_data)[1:5, 1:8]
# OS OS.time KLHL13 HSPB6 PDK4 MEOX1 SCN4A E2F2
# TCGA-C8-A275-01A 0 1 6.714246 6.409391 10.108524 8.339850 5.044394 10.727920
# TCGA-BH-A1F8-11B 1 1 11.071462 12.134105 13.323336 10.804131 9.079485 7.011227
# TCGA-BH-A1F8-01A 1 1 7.483816 7.228819 9.620220 5.459432 5.169925 8.791163
# TCGA-AC-A7VC-01A 0 1 7.098032 7.523562 9.552669 5.930737 4.169925 9.202124
# TCGA-AN-A0AM-01A 0 5 7.330917 6.491853 7.209453 4.392317 3.584963 10.099348
dim(lasso_data)
# [1] 1194 444
# 現(xiàn)在的數(shù)據(jù)有1194行(樣本)、444列(基因),第一列為OS(樣本結(jié)局,1為死亡,0為生存)
# 第二列為OS.time(生存時(shí)間,單位為day),后面的442列為基因的表達(dá)量。
# 現(xiàn)在正式開始lasso回歸!
# 把生存時(shí)間單位從day轉(zhuǎn)換為year
lasso_data$OS.time <- lasso_data$OS.time / 365
head(lasso_data)[1:5, 1:8]
# OS OS.time KLHL13 HSPB6 PDK4 MEOX1 SCN4A E2F2
# TCGA-C8-A275-01A 0 0.002739726 6.714246 6.409391 10.108524 8.339850 5.044394 10.727920
# TCGA-BH-A1F8-11B 1 0.002739726 11.071462 12.134105 13.323336 10.804131 9.079485 7.011227
# TCGA-BH-A1F8-01A 1 0.002739726 7.483816 7.228819 9.620220 5.459432 5.169925 8.791163
# TCGA-AC-A7VC-01A 0 0.002739726 7.098032 7.523562 9.552669 5.930737 4.169925 9.202124
# TCGA-AN-A0AM-01A 0 0.013698630 7.330917 6.491853 7.209453 4.392317 3.584963 10.099348
# 設(shè)置隨機(jī)種子,方便大家重復(fù)!
set.seed(1234)
# 設(shè)置自變量和因變量
x <- as.matrix(lasso_data[ , c(3:ncol(lasso_data))])
y <- as.matrix(Surv(lasso_data$OS.time, lasso_data$OS))
# x為自變量,矩陣中是特征數(shù)據(jù),用于構(gòu)建lasso回歸模型
# y為因變量,矩陣中中是生存分析中的生存對(duì)象
# 因?yàn)樵蹅円獦?gòu)建生存風(fēng)險(xiǎn)模型,如果只是普通的分類模型,y里面是樣本的分類信息就好啦!
# 假設(shè)我們以生存/死亡為分類依據(jù)(或者疾病/正常,都行?。?,那么y按照下面這樣寫就好啦!
# y <- as.matrix(lasso_data$OS)
# 使用glmnet()建模,參數(shù)解釋見下文
alpha1_fit <- glmnet(x, y, alpha = 1, family = 'cox', nlambda = 100)
解釋一下glmnet
函數(shù)各參數(shù)的意義:
x
: 是一個(gè)矩陣或數(shù)據(jù)框,包含自變量(特征)的觀測(cè)值。模型將使用這些自變量來預(yù)測(cè)因變量 y
。在Cox比例風(fēng)險(xiǎn)模型中,這些自變量通常是影響生存時(shí)間的因素,我們今天用的都是基因表達(dá)量。y
: 這是因變量的向量或矩陣,它包含了我們想要預(yù)測(cè)或建模的目標(biāo)變量的觀測(cè)值。我們這里用的是cox
,表示用于生存分析,因變量是生存時(shí)間(或事件發(fā)生時(shí)間)和事件狀態(tài)(如生存與死亡)。alpha
: 這是一個(gè)介于0和1之間的參數(shù),用于控制Lasso回歸和Ridge回歸之間的權(quán)衡。當(dāng) alpha = 1
時(shí),表示使用Lasso回歸(L1正則化),當(dāng) alpha = 0
時(shí),表示使用Ridge回歸(L2正則化)。如果alpha在0~1之間,則為彈性網(wǎng)絡(luò)。在這里,我們使用alpha = 1
,因?yàn)槲覀兿M褂肔asso正則化。family
參數(shù)用于指定模型所假設(shè)的概率分布族,這將影響到模型的形式和假設(shè)。不同的問題類型和數(shù)據(jù)類型需要選擇不同的分布族。下面咱們解釋幾種參數(shù)值:family = 'cox'
表示我們正在構(gòu)建一個(gè)Cox比例風(fēng)險(xiǎn)模型,用于生存分析。Cox模型是用于研究事件發(fā)生時(shí)間與各種協(xié)變量之間的關(guān)系的模型。'gaussian'
:用于連續(xù)的數(shù)值型因變量,假設(shè)因變量服從正態(tài)分布,這是線性回歸的經(jīng)典用途,也叫做普通最小二乘回歸(Ordinary Least Squares,OLS)。'binomial'
:用于二元分類問題,假設(shè)因變量是二元的,表示患病與正常或1與0的概率,這是邏輯回歸(logistic distribution)的常見用途。'poisson'
:用于計(jì)數(shù)數(shù)據(jù),假設(shè)因變量是計(jì)數(shù)數(shù)據(jù),且服從泊松分布,常用于描述事件在一段時(shí)間或空間內(nèi)的發(fā)生次數(shù)。'multinomial'
:用于多類別分類問題,假設(shè)因變量有多個(gè)類別,通常用于多類別邏輯回歸或多類別分類問題。'mgaussian'
:用于多元高斯分布,假設(shè)因變量是多維連續(xù)數(shù)值型數(shù)據(jù),且服從多元高斯分布。這可以用于多元回歸問題。
nlambda
: 這個(gè)參數(shù)表示要擬合的lambda(正則化參數(shù))的數(shù)量。Lambda是一個(gè)控制正則化強(qiáng)度的參數(shù),它控制著Lasso回歸中特征的稀疏性程度。nlambda
指定了在不同lambda值上進(jìn)行模型擬合的次數(shù),以便找到合適的正則化強(qiáng)度。一般我們使用默認(rèn)值100就可以啦!
# 繪圖
plot(alpha1_fit, xvar = 'lambda', label = TRUE)
這是回歸系數(shù)路徑圖。
圖中的每一條曲線上都有變量編號(hào),即每一條曲線代表了每一個(gè)自變量系數(shù)的變化軌跡,縱坐標(biāo)是系數(shù)的值,下橫坐標(biāo)是log λ,上橫坐標(biāo)是此時(shí)模型中非零系數(shù)的個(gè)數(shù)。
我們可以看到,隨著參數(shù)log λ增大,回歸系數(shù)(即縱坐標(biāo)值)不斷收斂,最終收斂成0,最終系數(shù)被壓縮為0的變量,說明它比較重要。例如,最上面那條代表的自變量135(綠色)在λ值很大時(shí)就有非零的系數(shù),然后隨著λ值變大不斷變小。
接下來,我們進(jìn)行交叉驗(yàn)證(cross-validation),默認(rèn)nfolds = 10
,也就是十折交叉驗(yàn)證,大家可以按需修改!
# 交叉驗(yàn)證,參數(shù)解釋見下文
set.seed(1234)
alpha1.fit.cv <- cv.glmnet(x, y, type.measure = 'deviance', alpha = 1, family = 'cox')
type.measure
參數(shù)用于指定在進(jìn)行交叉驗(yàn)證時(shí)用于度量模型性能的指標(biāo)。咱們的type.measure
被設(shè)置為 'deviance
',這是用于分類問題的一種常見指標(biāo),特別是對(duì)于邏輯回歸等模型。它表示對(duì)數(shù)似然的負(fù)二倍,可用于評(píng)估分類模型的擬合。除了 'deviance
',cv.glmnet
函數(shù)還支持其他一些常見的度量指標(biāo),具體取決于問題類型。下面解釋幾種 type.measure
的參數(shù)值:
- 'mae'(Mean Absolute Error): 這是用于回歸問題的指標(biāo),衡量了模型預(yù)測(cè)值與真實(shí)值之間的絕對(duì)差異的平均值。
- 'auc'(Area Under the ROC Curve): 這是用于二元分類問題的指標(biāo),度量了模型的分類性能,尤其是對(duì)正類別的分類準(zhǔn)確性。
- 'class':這是另一個(gè)用于分類問題的度量指標(biāo),通常與'deviance'或'auc'一起使用,以提供更多的分類性能信息。
- 'mse'(Mean Squared Error): 這是用于回歸問題的度量指標(biāo),衡量了模型預(yù)測(cè)值與真實(shí)值之間的平方差的平均值。
# 繪圖
plot(alpha1.fit.cv)
介個(gè)是Lasso回歸的交叉驗(yàn)證曲線。
X軸是懲罰系數(shù)的對(duì)數(shù) log λ,Y軸是似然偏差,Y軸越小說明方程的擬合效果越好。最上面的數(shù)字則為不同λ時(shí),方程剩下的變量數(shù)。兩條虛線代表兩個(gè)特殊的lambda(λ)值。
左邊虛線為λ min,意思是偏差最小時(shí)的λ ,代表在該lambda取值下,模型擬合效果最高。變量數(shù)是29,相比λ-se,保留下來的變量更多。
右邊虛線為λ-se,意思是最小λ右側(cè)的1個(gè)標(biāo)準(zhǔn)誤。在該λ取值下,構(gòu)建模型的擬合效果也很好,同時(shí)納入方程的個(gè)數(shù)更少,模型更簡(jiǎn)單。因此,臨床上一般會(huì)選擇右側(cè)的λ1-se作為最終方程篩選標(biāo)準(zhǔn)。我們今天選擇的數(shù)據(jù)不太好哈,變量數(shù)為100以后擬合效果就都差不多啦!不過沒事!咱主要是為了理解方法!不管結(jié)果好壞!
# 打印結(jié)果
print(alpha1.fit.cv)
# Call: cv.glmnet(x = x, y = y, type.measure = 'deviance', alpha = 1, family = 'cox')
#
# Measure: Partial Likelihood Deviance
#
# Lambda Index Measure SE Nonzero
# min 0.01289 18 11.84 0.245 29
# 1se 0.06267 1 12.03 0.238 0
# 這里我們看最后一列,29是λ min時(shí)的變量數(shù),那個(gè)0,理論上不能是0!
# 應(yīng)該是比29小的數(shù)字!這是個(gè)意外,要不咱們假裝那里是6!繼續(xù)往后講!
# 特征篩選,下面打印出來的不是我們的數(shù)據(jù)得到的結(jié)果,只是給大家展示一下
# 第一行顯示的442 x 1代表的就是442個(gè)變量,略多哈,放心!一般大家分析的時(shí)候應(yīng)該會(huì)很少!
coef(alpha1.fit.cv, s = alpha1.fit.cv$lambda.1se)
# 442 x 1 sparse Matrix of class 'dgCMatrix'
# s1
# (Intercept) -13.320908508
# CSMD3 .
# KCNIP3 0.302591830
# GNAL 0.518595837
# KCNB1 0.432713282
# SPHKAP .
# CRTAC1 .
# CDHR1 0.143984185
# TNR 0.009118150
# PTPRT .
# GABRG1 0.062742617
# ELFN2 .
# RTP5 .
# 點(diǎn)點(diǎn)就代表系數(shù)為0,系數(shù)不為0的就是最終納入的變量
# 提取特征
feature_all <- as.data.frame(as.matrix(coef(alpha1.fit.cv, s = alpha1.fit.cv$lambda.1se)))
colnames(feature_all) <- 'coff'
feature_opt <- feature_all %>% filter(abs(coff) > 0)
rownames(feature_opt)
# [1] 'NKAIN1' 'IL4I1' 'LEPR' 'FOXJ1' 'LYVE1' 'CRHBP' 'S100B'
# [8] 'FREM1' 'CST1' 'HTR1D' 'WT1' 'TNFRSF18'
# 然后我們就成功篩選出候選變量啦!接下來就可以進(jìn)行后續(xù)的分析啦!
# 比如進(jìn)一步篩選變量、構(gòu)建預(yù)測(cè)模型等等。
所以通常文章中會(huì)放入這兩張圖,表示進(jìn)行了Lasso回歸。
在這之后,我們就可以進(jìn)行后續(xù)的Cox回歸分析、列線圖、生存分析、風(fēng)險(xiǎn)因子三圖聯(lián)動(dòng)等等步驟啦!像下面這樣!
圖片來源:Li Y, Mo H, Wu S, Liu X, Tu K. A Novel Lactate Metabolism-Related Gene Signature for Predicting Clinical Outcome and Tumor Microenvironment in Hepatocellular Carcinoma. Front Cell Dev Biol. 2022 Jan 3;9:801959.這些我后面應(yīng)該都會(huì)陸續(xù)更新分享,大家可以慢慢學(xué)習(xí)喲!不過總會(huì)收到小伙伴詢問接不接生信分析的單子,主要是考慮到近期有點(diǎn)忙,沒精力接很多任務(wù),所以這里先回復(fù)一下 —— 暫時(shí)不接。我還是比較建議大家學(xué)習(xí)一下我的保姆級(jí)教程哈哈哈哈哈哈哈,感覺應(yīng)該可以包會(huì)!當(dāng)然,也有一些小伙伴們可能希望有人指導(dǎo)或者協(xié)助完成生信分析工作,如果急需的話,也可以找我,我們可以協(xié)商一下是否合適!最后,真心希望可以幫到大家! —— 一只小蠻要
結(jié)語
好啦!以上就是Lasso回歸的介紹啦!大家一定收獲滿滿吧!
其實(shí)我們?cè)诜治鲋锌赡軙?huì)先進(jìn)行單因素Cox回歸,然后依據(jù)剩下的特征數(shù)去選擇后續(xù)分析方法。如果單因素Cox回歸后剩下的特征還很多,那可以用Lasso進(jìn)行再次篩選,然后再進(jìn)行多因素Cox回歸分析,得到最終的特征用于構(gòu)建風(fēng)險(xiǎn)模型。當(dāng)然,如果剩下的特征非常少,一般就不進(jìn)行Lasso啦!大家可以依據(jù)實(shí)際情況靈活調(diào)整!
在之后的分享中,我們還會(huì)介紹單因素和多因素Cox回歸分析喲!敬請(qǐng)期待!