ixxmu
commented
1 year ago 看完还不会来揍/找我 | 基于 LASSO 回归筛选变量以构建预测模型 | 附完整代码 + 注释 by 生信小白要知道
基于 LASSO 回归筛选变量以构建预测模型
在构建疾病风险评估和预测模型时，我们会想要纳入尽可能多的指标去构建模型。这么做虽然会更准确，但却不实用。因为医生不可能让所有刚入院的患者把所有的检查都做一遍。我们更希望在尽可能知道较少指标的情况下而不损失风险评估准确性，这样才可以实现效益和效率最大化。又或者，在大量基因中选取可能对疾病有重要影响的基因（比如我们前面在看完还不会来揍我 | 差异分析三巨头 —— DESeq2、edgeR 和 limma 包 | 附完整代码 + 注释中得到的差异基因），用于构建风险预测模型。这个时候，我们就可以使用 LASSO 回归在众多变量中进行筛选，获取重要变量，也就是特征，用于构建模型啦！今天我们就来分享一下具体要怎么做喽！
我们先介绍一下Lasso回归的原理，大家如果不需要的话，也可以跳过这部分直接去看用法和代码部分哟！
Lasso回归（Least Absolute Shrinkage and Selection Operator）是一种用于线性回归问题的正则化方法，它可以用来降低模型的复杂性，防止过拟合，并选择重要的特征变量。
LASSO 回归的原理
Lasso回归的原理基于在线性回归（想深入了解线性回归的小伙伴们，详见跟我一起啃西瓜书系列中的第三章线性模型3.2小节）的基础上添加一个L1正则化项（不要慌！后面会讲！），以约束模型的参数估计。L1正则化项是参数向量中绝对值之和，它的存在迫使许多参数变为零，从而实现了变量选择（特征筛选）的功能。Lasso回归的优化目标是最小化以下损失函数：
其中：
是样本数量。 是特征变量的数量。 是第 个观测值的目标值。 是第 个观测值的第 个特征变量的值。 是第 个特征变量的系数，当其为 0 时， 就被踢出方程咯，这样就达到筛选变量的效果啦！ 是截距项。 是正则化参数，用于控制正则化的强度。较大的 值将导致更多的参数估计为零。
简单介绍一下正则化。
正则化（Regularization）是一种在机器学习和统计建模中常用的技术，旨在帮助模型更好地泛化到未见过的数据，避免过拟合问题。过拟合（overfitting）是指模型在训练数据上表现得非常好，但在未见过的测试数据上表现较差的情况。简单来说，模型过度“记忆”了训练数据的细节和噪声，导致在新数据上的预测表现不佳。过拟合通常出现在模型过于复杂或者训练数据较少的情况下。就好像是一个学生过度死记硬背了大量的题目答案，但实际上并没有理解问题的本质。因此，当遇到新的类似问题时，虽然能够准确回答原有题目，但对于新问题的推广能力较差。
正则化就好像是在机器学习中的模型训练过程中施加的一种规则，它会对模型进行“限制”或“惩罚”，以确保它不会变得过于复杂。
想象一下，你正在尝试通过拟合一条曲线来拟合一组散点数据，但你担心曲线会过于弯曲，导致在散点之间出现怪异的起伏。正则化就像是在你的曲线拟合问题中引入一个规则，告诉你不要让曲线变得太弯曲。岭回归（L2正则化）和Lasso回归（L1正则化）就是为了解决线性回归中的过拟合问题而出现的，弹性网络是L1和L2正则化的结合，可以平衡两者的影响。
L1正则化（也称为Lasso正则化）：它告诉你，如果曲线中有一些不太重要的部分，就将它们收缩到零，从而削减模型的复杂度。这相当于在模型中删除不必要的特征或参数，使模型更简单。 L2正则化（也称为岭正则化）：它告诉你，如果曲线的某些部分变得过于陡峭，就要降低它们的陡峭程度。这使得曲线的参数不能太大，从而控制了模型的复杂度（在本文的后半部分，L2正则化马上就要登场啦）。 正则化的目标是在模型的损失函数中添加一个额外的项，这个项与模型参数相关，通过这个项来限制参数的大小或稀疏性，从而防止模型在训练数据上过度拟合。这种“限制”或“惩罚”有助于模型更好地泛化到新数据，降低了过拟合的风险。
LASSO 回归与岭回归、弹性网络的区别
Lasso回归（Least Absolute Shrinkage and Selection Operator）：Lasso回归引入了L1正则化项，通过约束参数向量的绝对值之和来实现稀疏性。这使得Lasso能够执行特征选择，将许多参数估计为零，从而简化模型。Lasso回归倾向于在高维数据中选择少数关键特征，提高了模型的可解释性。
Lasso回归的优化目标如下：
岭回归（Ridge Regression）：与Lasso不同，岭回归使用L2正则化项，其目标是最小化参数的平方和，而不是绝对值和。这意味着岭回归不会将参数估计为零，但会减小参数的绝对值，从而减小参数的大小。因此岭回归适合在存在多个相关特征的情况下稳定模型，但不执行特征选择。
岭回归的优化目标如下：
弹性网络（Elastic Net）：弹性网络是L1和L2正则化的组合，综合了岭回归和Lasso的优点。它的损失函数包括L1和L2正则化项，可以同时控制参数的大小和稀疏性，在特征选择和参数缩减之间进行权衡。弹性网络适用于数据集中既存在高度相关特征又需要进行特征选择的情况。
弹性网络的优化目标如下：
为何选择 LASSO 回归进行特征筛选
选择在什么时候使用岭回归、Lasso回归或弹性网络来筛选变量，其实取决于具体情况和数据性质，有一些指导原则大家可以了解一下。
当你怀疑有多个特征对目标变量有影响，而且这些特征可能存在高度相关性时，可以考虑使用岭回归：岭回归通过L2正则化可以减小参数的绝对值，但不会将参数估计为零，因此适合在存在多个相关特征的情况下稳定模型。它不会执行严格的特征选择，而是缩小了参数的幅度。 当你怀疑只有少数几个特征对目标变量有显著影响，并且希望执行特征选择时，可以考虑使用Lasso回归：Lasso回归通过L1正则化倾向于将许多参数估计为零，从而实现了特征选择的功能。如果你想要简化模型并提高可解释性，Lasso通常是一个不错的选择。 当你既想要特征选择，又想要保留某些相关特征时，可以考虑使用弹性网络：弹性网络结合了L1和L2正则化，允许你在特征选择和参数缩减之间进行权衡。这在存在多个相关特征，但你仍然希望进行一定程度的特征选择时非常有用。通过调整两个正则化参数的值，你可以控制模型的行为。
总的来说，选择哪种正则化方法应该基于以下考虑：
数据性质：了解数据的特征和相关性，以确定哪种正则化方法更适合。 目标：确定你是否更关注特征选择、参数缩减或平衡两者。 交叉验证：通常需要使用交叉验证来选择适当的正则化参数，以确保模型的性能。
通过上面的介绍，大家应该对我们选择Lasso进行特征筛选有了一定了解。一般我们在临床预测模型构建过程中，都会期望找到相对较少的特征数量，希望用更少的特征就可以对目标变量产生较大影响，且希望模型更具可解释性。在临床预测中，可能存在大量的生物医学特征，其中一些可能不是很重要或相关。Lasso可以自动将某些特征的系数估计为零，从而实现特征选择，帮助减少模型中的不必要的特征，这样生成的模型更容易解释，只包含最重要的特征，更便于医生或临床决策者理解模型的预测结果。而且Lasso在样本数量相对较少的情况下也可以工作良好，尤其是在高维数据中，它有助于防止过拟合，即使数据点有限（试问大家能拿到多少有效样本哈哈哈哈哈哈哈呜呜呜呜呜呜呜呜呜呜，帮你们哭一下）。所以，我们见到的临床模型多数是基于Lasso回归构建的！
当然！咱也不能无脑Lasso，还是要多尝试滴！万一其他模型更合适呢！
这个时候，可能会有小伙伴有疑问，那既然弹性网络是L1和L2正则化的组合，看起来似乎它考虑更全面呐，我们为什么不直接用它来进行特征筛选呢？
其实，当数据集中存在高度相关的特征时，Lasso有可能随机选择其中的一个特征并将其他特征的系数设为零，这可能导致结果不稳定，因为具体选择哪个特征取决于模型的随机性。这在高维数据中尤其成问题。而弹性网络可以平衡这两种正则化的影响，既能够进行特征选择，又能够处理多重共线性问题。它不会像Lasso那样过于随机地选择特征，所以更加稳定。
但是吧，弹性网络也不一定在所有数据中都表现最好。当数据量较大、质量较好，且问题适用于线性模型时，弹性网络可能表现良好。当数据量较小、数据质量较差，或问题具有高度非线性性质时，弹性网络的效果可能就会表现不如其他模型。
在实际应用中，我们还是要尝试不同的模型（包括弹性网络、岭回归、Lasso等）并使用交叉验证来确定哪种模型在给定问题上表现最佳。
LASSO 回归在预测模型中的用法
用法一
这也是最常见的用法，采用Lasso进行回归，根据十折交叉验证，筛选得到候选预测因子。然后用候选预测因子继续做后续的多因素Logsitic回归或者多因素Cox回归，得到最终的临床预测模型，用于后续的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回归，对19个预测因子进行筛选，最终选得7个预测因子。然后对7个预测因子构建多因素Logistic回归模型。
然后作者就基于这7个因素构建列线图（Nomogram），并进行了后续3个度的评价。
这些图的绘制，我们后续都会陆续更新哟！
用法二
当我们研究的因素较多，可以先进行单因素Logistic或Cox回归，先筛选一批可能的预测因子;然后再采用Lasso回归进行筛选，将筛选得到的预测因子，再次进行多因素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.
这篇文章通过单变量分析发现59个因素，基于专业背景知识又加了5个 P > 0.05 的因素（有的时候我们建模，发现某个专业背景上有意义的变量，却没有进入多因素分析阶段，就可以通过这样的表述把它加进去），共64个因素。然后对这个64因素进行Lasso筛选，发现了6个系数不为零，也就是有意义的预测因子。
作者对这6个预测因子，做了多因素Cox回归，最终发现4个因子，并使用这4个因子构建最终的最优模型。
然后构建了4个因子的Nomo图以及后续的3个度的验证等等！
用法三
这种用法相对少见，就是直接用Lasso-Logit或者Lasso-Cox进行变量筛选，筛选后，用最小误差解构建模型，无需再对Lasso筛选得到的多个因子，再进行多因素的Logistic或Cox进行分析。
LASSO 回归存在的问题
尽管Lasso回归应用广泛且优点多多，但不可避免的，它也存在一些问题和局限性：
参数选择问题：选择合适的正则化参数 是一个挑战，需要进行交叉验证或其他模型选择方法。 估计偏差：Lasso倾向于对参数进行严格的稀疏化，可能导致估计偏差，尤其是在特征之间存在高度相关性时。 难以处理大规模数据集：对于大规模数据集，Lasso的计算成本可能很高，因为要考虑所有特征的组合。
LASSO 回归代码实现及结果解读
数据用到的是我们之前看完还不会来揍我 | 差异分析三巨头 —— DESeq2、edgeR 和 limma 包 | 附完整代码 + 注释中得到的差异基因表达矩阵，进行Lasso回归需要样本分类信息，我们使用生存数据，大家可以通过 UCSC Xena（http://xena.ucsc.edu/）进行下载。从 UCSC Xena 下载数据的方法详见：看完还不会来揍/找我 | TCGA 与 GTEx 数据库联合分析 | 附完整代码 + 注释的数据准备部分。
如果大家想快快实战，不想耗费时间去运行前面的差异分析步骤得到差异分析表达矩阵或者下载生存数据，完全没问题！我已经把用到的数据上传到了GitHub，大家可以在公众号后台回复
lasso回归，即可得到存放这些数据的链接。
咱们正式开始！
代码看起来有点长，但其实正式的Lasso回归部分一点都不长！前面大多是数据处理部分，大家如果有现成的数据，直接从”现在正式开始lasso回归！”这里开始看就好啦！
################################# Lasso 回归 ###################################
# 加载包
# 如果大家在运行的时候发现某个函数不存在或者不能用，可能是忘记加载包啦
# 我这里如果没有写全的话，先给大家说声抱歉，因为有时候我的环境里可能已经加载了就直接用了
# 所以可能会忘记写某个包，大家不要慌，如果遇到了且不知道函数来自哪个包的时候
# 咱们去百度谷歌噼里啪啦查一下就好！不怕！你是最棒的！
library(glmnet)
library(survival)
library(dplyr)
### 加载数据
# 咱就今天使用的数据是BRCA数据集中的癌症样本，如果大家想了解这个数据是如何得出的，
# 可以查看：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
# 定义函数挑选差异基因
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
# 三个包的结果取交集得到442个差异基因
# 获取取交集得到的差异基因的表达矩阵
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
# 加载生存数据
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
# 我们再看一下表达数据
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
# 咱们需要转置一下，然后合并！
# 整合表达数据和生存数据
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
# 咱们整理一下下，样本名设为行名，剔除没用的列
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
# 现在的数据有1194行（样本）、444列（基因），第一列为OS（样本结局，1为死亡，0为生存）
# 第二列为OS.time（生存时间，单位为day），后面的442列为基因的表达量。
# 现在正式开始lasso回归！
# 把生存时间单位从day转换为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
# 设置随机种子，方便大家重复！
set.seed(1234)
# 设置自变量和因变量
x <- as.matrix(lasso_data[ , c(3:ncol(lasso_data))])
y <- as.matrix(Surv(lasso_data$OS.time, lasso_data$OS))
# x为自变量，矩阵中是特征数据，用于构建lasso回归模型
# y为因变量，矩阵中中是生存分析中的生存对象
# 因为咱们要构建生存风险模型，如果只是普通的分类模型，y里面是样本的分类信息就好啦！
# 假设我们以生存/死亡为分类依据（或者疾病/正常，都行！），那么y按照下面这样写就好啦！
# y <- as.matrix(lasso_data$OS)
# 使用glmnet()建模，参数解释见下文
alpha1_fit <- glmnet(x, y, alpha = 1, family = "cox", nlambda = 100)
解释一下glmnet函数各参数的意义：
x: 是一个矩阵或数据框，包含自变量（特征）的观测值。模型将使用这些自变量来预测因变量y。在Cox比例风险模型中，这些自变量通常是影响生存时间的因素，我们今天用的都是基因表达量。y: 这是因变量的向量或矩阵，它包含了我们想要预测或建模的目标变量的观测值。我们这里用的是cox，表示用于生存分析，因变量是生存时间（或事件发生时间）和事件状态（如生存与死亡）。alpha: 这是一个介于0和1之间的参数，用于控制Lasso回归和Ridge回归之间的权衡。当alpha = 1时，表示使用Lasso回归（L1正则化），当alpha = 0时，表示使用Ridge回归（L2正则化）。如果alpha在0~1之间，则为弹性网络。在这里，我们使用alpha = 1，因为我们希望使用Lasso正则化。family参数用于指定模型所假设的概率分布族，这将影响到模型的形式和假设。不同的问题类型和数据类型需要选择不同的分布族。下面咱们解释几种参数值：family = "cox"表示我们正在构建一个Cox比例风险模型，用于生存分析。Cox模型是用于研究事件发生时间与各种协变量之间的关系的模型。"gaussian"：用于连续的数值型因变量，假设因变量服从正态分布，这是线性回归的经典用途，也叫做普通最小二乘回归（Ordinary Least Squares，OLS）。"binomial"：用于二元分类问题，假设因变量是二元的，表示患病与正常或1与0的概率，这是逻辑回归（logistic distribution）的常见用途。"poisson"：用于计数数据，假设因变量是计数数据，且服从泊松分布，常用于描述事件在一段时间或空间内的发生次数。"multinomial"：用于多类别分类问题，假设因变量有多个类别，通常用于多类别逻辑回归或多类别分类问题。"mgaussian"：用于多元高斯分布，假设因变量是多维连续数值型数据，且服从多元高斯分布。这可以用于多元回归问题。nlambda: 这个参数表示要拟合的lambda（正则化参数）的数量。Lambda是一个控制正则化强度的参数，它控制着Lasso回归中特征的稀疏性程度。nlambda指定了在不同lambda值上进行模型拟合的次数，以便找到合适的正则化强度。一般我们使用默认值100就可以啦！
# 绘图
plot(alpha1_fit, xvar = "lambda", label = TRUE)
这是回归系数路径图。
图中的每一条曲线上都有变量编号，即每一条曲线代表了每一个自变量系数的变化轨迹，纵坐标是系数的值，下横坐标是log λ，上横坐标是此时模型中非零系数的个数。
我们可以看到，随着参数log λ增大，回归系数（即纵坐标值）不断收敛，最终收敛成0，最终系数被压缩为0的变量，说明它比较重要。例如，最上面那条代表的自变量135（绿色）在λ值很大时就有非零的系数，然后随着λ值变大不断变小。
接下来，我们进行交叉验证（cross-validation），默认nfolds = 10，也就是十折交叉验证，大家可以按需修改！
# 交叉验证，参数解释见下文
set.seed(1234)
alpha1.fit.cv <- cv.glmnet(x, y, type.measure = "deviance", alpha = 1, family = "cox")
type.measure 参数用于指定在进行交叉验证时用于度量模型性能的指标。咱们的type.measure 被设置为 "deviance"，这是用于分类问题的一种常见指标，特别是对于逻辑回归等模型。它表示对数似然的负二倍，可用于评估分类模型的拟合。除了 "deviance"，cv.glmnet 函数还支持其他一些常见的度量指标，具体取决于问题类型。下面解释几种 type.measure的参数值：
"mae"（Mean Absolute Error）: 这是用于回归问题的指标，衡量了模型预测值与真实值之间的绝对差异的平均值。 "auc"（Area Under the ROC Curve）: 这是用于二元分类问题的指标，度量了模型的分类性能，尤其是对正类别的分类准确性。 "class"：这是另一个用于分类问题的度量指标，通常与"deviance"或"auc"一起使用，以提供更多的分类性能信息。 "mse"（Mean Squared Error）: 这是用于回归问题的度量指标，衡量了模型预测值与真实值之间的平方差的平均值。
# 绘图
plot(alpha1.fit.cv)
介个是Lasso回归的交叉验证曲线。
X轴是惩罚系数的对数 log λ，Y轴是似然偏差，Y轴越小说明方程的拟合效果越好。最上面的数字则为不同λ时，方程剩下的变量数。两条虚线代表两个特殊的lambda（λ）值。
左边虚线为λ min，意思是偏差最小时的λ ，代表在该lambda取值下，模型拟合效果最高。变量数是29，相比λ-se，保留下来的变量更多。
右边虚线为λ-se，意思是最小λ右侧的1个标准误。在该λ取值下，构建模型的拟合效果也很好，同时纳入方程的个数更少，模型更简单。因此，临床上一般会选择右侧的λ1-se作为最终方程筛选标准。我们今天选择的数据不太好哈，变量数为100以后拟合效果就都差不多啦！不过没事！咱主要是为了理解方法！不管结果好坏！
# 打印结果
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时的变量数，那个0，理论上不能是0！
# 应该是比29小的数字！这是个意外，要不咱们假装那里是6！继续往后讲！
# 特征筛选，下面打印出来的不是我们的数据得到的结果，只是给大家展示一下
# 第一行显示的442 x 1代表的就是442个变量，略多哈，放心！一般大家分析的时候应该会很少！
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 .
# 点点就代表系数为0，系数不为0的就是最终纳入的变量
# 提取特征
feature_all <- as.data.frame(as.matrix(coef(alpha1.fit, s = alpha1.fit$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"
# 然后我们就成功筛选出候选变量啦！接下来就可以进行后续的分析啦！
# 比如进一步筛选变量、构建预测模型等等。
所以通常文章中会放入这两张图，表示进行了Lasso回归。
在这之后，我们就可以进行后续的Cox回归分析、列线图、生存分析、风险因子三图联动等等步骤啦！像下面这样！
这些我后面应该都会陆续更新分享，大家可以慢慢学习哟！不过总会收到小伙伴询问接不接生信分析的单子，主要是考虑到近期有点忙，没精力接很多任务，所以这里先回复一下 —— 暂时不接。我还是比较建议大家学习一下我的保姆级教程哈哈哈哈哈哈哈，感觉应该可以包会！当然，也有一些小伙伴们可能希望有人指导或者协助完成生信分析工作，如果急需的话，也可以找我，我们可以协商一下是否合适！最后，真心希望可以帮到大家！ —— 一只小蛮要
结语
好啦！以上就是Lasso回归的介绍啦！大家一定收获满满吧！
其实我们在分析中可能会先进行单因素Cox回归，然后依据剩下的特征数去选择后续分析方法。如果单因素Cox回归后剩下的特征还很多，那可以用Lasso进行再次筛选，然后再进行多因素Cox回归分析，得到最终的特征用于构建风险模型。当然，如果剩下的特征非常少，一般就不进行Lasso啦！大家可以依据实际情况灵活调整！
在之后的分享中，我们还会介绍单因素和多因素Cox回归分析哟！敬请期待！
如果大家觉得我分享的内容对你有一丢丢帮助的话，可以关注点赞收藏分享阿巴阿巴来一波哟！
参考资料
https://zhuanlan.zhihu.com/p/426162272 https://zhuanlan.zhihu.com/p/509021086 https://www.jianshu.com/p/617db057df37 https://www.cnblogs.com/zingp/p/10375691.html https://www.bilibili.com/video/BV1kX4y1g7jp https://mp.weixin.qq.com/s/aabknP4xq2j5JyKWCSEiUQ https://www.hindawi.com/journals/jdr/2020/7261047/ https://onlinelibrary.wiley.com/doi/10.1002/cam4.2896 https://mp.weixin.qq.com/s/EpVcu7x-EI2bzGjH5lG19Q https://mp.weixin.qq.com/s/sCDEP0QmdYRp70LUMPcXIw https://mp.weixin.qq.com/s/5E1WXTawQPd4JsaShAWrxA
https://mp.weixin.qq.com/s/kSrr6regfAtX4Bw6gSvmgw