回歸我們并不陌生,線性回歸和最小二乘法,邏輯回歸和最大似然法,這些都是我們耳熟能詳?shù)氖挛?,在生物信息學(xué)中的應(yīng)用也比較廣泛, 回歸中經(jīng)常出現(xiàn)兩類問題,欠擬合和過擬合。
對于欠擬合,簡單而言就是我們考慮的少了,一般通過在回歸模型中增加自變量或者擴大樣本數(shù)量來解決;對于過擬合,簡單而言就是考慮的太多了,模型過于復(fù)雜了,這時候可以對已有的自變量進行篩選,在代價函數(shù)中增加懲罰項來限制模型的復(fù)雜度,增加的懲罰項我們稱之為正則化,正則化常用的有L1正則化和L2正則化,
所謂正則化Regularization, 指的是在回歸模型代價函數(shù)后面添加一個約束項, 在線性回歸模型中,有兩種不同的正則化項
1. 所有參數(shù)絕對值之和,即L1范數(shù),對應(yīng)的回歸方法叫做Lasso回歸
2. 所有參數(shù)的平方和,即L2范數(shù),對應(yīng)的回歸方法叫做Ridge回歸,嶺回歸
lasso回歸對應(yīng)的代價函數(shù)如下
嶺回歸對應(yīng)的代價函數(shù)如下
紅框標記的就是正則項,需要注意的是,正則項中的回歸系數(shù)為每個自變量對應(yīng)的回歸系數(shù),不包含回歸常數(shù)項。
在預(yù)后建模的文章中,我們需要針對多個marker基因的表達量匯總形成一個指標,使用該指標來作為最終的maker, 而這個指標在文章中被稱之為各種risk score, 比如NAD+基因的預(yù)后模型,構(gòu)建的maker就叫做NPRS, 全稱的解釋如下
The NAD+ metabolism-related prognostic risk score (NPRS) of each sample was calculated using the formula: NPRS = ΣExp (mRNA?) × Coefficient (mRNA?)
所以各種的預(yù)后建模,其實都是lasso回歸技術(shù)在生物信息學(xué)領(lǐng)域的應(yīng)用。注意觀察上述的Lasso回歸代價函數(shù),,可以看到有一個未知數(shù)λ, 這個參數(shù)是一個懲罰項的系數(shù),數(shù)值越大,懲罰項對應(yīng)的影響就越大,我們求解的目標是代價函數(shù)值最小,λ = 0時,懲罰項失去意義,代價函數(shù)變成了普通的線性回歸,而λ過大,懲罰項的影響被放的過大,過小時,懲罰項又失去了原本的意義,所以使用lasso回歸,第一個問題是設(shè)置合理的λ 值。
這個λ 值 如何設(shè)置呢?最簡單的辦法是找到兩個隊列,訓(xùn)練集和驗證集,適應(yīng)一系列的λ值對訓(xùn)練集進行建模,觀察模型在驗證集上的表現(xiàn),然后選擇在驗證集上表現(xiàn)最佳模型的λ值,當沒有額外的驗證集時,就只能通過交叉驗證的方式將數(shù)據(jù)集人工劃分為訓(xùn)練集和驗證集,然后進行分析。在NAD+的文獻中,也是采用了10折交叉驗證的方式
In the training cohort, using the Least Absolute Shrinkage And Selection Operator (LASSO) regression with 10-fold cross-validated to screen out NMRGs associated with survival in ALS patients.
具體到實際操作,使用的是glmnet這個R包
Here, the glmnet package was applied to determine the optimal lambda value corresponding to the minimum of the error mean via cross-validation.
官方鏈接如下
https://glmnet.stanford.edu/
正則項本身只是一個代價函數(shù)中的添加項,所以其應(yīng)用范圍不僅局限于線性回歸,邏輯回歸,cox回歸都支持,所以glmnet這個R包也支持多種回歸模型的正則化處理。對于cox回歸而言,其用法可以參考如下鏈接
https://glmnet.stanford.edu/articles/Coxnet.html
基本的操作步驟如下
1. 準備輸入文件
包括自變量和因變量,自變量是一個矩陣,每一行表示一個患者,每一列表示一個自變量;因變量也是一個矩陣,共兩列,分別為代表生存信息的time加status, 代碼如下
> library(glmnet)
載入需要的程輯包:Matrix
Loaded glmnet 4.1-2
> library(survival)
> data(CoxExample)
> x <- CoxExample$x
> y <- CoxExample$y
# 自變量數(shù)據(jù),每一行表示一個患者,每一列表示一個自變量
> head(x[, 1:5])
[,1] [,2] [,3] [,4] [,5]
[1,] -0.8767670 -0.6135224 -0.56757380 0.6621599 1.82218019
[2,] -0.7463894 -1.7519457 0.28545898 1.1392105 0.80178007
[3,] 1.3759148 -0.2641132 0.88727408 0.3841870 0.05751801
[4,] 0.2375820 0.7859162 -0.89670281 -0.8339338 -0.58237643
[5,] 0.1086275 0.4665686 -0.57637261 1.7041314 0.32750715
[6,] 1.2027213 -0.4187073 -0.05735193 0.5948491 0.44328682
# 因變量數(shù)據(jù),生存數(shù)據(jù)的因變量為time加status
> head(y)
time status
[1,] 1.76877757 1
[2,] 0.54528404 1
[3,] 0.04485918 0
[4,] 0.85032298 0
[5,] 0.61488426 1
[6,] 0.29860939 0
2. 交叉驗證
通過交叉驗證,選擇最佳的λ值。在選擇λ值時,我們需要指定評價指標,就是根據(jù)評價指標的值來選擇最佳模型和最佳λ值,對應(yīng)的是typpe.measure參數(shù),對于cox模型而言,只支持以下兩種指標
1. deviance
2. C-index
評價指標c-index的代碼如下
> cvfit <- cv.glmnet(x, y, family = "cox", type.measure = "C", nfolds = 10)
> plot(cvfit)
輸出如下
評價指標deviance的代碼如下
> cv.glmnet(x, y, family = "cox", type.measure = "deviance", nfolds = 10)
> plot(cvfit)
輸出如下
在上述圖片中,橫坐標為log λ值,縱坐標為每個λ值對應(yīng)的評價指標,用error bar的形式展現(xiàn)了多個模型評價指標的均值+標準誤,可以看到在圖中有兩條垂直的虛線,左邊的虛線對應(yīng)評價指標最佳的λ值,即lambda.min, c-index值越大越好,deviance值越小越好;右邊的虛線表示評價指標在最佳值1個標準誤范圍的模型的λ值,即lambda.1se, 通過以下方式可以提取對應(yīng)的值
> cvfit$lambda.min
[1] 0.01749823
> cvfit$lambda.1se
[1] 0.04868986
通過print函數(shù)可以看到交叉驗證的關(guān)鍵信息
> print(cvfit)
Call: cv.glmnet(x = x, y = y, type.measure = "deviance", nfolds = 10, family = "cox")
Measure: Partial Likelihood Deviance
Lambda Index Measure SE Nonzero
min 0.01750 29 13.08 0.06221 15
1se 0.04869 18 13.14 0.05369 10
通過coef函數(shù)可以顯示自變量的回歸系數(shù),可以看到很多自變量的回歸系數(shù)都是0,就意味著這些自變量被過濾掉了
> coef(cvfit, s = cvfit$lambda.1se)
30 x 1 sparse Matrix of class "dgCMatrix"
1
V1 0.38108115
V2 -0.09838545
V3 -0.13898708
V4 0.10107014
V5 -0.11703684
V6 -0.39278773
V7 0.24631270
V8 0.03861551
V9 0.35114295
V10 0.04167588
V11 .
V12 .
V13 .
V14 .
V15 .
V16 .
V17 .
V18 .
V19 .
V20 .
V21 .
V22 .
V23 .
V24 .
V25 .
V26 .
V27 .
V28 .
通過交叉驗證,在選擇最佳λ值的同事,也確定了最佳的回歸模型,通過coef提取回歸系數(shù),我們就得到了最終的回歸模型。
聯(lián)系客服