临床预测模型概述6-统计模型实操-Lasso回归
原创
基础知识回顾：
tps://mp.weixin.qq.com/s/pXRZ1rYUr3lwH5OlDeB0_Q
https://mp.weixin.qq.com/s/UVR6ZHCwhWqTfFBmPYPV9Q
https://mp.weixin.qq.com/s/NzLmmYReekA3E-EIonP22w
简单回顾一下 Lasso回归（Least Absolute Shrinkage and Selection Operator) ， Lasso是一种回归分析方法，通过引入L1正则化(惩罚)项对模型参数进行约束，从而实现变量选择和模型的正则化。Lasso回归通过最小化预测误差和惩罚项的和，能够将不重要的特征系数缩减为零，适用于高维数据分析，帮助防止模型过拟合。其惩罚强度由参数λ控制，λ值越大，模型越简单，选择的变量越少。通常该方法用于筛选自变量(大量的基因数据/临床参数等)，有时候也可以用于获取建模前自变量的系数。
Lasso回归可以使用glmnet包实现，研究者对该包的介绍为：Glmnet 是一个用于拟合广义线性模型和类似模型的R语言包，通过带有惩罚项的最大似然估计来实现。这种方法会在一系列（对数尺度上）的惩罚参数 lambda 值上计算 Lasso 或 Elastic Net 的正则化路径。它的算法非常快速，能有效利用输入矩阵的稀疏性。Glmnet 可以拟合线性回归、逻辑回归、多分类回归、泊松回归以及Cox回归模型，还可以处理多响应线性回归、自定义族的广义线性模型，以及Lasso回归模型。这个包还包括用于预测、绘图的函数，以及交叉验证的功能。
此外，需要知道的是除了L1正则化，还有L2正则化和弹性网络分析，如果是L1正则化就是lasso回归，L2正则化就是岭回归，弹性网络是L1和L2正则化的结合。
接下来进行Lasso回归模型筛选自变量的代码演示，其中最佳模型一般会采用10乘交叉验证法确定。
详细流程
1、加载R包和导入数据
rm(list = ls())
library(glmnet)
library(survival)
load("TCGA-LIHC_dat.Rdata")2、check数据
head(exprSet)[1:4,1:4]
# TCGA-FV-A495-01A TCGA-ED-A7PZ-01A TCGA-ED-A97K-01A TCGA-ED-A7PX-01A
# WASH7P 1.913776 1.2986076 1.967382 1.586170
# AL627309.6 3.129116 0.5606928 3.831265 1.363539
# WASH9P 2.490476 2.8140204 2.960338 2.106464
# MTND1P23 2.773335 3.4114257 2.591028 3.353850
head(meta)
# ID OS OS.time race age gender stage T N M
# TCGA-FV-A495-01A TCGA-FV-A495-01A 0 0.03333333 WHITE <=60 FEMALE II 2 <NA> 0
# TCGA-ED-A7PZ-01A TCGA-ED-A7PZ-01A 0 0.20000000 ASIAN >60 MALE II 2 <NA> 0
# TCGA-ED-A97K-01A TCGA-ED-A97K-01A 0 0.20000000 ASIAN <=60 MALE III 3 0 0
# TCGA-ED-A7PX-01A TCGA-ED-A7PX-01A 0 0.20000000 ASIAN <=60 FEMALE II 2 <NA> 0
# TCGA-BC-A3KF-01A TCGA-BC-A3KF-01A 0 0.26666667 WHITE >60 FEMALE I 1 <NA> 0
# TCGA-DD-A4NR-01A TCGA-DD-A4NR-01A 1 0.30000000 WHITE >60 FEMALE I 1 0 0
str(meta)
# 'data.frame': 365 obs. of 10 variables:
# $ ID : chr "TCGA-FV-A495-01A" "TCGA-ED-A7PZ-01A" "TCGA-ED-A97K-01A" "TCGA-ED-A7PX-01A" ...
# $ OS : Factor w/ 2 levels "0","1": 1 1 1 1 1 2 1 1 1 2 ...
# $ OS.time: num 0.0333 0.2 0.2 0.2 0.2667 ...
# $ race : Factor w/ 4 levels "ASIAN","AMERICAN INDIAN OR ALASKA NATIVE",..: 4 1 1 1 4 4 1 4 4 4 ...
# $ age : Factor w/ 2 levels "<=60",">60": 1 2 1 1 2 2 2 2 2 2 ...
# $ gender : Factor w/ 2 levels "FEMALE","MALE": 1 2 2 1 1 1 2 2 1 2 ...
# $ stage : Factor w/ 4 levels "I","II","III",..: 2 2 3 2 1 1 1 2 1 2 ...
# $ T : Factor w/ 4 levels "1","2","3","4": 2 2 3 2 1 1 1 2 1 2 ...
# $ N : Factor w/ 2 levels "0","1": NA NA 1 NA NA 1 1 NA NA 1 ...
# $ M : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...3、构建模型-二项logistic
# 假定有一个基因集
library(openxlsx)
a <- read.xlsx("~/Desktop/axlsx",colNames = T)
genes <- a$b
length(genes)
#与表达矩阵取交集后构建模型
x = t(exprSet[rownames(exprSet) %in% genes,]) #转置后的基因表达矩阵
# 请注意lasso回归时需要把事件发生情况设置成数值
meta$OS <- as.numeric(as.character(meta$OS))
y = meta$OS #提取生存状态信息
set.seed(123)
# 对于二分类数据应该选择binomial,除此之外还有“poisson”, “multinomial”, “cox”, “mgaussian”
# 当alpha设置为0则为ridge回归，将alpha设置为0和1之间则为elastic net
cvfit = cv.glmnet(x, y, alpha = 1, family="binomial") #10折交叉验证，用于调优参数，默认nfolds = 10
plot(cvfit)
fit = glmnet(x, y, alpha = 1,family = "binomial") #建模
plot(fit,label = T)
plot(fit,xvar="lambda",label = F)
co = coef(fit,s = cvfit$lambda.min) #提取最优lamda参数对应的模型的基因系数
class(co) #稀疏矩阵
coef = data.frame(gene = rownames(co),
coefficient = as.numeric(co[,1]))
coef
# gene coefficient
# 1 (Intercept) -3.47047362
# 2 CASP9 0.00000000
# 3 NLRP3 0.00000000
# 4 NLRC4 0.19304233
# 5 CHMP3 0.27393827
# 6 IL1B 0.00000000
# 7 PJVK -0.18704218
# 8 CASP8 0.00000000
# 9 CHMP2B 0.00000000
# 10 TP63 0.00000000
# 11 CASP6 0.00000000
# 12 IRF2 0.00000000
# 13 CASP3 0.00000000
# 14 IRF1 0.00000000
# 15 BAK1 0.00000000
# 16 GSDME 0.07706212
# 17 CYCS 0.06668625
# 18 NOD1 0.00000000
# 19 CHMP7 0.00000000
# 20 CHMP4C -0.09313318
# 21 GSDMD 0.00000000
# 22 NLRP6 -0.14777049
# 23 CASP4 0.00000000
# 24 CASP1 0.00000000
# 25 IL18 0.00000000
# 26 TIRAP 0.00000000
# 27 SCAF11 0.00000000
# 28 HMGB1 0.00000000
# 29 CHMP4A 0.00000000
# 30 GZMB -0.11255841
# 31 PYCARD 0.00000000
# 32 NOD2 0.37458658
# 33 NLRP1 0.00000000
# 34 TP53 -0.04195642
# 35 GSDMB 0.00000000
# 36 CHMP6 0.00000000
# 37 GPX4 0.42890846
# 38 PRKACA -0.21837358
# 39 BAX 0.00000000
# 40 NLRP2 -0.02085698
# 41 CHMP2A -0.07132867
# 42 CHMP4B 0.00000000
# 43 PLCG1 0.00000000
coef = coef[coef$coefficient!=0,]
nrow(coef)
lassoGene = coef$gene
lassoGene
# [1] "(Intercept)" "NLRC4" "CHMP3" "PJVK" "GSDME" "CYCS" "CHMP4C"
# [8] "NLRP6" "GZMB" "NOD2" "TP53" "GPX4" "PRKACA" "NLRP2"
# [15] "CHMP2A"
这张图是Lasso回归模型的系数路径图，与后面的系数路径图类似，但X轴使用的是L1范数（L1 Norm），而不是λ的对数值。
1. 轴（L1 Norm）:
● 代表模型的L1范数，即所有非零系数的绝对值之和。横坐标从左到右，表示L1范数的逐渐增大。
● 在图的顶部，标出了在对应的L1范数下，模型中有多少个变量的系数是非零的。比如在L1范数为15时，对应有42个非零系数（即42个变量被保留下来）。
2. Y轴（Coefficients）:
● 代表各个变量的回归系数值。系数值的范围可能在-1到1之间，表示每个变量在模型中的影响方向和大小。
● 每条线代表一个变量的系数。随着L1范数的增大（即正则化的减弱），一些系数逐渐从0开始增大或减小，表示这些变量被逐渐纳入模型。
3. 各条线的变化:
● 每一条线代表一个变量的系数变化情况。在线条的起点（L1范数接近0）时，大多数系数都是0，这表示强正则化使得模型变得非常简单，几乎所有变量的系数都被压缩为零。 ● 随着L1范数的增加（横坐标向右移动），正则化的力度减弱，越来越多的变量系数变得不为零，这表示更多的变量开始被纳入模型中。
4. 左侧的情况: 当L1范数较小（接近0）时，模型施加了强烈的正则化，大多数变量的系数被压缩为零。此时，模型只包含了少数几个对预测最重要的变量。
5. 右侧的情况: 随着L1范数增大，正则化减弱，更多的变量开始进入模型，系数也逐渐远离零值。L1范数越大，模型越复杂，包含的变量越多。
1. X轴（Log(λ)）：
● 横轴表示的是λ的对数值（Log(λ)）。λ是Lasso正则化中的惩罚参数，它控制了模型的稀疏性。较大的λ值意味着更强的正则化，可能会导致更多的特征系数被压缩为零。
2. Y轴（Binomial Deviance）：
● 纵轴表示的是二项分布偏差（Binomial Deviance），它是用来衡量模型在交叉验证中的表现的一个误差度量。偏差越小，模型的表现越好。
3. 红色点和灰色误差条：
● 红色点代表每个λ值下的平均二项分布偏差。
● 灰色误差条表示偏差的标准误差范围。误差条越短，说明该λ值下的模型结果越稳定。
4. 垂直虚线：
● 左侧虚线对应的是最小偏差点（min λ），即使模型误差最小的λ值。此时模型的预测性能最佳。
● 右侧虚线对应的是一个较大的λ值（λ.1se），它位于最优λ的一个标准误差范围内。选择这个λ值通常会得到一个更加稀疏的模型，同时还能保证性能不显著下降。
5.选择最优λ值：
● 最小偏差点对应的λ值通常是选择的默认值，因为它在交叉验证中表现最好。
● 如果你希望一个更加简单且稀疏的模型，可以选择右侧虚线对应的λ值（λ.1se），因为它会导致更多的系数变为零，从而简化模型。
这张图是Lasso回归模型的系数路径图（Coefficient Path Plot），展示了不同的正则化参数λ值下，每个特征变量的系数如何变化。
1. X轴（Log Lambda）：
● 横轴表示的是λ的对数值（Log Lambda）。随着λ值的变化，Lasso正则化对模型施加的惩罚力度也在变化。
● 当Log Lambda值较大（即λ值较大）时，正则化强度更大，模型会倾向于压缩更多的特征系数为零。
2. Y轴（Coefficients）：
● 纵轴表示模型中每个特征变量的系数值。不同颜色的曲线代表不同的特征变量。
● 当Log Lambda值很大时，大部分系数都被压缩为零。随着λ值减小，模型对特征的惩罚减少，更多的系数开始从零值增长。
3. 曲线解读：
● 每一条曲线代表一个特征变量的系数。曲线从左到右表示在λ值逐渐减小时，这个特征系数的变化情况。
● 曲线从接近于0（所有特征系数接近0）逐渐增长，表示随着λ值减小，特征逐渐被纳入模型。
4. 上方的数字：
● X轴顶部的数字表示对应λ值下被选中的特征数目。随着λ值减小，选择进入模型的特征数量逐渐增加。
5. 特征选择：
● 通过观察曲线的走势，你可以了解哪些特征在较大的λ值下被纳入模型，并且随着λ的减小，哪些特征系数逐渐增大或减小。
● 可以结合前面交叉验证图中的最优λ值（min λ或λ.1se）来选择特征。路径图能够直观地展示在最优λ值下哪些特征被纳入模型以及其系数大小。
print(fit）查看内部结果
print(fit)
# Call: glmnet(x = x, y = y, family = "binomial")
#
# Df %Dev Lambda
# 1 0 0.00 0.081590
# 2 1 0.38 0.074340
# 3 2 0.71 0.067740
# 4 3 1.15 0.061720
# 5 5 1.68 0.056240
# 6 5 2.26 0.051240
# 7 5 2.73 0.046690
# 8 6 3.15 0.042540
# 9 7 3.59 0.038760
# 10 9 4.05 0.035320
# 11 9 4.53 0.032180
# 12 10 4.99 0.029320
# 13 10 5.55 0.026720
# 14 12 6.16 0.024340
# 15 13 6.86 0.022180
# 16 13 7.46 0.020210
# 17 14 8.06 0.018420
# 18 14 8.58 0.016780
# 19 18 9.07 0.015290
# 20 20 9.65 0.013930
# 21 23 10.20 0.012690
# 22 25 10.82 0.011570
# 23 26 11.37 0.010540
# 24 28 11.96 0.009602
# 25 29 12.51 0.008749
# 26 30 13.02 0.007972
# 27 31 13.45 0.007264
# 28 32 13.82 0.006618
# 29 34 14.15 0.006030
# 30 34 14.44 0.005495
# 31 35 14.70 0.005007
# 32 36 14.92 0.004562
# 33 36 15.10 0.004156
# 34 36 15.26 0.003787
# 35 37 15.40 0.003451
# 36 37 15.51 0.003144
# 37 39 15.62 0.002865
# 38 39 15.71 0.002610
# 39 39 15.79 0.002378
# 40 39 15.86 0.002167
# 41 40 15.92 0.001975
# 42 41 15.97 0.001799
# 43 41 16.01 0.001639
# 44 41 16.05 0.001494
# 45 41 16.07 0.001361
# 46 41 16.10 0.001240
# 47 41 16.12 0.001130
# 48 41 16.14 0.001030
# 49 41 16.15 0.000938
# 50 42 16.17 0.000855
# 51 42 16.18 0.000779
# 52 42 16.19 0.000710
# 53 42 16.19 0.000647
# 54 42 16.20 0.000589
# 55 42 16.21 0.000537
# 56 42 16.21 0.000489
# 57 42 16.22 0.000446
# 58 42 16.22 0.000406
# 59 42 16.22 0.000370
# 60 42 16.22 0.000337
# 61 42 16.23 0.000307
# 62 42 16.23 0.000280
# 63 42 16.23 0.000255
# 64 42 16.23 0.000232在使用glmnet进行Lasso回归建模后，打印出的模型结果展示了不同λ值（Lambda）对应的模型信息，包括选择的特征数量（Df）、偏差解释率（%Dev）和λ值本身。这些信息可以帮助理解模型在不同的正则化强度下的表现。
1. Df（Degrees of Freedom）：
● 这一列显示的是在每个λ值下，模型中非零系数（特征）的数量。随着λ值减小，Lasso正则化的强度减弱，模型中纳入的特征数量增加。
● 例如，在λ = 0.081590时，模型中没有纳入任何特征（Df = 0），但在λ = 0.001361时，模型纳入了41个非零系数的特征。
2. %Dev（Percent Deviance Explained）：
● 这一列表示模型在该λ值下解释的偏差比例，即模型的拟合优度。通常，随着λ值减小（正则化减弱），模型解释的偏差比例会上升，因为更多的特征被纳入模型。
● 例如，在λ = 0.081590时，模型几乎没有解释任何偏差（%Dev = 0.00），而在λ = 0.001361时，模型解释了16.07%的偏差。
3. Lambda：
● 这一列显示的是不同的λ值。λ值越大，Lasso正则化的强度越大，导致更多的特征系数被压缩为零；λ值越小，正则化强度减弱，更多的特征被纳入模型中。
● 例如，在λ = 0.081590时，模型的正则化强度较大，没有特征被纳入模型，而在λ = 0.001361时，正则化强度较小，模型纳入了41个特征。
4. 模型选择：
● 可以根据%Dev（解释的偏差比例）选择一个合适的λ值。通常，会希望选择一个能够解释足够多偏差（%Dev较高），同时避免过多特征（Df较低）的λ值。
5. 交叉验证：
● 通常会使用交叉验证来选择一个最优的λ值。交叉验证会给研究者提供两个有用的λ值：lambda.min（使交叉验证误差最小的λ值）和lambda.1se（在最优误差内的最大λ值，通常会得到更稀疏的模型）。
6. 绘制路径图：
● 可以绘制系数路径图或交叉验证曲线来直观地查看模型在不同λ值下的表现，从而更好地选择合适的λ值。
4. 构建模型-cox
# 假定有一个基因集
library(openxlsx)
a <- read.xlsx("~/Desktop/a.xlsx",colNames = T)
genes <- a$b
length(genes)
#与表达矩阵取交集后构建模型
x = t(exprSet[rownames(exprSet) %in% genes,]) #转置后的基因表达矩阵
# 请注意lasso回归时需要把事件发生情况设置成数值
meta$OS <- as.numeric(as.character(meta$OS))
y = Surv(meta$OS.time,meta$OS) #生存信息
set.seed(10210)
cvfit = cv.glmnet(x, y, alpha = 1,family="cox") #10折交叉验证，用于调优参数
plot(cvfit)
fit = glmnet(x, y, alpha = 1,family = "cox") #建模
#plot(fit,label = T)
plot(fit,xvar="lambda",label = F)
co = coef(fit,s = cvfit$lambda.min) #提取最优lamda参数对应的模型的基因系数
class(co) #稀疏矩阵
coef = data.frame(gene = rownames(co),
coefficient = as.numeric(co[,1]))
head(coef)
# gene coefficient
# 1 CASP9 0.0000000
# 2 NLRP3 0.0000000
# 3 NLRC4 0.1719762
# 4 CHMP3 0.0000000
# 5 IL1B 0.0000000
# 6 PJVK -0.1242998
coef = coef[coef$coefficient!=0,]
nrow(coef)
lassoGene = coef$gene
lassoGene
# [1] "NLRC4" "PJVK" "CASP8" "BAK1" "GSDME" "NLRP6" "SCAF11" "GZMB" "NOD2" "TP53" "GPX4"
# [12] "CHMP2A"cox-Lasso跟logistics-Lasso十分相似修改其中的几条代码即可。
参考资料：
1、glmnet：https://glmnet.stanford.edu/articles/glmnet.html
2、阿越老师：https://ayueme.github.io/R_clinical_model/feature-selection_lasso.html
3、生信小白要知道：https://mp.weixin.qq.com/s/kSrr6regfAtX4Bw6gSvmgw
致谢：感谢曾老师以及生信技能树团队全体成员。
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟
- END -
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
