从原理到实战，一文带你掌握LASSO回归在临床数据中的应用｜实战分享·24-06-14

最新推荐文章于 2024-08-14 17:31:38 发布

罗小罗同学

最新推荐文章于 2024-08-14 17:31:38 发布

阅读量2.3k

点赞数 28

分类专栏：项目复现文章标签：回归 kotlin 数据挖掘人工智能医学人工智能精准肿瘤学 LASSO

本文链接：https://blog.csdn.net/qq_45404805/article/details/139677550

版权

项目复现专栏收录该内容

12 篇文章 9 订阅

订阅专栏

小罗碎碎念

今天这期文章带大家系统的解决一个问题——LASSO回归，从原理到实战，再到实际文献的应用，我可以说这应该是目前对LASSO最系统、最详细、排版最好的一篇推文，没有之一。

小罗是工科出身，其实写这一方面的文章才真正是我自己所属的领域（but太耗时了，大家不到需要的时候基本没耐心看）。

另外，每一行代码都做了解析，并且配了非常多的图进行注解，只要你跟着一步一步走，一定能彻底掌握这一部分的内容！！

一、什么是LASSO？

Lasso回归，全称最小绝对收缩和选择算子（Least Absolute Shrinkage and Selection Operator），是一种用于估计稀疏系数的线性回归模型。

它由Robert Tibshirani在1996年提出，作为对普通最小二乘法（Ordinary Least Squares, OLS）的一种改进。Lasso通过在损失函数中添加一个绝对值惩罚项，来促使一些系数变为零，从而实现变量选择和正则化的效果。

1-1：发展历史

Lasso回归由Robert Tibshirani在1996年提出，是对线性回归模型的一种改进。Lasso是“Least Absolute Shrinkage and Selection Operator”的缩写，意为“最小绝对收缩和选择算子”。它通过引入L1正则化项来实现变量的选择和系数的压缩。

在Lasso之前，1982年Cox和Olshen提出了一种类似的方法，称为Non-negative Garrote（非负加勒特）。这种方法在正则化项中使用了一个非负性约束，而Lasso则放宽了这一约束，允许系数为零，从而实现变量的选择。

Lasso算法的发展受到了斯坦福大学统计学教授Bradley Efron提出的bootstrap方法的影响。Efron是Tibshirani的导师，他的bootstrap方法为Lasso的提出提供了理论基础。

Lasso回归的发展是统计学习领域中的一个重要里程碑，它的出现是为了解决普通最小二乘法在数据特征多且存在多重共线性时表现不佳的问题。

在Lasso之前，岭回归（Ridge Regression）是解决多重共线性问题的常用方法，但它不会产生稀疏解，即不会自动进行特征选择。Lasso则通过引入L1惩罚项（即系数的绝对值之和）来实现稀疏性，这使得Lasso在特征选择方面非常有用。

知识点补充：`L1正则化项`

L1正则化项是一种在机器学习和统计建模中常用的正则化技术，特别是在处理具有大量特征的数据集时。它通过向损失函数添加一个基于模型参数绝对值的惩罚项来实现。

L1正则化项的主要目的是促进模型的稀疏性，即在模型的参数中产生更多的零值，从而实现特征选择。

L1正则化项的定义

L1正则化项可以定义为模型参数的绝对值之和，数学表达式如下：
$\text{L1 Regularization Term} = \lambda \|w\|_1 = \lambda \sum_{j=1}^n |w_j|$
其中：

( $w$ ) 是模型参数向量，( $w_j $) 是向量中的第 ( $j $) 个参数。
( $\lambda$ ) 是正则化系数，一个非负实数，用于控制正则化项的强度。
( $w\|_1$ ) 是向量 ( $w$ ) 的L1范数，即向量中所有元素绝对值的总和。

L1正则化项的作用

特征选择：L1正则化倾向于将不重要的特征参数压缩至零，从而实现自动特征选择，这在特征数量很多时特别有用。

防止过拟合：通过减少模型复杂度，L1正则化有助于防止模型在训练数据上过拟合。
提高模型解释性：由于L1正则化能够产生稀疏的模型参数，它可以帮助我们更容易地理解哪些特征对模型的预测结果有实质性的影响。

L1正则化与其他正则化技术的比较

L2正则化：L2正则化项是参数平方和的惩罚，它鼓励模型参数的值接近零但不完全为零，从而使得模型参数不会过于集中于某些特定的特征。L2正则化不会导致稀疏解。
Elastic Net：结合了L1和L2正则化的特点，通过调整两者的比例，可以在特征选择和参数收缩之间取得平衡。

L1正则化项在优化中的应用

在优化问题中，L1正则化项通常与损失函数结合使用，形成一个新的目标函数，如下所示：
$\text{Cost Function with L1 Regularization} = \text{Loss}(y, \hat{y}) + \lambda \|w\|_1$
其中， $\text{Loss}(y, \hat{y})$ 是模型的损失函数，$\hat{y} $是模型预测值。优化算法需要在最小化损失函数的同时，考虑L1正则化项的影响。

L1正则化项的引入增加了优化问题的复杂性，因为非线性和不可微的特性使得传统的优化方法不再适用。常用的求解方法包括坐标下降法、梯度下降法的变种（如近端梯度下降法）和最小角回归法等。

1-2：公式推导

先来回顾一下岭回归的代价函数，在原来标准线性回归代价函数上加上了一个带惩罚系数 λ 的 w 向量的L2-范数的平方：
$\operatorname{Cost}(w) = \sum_{i = 1}^N \left( y_i - w^Tx_i \right)^2 + \lambda\|w\|_{2}^{2}$
Lasso回归算法也同岭回归一样加上了正则项，只是改成加上了一个带惩罚系数 λ 的 w 向量的L1-范数作为惩罚项（L1-范数的含义为向量 w 每个元素绝对值的和），所以这种正则化方式也被称为L1正则化。
$\operatorname{Cost}(w) = \sum_{i = 1}^N \left( y_i - w^Tx_i \right)^2 + \lambda\|w\|_1$
同样是求使得代价函数最小时 w 的大小：
$\underset{w}{\operatorname{argmin}}\left(\sum_{i = 1}^N \left( y_i - w^Tx_i \right)^2 + \lambda\|w\|_1\right)$
由于加入的是向量的L1-范数，其中存在绝对值，导致其代价函数不是处处可导的，所以就没办法通过直接求导的方式来直接得到 w 的解析解。

下面介绍两种求解权重系数 w 的方法：坐标下降法（coordinate descent）、最小角回归法（Least Angle Regression，LARS）

坐标下降法

坐标下降法的核心与它的名称一样，就是沿着某一个坐标轴方向，通过一次一次的迭代更新权重系数的值，来渐渐逼近最优解。

具体步骤：

（1）初始化权重系数 w，例如初始化为零向量。

（2）遍历所有权重系数，依次将其中一个权重系数当作变量，其他权重系数固定为上一次计算的结果当作常量，求出当前条件下只有一个权重系数变量的情况下的最优解。

在第 k 次迭代时，更新权重系数的方法如下：
$\begin{matrix} w_m^k 表示第k次迭代，第m个权重系数 \\ w_1^k = \underset{w_1}{\operatorname{argmin}} \left( \operatorname{Cost}(w_1, w_2^{k-1}, \dots, w_{m-1}^{k-1}, w_m^{k-1}) \right) \\ w_2^k = \underset{w_2}{\operatorname{argmin}} \left( \operatorname{Cost}(w_1^{k}, w_2, \dots, w_{m-1}^{k-1}, w_m^{k-1}) \right) \\ \vdots \\ w_m^k = \underset{w_m}{\operatorname{argmin}} \left( \operatorname{Cost}(w_1^{k}, w_2^{k}, \dots, w_{m-1}^{k}, w_m) \right) \\ \end{matrix}$
（3）步骤（2）为一次完整迭代，当所有权重系数的变化不大或者到达最大迭代次数时，结束迭代。

preview

如上图所示，每次迭代固定其他的权重系数，只朝着其中一个坐标轴的方向更新，最后到达最优解。

最小角回归法

最小角回归法是一种特征选择的方法，计算出每个特征的相关性，通过数学公式逐渐计算出最优解。

具体步骤：

（1）初始化权重系数 w，例如初始化为零向量。

（2）初始化残差向量 residual 为目标标签向量 $y - Xw$ ，由于此时 w 为零向量，所以残差向量等于目标标签向量 y

（3）选择一个与残差向量相关性最大的特征向量 $x_i$ ，沿着特征向量 $x_i $的方向找到一组权重系数 w，出现另一个与残差向量相关性最大的特征向量 $x_j$ 使得新的残差向量 residual 与这两个特征向量的相关性相等（即残差向量等于这两个特征向量的角平分向量上），重新计算残差向量。

（4）重复步骤（3），继续找到一组权重系数 w，使得第三个与残差向量相关性最大的特征向量 $x_k$ 使得新的残差向量 residual 与这三个特征向量的相关性相等（即残差向量等于这三个特征向量的等角向量上），以此类推。

（5）当残差向量 residual 足够小或者所有特征向量都已被选择，结束迭代。

1-3：应用场景

在医学数据分析中，Lasso回归可以用于基因表达数据分析，通过选择与疾病相关的基因，帮助研究者识别生物标志物。例如，在癌症研究中，通过分析肿瘤样本的基因表达数据，Lasso回归可以帮助识别与癌症发生、发展和预后相关的基因。此外，在药物开发中，Lasso回归可以用于筛选影响药物效果的生物标志物，从而指导个性化治疗。
Lasso回归的稀疏性使得它特别适合处理高维数据，在医学研究中，数据往往具有成千上万的特征（如基因表达水平），而实际与疾病相关的因素可能只占少数。通过Lasso回归，可以在不牺牲预测准确性的同时，筛选出重要的特征，这对于理解疾病的生物学机制和开发新的治疗方法具有重要意义。

LASSO，全称Least absolute shrinkage and selection operator，是一种数据挖掘方法，即在常用的多元线性回归中，添加惩罚函数，不断压缩系数，从而达到精简模型的目的，以避免共线性和过拟合。当系数为0时，同时达到筛选变量的效果。（以下是一个不严谨的示意图）

所以，LASSO回归高效解决了筛选变量的难题：区别于传统的逐步回归stepwise前进、后退变量筛选方法，LASSO回归可以利用较少样本量，高效筛选较多变量。比如在基因组学、影像学、以及其他小样本分析中，LASSO回归都可以派上大用场。

二、教你用LASSO

参考资料

R语言glmnet包lasso回归中分类变量的处理

在第一部分小罗阐述了LASSO的原理，还讲了一堆公式，其实我们实际要用它并没有那么麻烦，已经有现成的工具可以供我们使用，我们只需要懂得如何把方法迁移到我们自己的数据即可。

巧妇难为无米之炊，我们开始演练之前第一步就应该先准备好数据集，建议大家先跟着小罗用同样的数据集，然后再去套自己的数据集。

2-1：数据集获取

后台回复“LASSO”即可。

2-2：glmnet包

R语言的glmnet包是由Jerome Friedman、Robert Tibshirani和Trevor Hastie共同开发的，用于拟合具有L1和L2正则化项的广义线性模型。这个包提供了一种高效的算法来估计回归模型的参数，特别是当特征数量接近或超过样本数量时。glmnet包支持多种类型的响应变量，包括连续的、二元的和计数数据。

主要特点

L1和L2正则化：glmnet实现了Lasso回归（L1正则化）和Ridge回归（L2正则化），以及它们的组合，即Elastic Net回归。
交叉验证：包内提供了交叉验证功能，用于选择最佳的正则化参数。
多种分布：glmnet支持多种概率分布，包括高斯分布（用于线性回归）、二项分布（用于逻辑回归）、泊松分布等。
并行计算：glmnet可以利用多核处理器进行并行计算，提高计算效率。
模型路径：glmnet能够计算和绘制模型路径，即随着正则化参数变化的模型系数路径。
稀疏性：glmnet生成的模型是稀疏的，即许多系数为零，有助于特征选择。

2-3：使用示例（R语言）

安装和加载glmnet包

library(glmnet)
library(foreign)
library("survival")
bc <- read.spss("/Users/luoxiaoluotongxue/Desktop/硕士课题进展记录/24年文献阅读/文献日推/24-06｜文献日推/24-06-14/Breast cancer survival agec.sav",
                use.value.labels=F, to.data.frame=T)
bc <- na.omit(bc)

这五行R代码用于加载必要的R包，读取SPSS格式的数据文件，并将其导入R环境中，同时排除含有缺失值的观测。

library(glmnet): 这行代码用于加载glmnet包，这是一个提供岭回归（Ridge Regression）和套索回归（Lasso Regression）函数的R包，常用于变量选择和高维数据分析。
library(foreign): 这行代码用于加载foreign包，该包提供了一些函数，用于读取和写入不同统计软件的数据格式，例如SPSS、SAS和Stata等。
library("survival"): 这行代码用于加载survival包，这是一个专门用于生存分析和生存数据的R包，提供了处理生存数据的工具和函数。
bc <- read.spss("/Users/luoxiaoluotongxue/Desktop/硕士课题进展记录/24年文献阅读/文献日推/24-06｜文献日推/24-06-14/Breast cancer survival agec.sav", use.value.labels=F, to.data.frame=T): 这行代码使用foreign包中的read.spss函数来读取位于指定路径的SPSS数据文件。参数use.value.labels=F指示函数不要使用SPSS数据中的值标签（value labels），而参数to.data.frame=T则指示函数将读取的数据转换为一个数据框（data.frame），这是R中用于存储表格数据的一种结构。
bc <- na.omit(bc): 这行代码使用na.omit函数从数据框bc中移除任何含有缺失值（NA）的行。na.omit函数返回一个没有缺失值的新数据框，只包含完全的观测。
综上所述，这些代码用于准备R环境，以便进行数据分析，特别是涉及到生存分析和需要使用glmnet包进行回归分析时。

数据分析

点开右侧工作区的bc（breast cancer），可以看到下图所示的数据。

age表示年龄
pathsize表示病理肿瘤大小（厘米）
lnpos表示腋窝淋巴结阳性
histgrad表示病理组织学等级
er表示雌激素受体状态
pr表示孕激素受体状态
status结局事件是否死亡
pathscat表示病理肿瘤大小类别（分组变量）
ln_yesno表示是否有淋巴结肿大
time是生存时间
后面的agec是自己设定的，不用管它

数据转换

bc$er<-as.factor(bc$er)
bc$pr<-as.factor(bc$pr)
bc$ln_yesno<-as.factor(bc$ln_yesno)
bc$histgrad<-as.factor(bc$histgrad)
bc$pathscat<-as.factor(bc$pathscat)

这五行R代码的目的是将数据框bc中的特定列转换为因子变量（factor variables）。

在R中，因子变量是一种用于存储分类数据的特殊类型，它可以表示有限数量的不同类别。将变量转换为因子对于后续的统计分析非常重要，因为它确保R正确地处理分类数据，特别是在进行回归分析或生成模型矩阵时。

在R中进行回归分析时，特别是当使用glmnet包进行岭回归或套索回归时，通常需要将分类变量转换为因子，以便模型可以正确地处理它们。如果分类变量不被转换为因子，R可能会错误地将它们视为连续变量，这可能会导致分析结果的错误。

构建模型

y<-bc$status
time<-bc$time

这两行R代码从数据框bc中选择了特定的列，并将它们分别赋值给变量y和time。在生存分析和cox回归模型中，y通常代表状态变量，表示观察结果（例如，事件发生或未发生），而time代表时间变量，表示从观察开始到事件发生或截止的时间长度。

y <- bc$status: 这行代码选择了bc数据框中的status列，并将其赋值给变量y。在生存分析中，status列通常是一个二进制变量，其中1表示事件发生（例如，死亡或疾病复发），0表示在观察期内未发生事件。y将作为生存分析或cox回归模型中的响应变量。
time <- bc$time: 这行代码选择了bc数据框中的time列，并将其赋值给变量time。在生存分析中，time列表示从观察开始到事件发生或截止的时间长度。如果是一个右截尾的数据（即，在观察期结束时事件尚未发生），则time列将表示从观察开始到截止时间的时间长度。time将作为生存分析或cox回归模型中的时间至事件发生变量。

剔除多余变量

data1<-bc[,-c(1,8,11,12)]

这行R代码从数据框bc中选择了一部分列，并创建了一个新的数据框data1，其中排除了指定的列。

bc[,-c(1,8,11,12)]: 这个表达式从bc数据框中选择所有行（，之前的部分是空的，表示选择所有行），并排除列索引为1、8、11和12的列。

，之前的部分表示行的选择，这里使用了空值，表示选择所有行。
，之后的部分表示列的选择，使用了-号，表示排除这些列。

具体来说，这行代码将排除bc数据框中的第1、第8、第11和第12列。在R中，列的索引是从1开始的，所以第1列是bc中的第一列，第8列是bc中的第八列，依此类推。

结果，data1将包含bc中除了第1、第8、第11和第12列之外的所有列。这种操作通常用于数据预处理，以便从分析中排除不必要或不需要的变量。

把分类变量变成哑变量矩阵形式

model_mat <-model.matrix(~ +er+pr+ln_yesno+histgrad+pathscat-1,data1)

这行R代码使用model.matrix函数从数据框data1中创建一个模型矩阵model_mat，该矩阵包含了分类变量（factor variables）的哑变量（dummy variables）表示。

model.matrix函数用于将模型公式和一个数据框作为输入，输出一个设计矩阵，这个矩阵包含了公式中所有变量的预测值，其中分类变量被转换为哑变量。

~: 这是模型公式的分隔符，左边是响应变量（在这个例子中是空的，因为没有指定），右边是解释变量。
+er+pr+ln_yesno+histgrad+pathscat: 这部分指定了模型中包含的变量。er、pr、ln_yesno、histgrad和pathscat是data1数据框中的列名（这些列是分类变量）。
-1: 这个符号告诉model.matrix不要包括截距项（常数项）在模型矩阵中。这意味着每个分类变量的哑变量表示将不包括参照类别（通常第一个类别，如果没有指定的话）。

model.matrix函数将自动处理分类变量，为每个类别创建一个哑变量，除了参照类别。参照类别被编码为所有哑变量都是0的情况。例如，如果一个分类变量有3个类别，那么将创建2个哑变量来表示这个变量。

这个模型矩阵model_mat可以用于统计模型中，如线性回归、逻辑回归等，因为它们期望输入数据是这种格式。
总之，这行代码的目的是将数据框中的分类变量转换为哑变量表示，以便进行统计分析。这种转换是必要的，因为统计模型不能直接处理分类数据，而是需要将它们转换为数值型数据。

知识点补充：哑变量

哑变量（dummy variable），也称为指示变量（indicator variable），是在统计学和机器学习中使用的一种特殊类型的变量，用于表示分类数据中的各个类别。哑变量通常用于回归分析和分类问题中，以便模型可以理解和解释分类数据。

哑变量的基本思想是，将一个具有多个类别的分类变量转换为多个二进制（0和1）变量，每个类别对应一个二进制变量。这样，原始的分类变量就可以表示为一系列的哑变量，每个哑变量代表一个可能的类别。

再举一个例子，假设有一个分类变量“颜色”，它有三个类别：红色、绿色和蓝色。为了在统计分析中使用这个变量，我们可以创建两个哑变量：

颜色_红色（1表示红色，0表示非红色）
颜色_绿色（1表示绿色，0表示非绿色）

在这种情况下，蓝色类别不需要单独的哑变量，因为当红色和绿色的哑变量都为0时，就隐含地表示了蓝色。

哑变量的使用有以下几个目的：

允许分类数据在回归模型中使用：回归模型通常适用于数值数据，哑变量允许我们将分类数据转换为数值格式，以便模型可以正确地处理它们。
捕捉类别效应：通过使用哑变量，模型可以捕捉到不同类别之间的差异，并估计每个类别的效应。
避免多重共线性：在回归分析中，如果直接使用分类变量的原始编码（如1, 2, 3等），可能会导致多重共线性问题。哑变量通过为每个类别创建一个独立的变量来避免这个问题。
简化模型解释：哑变量使得模型的结果更易于解释，因为每个哑变量的系数可以直接解释为对应类别的效应。
在R中，model.matrix函数可以自动创建哑变量，这通常是在拟合统计模型之前的步骤。通过这种方式，分类数据可以有效地纳入到统计分析中。

重新组合成数据

x<-as.matrix(data.frame(age=data1$age,
                        pathsize=data1$pathsize,lnpos=data1$lnpos,model_mat))

这行R代码的目的是将几个变量从数据框data1中选取出来，并重新组合成一个矩阵x。这个新的矩阵将包含age、pathsize、lnpos以及model_mat这几个变量，这些变量会用于后续的统计分析或建模。

age=data1$age: 这里从data1数据框中选取了名为age的列。
pathsize=data1$pathsize: 这里从data1数据框中选取了名为pathsize的列。
lnpos=data1$lnpos: 这里从data1数据框中选取了名为lnpos的列。
model_mat: 这个矩阵在前面已经定义了，所以直接加上就可以，你从工作环境的列数对比就可以发现了。

这些选定的变量被传递给data.frame函数，该函数将它们组合成一个数据框。然后，as.matrix函数将这个数据框转换为一个矩阵。

整个表达式可以理解为：创建一个数据框，其中包含data1中的age、pathsize和lnpos列，以及预先定义的model_mat矩阵，然后将这个数据框转换为矩阵x。

这种操作的常见场景是在进行统计分析或机器学习模型训练之前，需要从原始数据中提取相关的特征变量。通过重新组合数据，可以确保模型使用了所有感兴趣的变量，并且它们是以正确的格式（矩阵）呈现的。

交叉验证

set.seed(123)
cv.fit <- cv.glmnet(x,Surv(time,y),family="cox", maxit = 1000)
plot(cv.fit)

这三行R代码用于设置随机种子，以确保随机过程的可重复性，并使用交叉验证来选择最佳的λ值，然后绘制交叉验证的结果。

set.seed(123): 这行代码设置了一个随机种子，值为123。在R中，set.seed函数用于设置随机数生成器的种子，从而确保每次运行代码时生成相同的随机数序列。这对于需要可重复性的实验或分析非常有用。例如，如果你想要重复一个分析并得到相同的结果，你可以设置一个种子，然后每次运行代码时使用相同的种子值。
cv.fit <- cv.glmnet(x,Surv(time,y),family="cox", maxit = 1000): 这行代码使用cv.glmnet函数来执行交叉验证，以选择最佳的λ值。cv.glmnet是glmnet包中的一个函数，它结合了glmnet和交叉验证功能。参数解释如下：
- x: 是一个矩阵，包含了所有的自变量或解释变量。
- Surv(time,y): 这是一个生存分析的响应变量，其中time是生存时间，y是状态变量（1表示事件发生，0表示未发生）。
- family="cox": 指定模型使用的分布族。在这里，它告诉cv.glmnet使用Cox比例风险模型。
- maxit = 1000: 指定最大迭代次数。
  cv.fit是一个包含交叉验证结果的对象，它包含了最佳λ值、每个λ值对应的模型系数等。
plot(cv.fit): 这行代码使用plot函数来绘制cv.fit对象的交叉验证结果。这将显示不同λ值下模型的预测性能，通常包括交叉验证的误差曲线和λ路径图。这些图表有助于理解和选择最佳的λ值。

拟合Cox比例风险模型，并绘制模型的预测性能图

fit <- glmnet(x, Surv(time,y), family =  "cox", maxit = 1000)
plot(fit)

fit <- glmnet(x, Surv(time,y), family = "cox", maxit = 1000): 这行代码调用了glmnet函数来拟合一个广义线性模型（GLM），这里使用的是Cox比例风险模型（family="cox"），适用于生存分析。参数解释如下：
- x: 是一个矩阵，包含了所有的自变量或解释变量。
- Surv(time,y): 这是一个生存分析的响应变量，其中time是生存时间，y是状态变量（1表示事件发生，0表示未发生）。
- family="cox": 指定模型使用的分布族。在这里，它告诉glmnet拟合一个Cox比例风险模型。
- maxit = 1000: 指定最大迭代次数。
  fit是一个包含模型结果的对象，它包含了模型系数、λ值、残差等。
plot(fit): 这行代码使用plot函数来绘制fit对象的预测性能图。这将显示不同λ值下模型的预测性能，通常包括交叉验证的误差曲线和λ路径图。这些图表有助于理解和选择最佳的λ值。

总的来说，这两行代码用于拟合一个Cox比例风险模型，并绘制模型的预测性能图，以便更好地理解和选择最佳的λ值。

查看和提取系数

Coefficients <- coef(fit, s = cv.fit$lambda.min)
Active.Index <- which(Coefficients != 0)
Active.Coefficients <- Coefficients[Active.Index]
Active.Index
Active.Coefficients

这几行R代码是从拟合的Cox比例风险模型fit中提取活动（非零）系数，并返回活动系数的索引和值。

Coefficients <- coef(fit, s = cv.fit$lambda.min): 这行代码使用coef函数从fit对象中提取系数。fit是一个glmnet包中的对象，它包含了拟合的Cox模型。s = cv.fit$lambda.min参数指定使用交叉验证结果中的最小λ值（cv.fit$lambda.min是cv.fit对象中的一个元素，它包含了交叉验证过程中每个λ值对应的系数）。
coef函数返回一个系数向量，其中包含了模型中每个变量的系数。
Active.Index <- which(Coefficients != 0): 这行代码使用which函数找到系数向量中不等于0的系数的索引。这些索引对应于模型中具有非零系数的变量。
Active.Coefficients <- Coefficients[Active.Index]: 这行代码使用索引Active.Index从系数向量Coefficients中提取活动系数。这些活动系数是模型中显著的变量，即那些对生存时间有显著影响的变量。
Active.Index: 这是一个赋值表达式，它将活动系数的索引存储在一个新的变量Active.Index中。
Active.Coefficients: 这是一个赋值表达式，它将活动系数存储在一个新的变量Active.Coefficients中。

总的来说，这些代码行用于从拟合的Cox模型中提取活动系数，并返回这些系数的索引和值。这些活动系数可以用来理解哪些变量对生存时间有显著影响，这对于后续的分析或解释模型结果非常有用。

拟合二元逻辑回归模型

fit1 = glmnet(x, y, family = "binomial")
plot(fit1, xvar = "dev", label = TRUE)

这两行R代码用于拟合一个二元逻辑回归模型，并绘制模型性能的图。

fit1 = glmnet(x, y, family = "binomial"): 这行代码调用了glmnet函数来拟合一个二元逻辑回归模型。glmnet是一个R包，提供了用于拟合岭回归（Ridge Regression）和套索回归（Lasso Regression）的函数，也支持二元逻辑回归模型。参数解释如下：
- x: 是一个矩阵，包含了所有的自变量或解释变量。
- y: 是一个矩阵或向量，包含了因变量或响应变量。在这个例子中，y应该是一个二元响应变量。
- family="binomial": 指定模型使用的分布族。在这里，它告诉glmnet拟合一个二元逻辑回归模型。
  fit1是一个包含模型结果的对象，它包含了模型系数、λ值、残差等。
plot(fit1, xvar = "dev", label = TRUE): 这行代码使用plot函数来绘制fit1对象的预测性能图。xvar = "dev"参数指定x轴的变量为模型开发函数（deviance），这通常用于评估逻辑回归模型的性能。label = TRUE参数指示添加图例。
plot函数会绘制模型开发函数随λ值变化的曲线图，以及不同λ值下模型的预测性能。这些图表有助于理解和选择最佳的λ值。

总的来说，这两行代码用于拟合一个二元逻辑回归模型，并绘制模型性能的图，以便更好地理解和选择最佳的λ值。

现在替换X轴变量为lambda。

plot(fit1, xvar="lambda", label=TRUE)

交叉验证

set.seed(999)
cvfit=cv.glmnet(x,y, family = "binomial")
plot(cvfit)

这三行R代码用于设置随机种子以确保随机过程的可重复性，使用交叉验证来选择最佳的λ值，并绘制交叉验证的结果。

set.seed(999): 这行代码设置了一个随机种子，值为999。在R中，set.seed函数用于设置随机数生成器的种子，从而确保每次运行代码时生成相同的随机数序列。这对于需要可重复性的实验或分析非常有用。例如，如果你想要重复一个分析并得到相同的结果，你可以设置一个种子，然后每次运行代码时使用相同的种子值。
cvfit=cv.glmnet(x,y, family = "binomial"): 这行代码使用cv.glmnet函数来执行交叉验证，以选择最佳的λ值。cv.glmnet是glmnet包中的一个函数，它结合了glmnet和交叉验证功能。参数解释如下：
- x: 是一个矩阵，包含了所有的自变量或解释变量。
- y: 是一个矩阵或向量，包含了因变量或响应变量。在这个例子中，y应该是一个二元响应变量。
- family="binomial": 指定模型使用的分布族。在这里，它告诉cv.glmnet使用二元逻辑回归模型。
  cvfit是一个包含交叉验证结果的对象，它包含了最佳λ值、每个λ值对应的模型系数等。
plot(cvfit): 这行代码使用plot函数来绘制cvfit对象的交叉验证结果。这将显示不同λ值下模型的预测性能，通常包括交叉验证的误差曲线和λ路径图。这些图表有助于理解和选择最佳的λ值。

总的来说，这三行代码用于设置随机种子以确保实验的可重复性，使用交叉验证来选择最佳的λ值，并绘制交叉验证的结果以帮助选择最佳的λ值。

求出最小值

cvfit$lambda.min#求出最小值
cvfit$lambda.1se#求出最小值一个标准误的λ值

在glmnet包中，cvfit对象包含了交叉验证过程中得到的多个λ值及其对应的性能指标。这两个属性lambda.min和lambda.1se提供了特定的λ值，它们在模型选择和超参数优化中非常有用。

cvfit$lambda.min:
- 这表示在交叉验证过程中，模型性能最接近于全模型（即没有L1惩罚）的λ值。
- 在glmnet中，这通常对应于具有最小交叉验证误差（CV误差）的λ值。
- 选择lambda.min对应的模型可以平衡模型的预测性能和复杂度，通常被认为是一个好的模型选择。
cvfit$lambda.1se:
- 这表示在交叉验证过程中，模型性能与lambda.min对应的模型性能相差一个标准误差的λ值。
- 选择lambda.1se对应的模型可以提供比lambda.min模型稍差但更稳健的模型选择。
- 在实践中，lambda.1se模型通常被认为是一个比lambda.min模型更保守的选择，因为它包含了更多的模型复杂度，但仍然保持了与lambda.min模型相近的性能。

在实际应用中，这两个值通常用于模型选择和超参数优化。选择lambda.min模型可以得到一个性能较好的模型，而选择lambda.1se模型可以得到一个更稳健的选择。这两个值可以帮助你根据模型的复杂度和预测性能来做出决策。

不同λ值下的模型系数

coef1<-coef(cvfit, s = "lambda.min")
coef2<-coef(cvfit, s = "lambda.1se")
coef1
coef2

这几行R代码用于从交叉验证后的cvfit对象中提取两个不同λ值下的模型系数，并分别赋值给变量coef1和coef2。

coef1 <- coef(cvfit, s = "lambda.min"):
- 这行代码使用coef函数从cvfit对象中提取模型系数。
- 参数s = "lambda.min"告诉coef函数使用cvfit对象中的最小λ值对应的系数。
- coef1是一个包含lambda.min对应的模型系数的向量。
coef2 <- coef(cvfit, s = "lambda.1se"):
- 这行代码使用coef函数从cvfit对象中提取模型系数。
- 参数s = "lambda.1se"告诉coef函数使用cvfit对象中与lambda.min相差一个标准误差的λ值对应的系数。
- coef2是一个包含lambda.1se对应的模型系数的向量。

这些系数可以用来理解不同模型复杂度下的变量重要性，以及选择最合适的模型。通常，lambda.min对应的模型系数被认为是最优的，因为它在预测性能和模型复杂度之间达到了一个平衡。而lambda.1se对应的模型系数可以作为一个保守的选择，因为它提供了更多的模型复杂度，但仍然保持了与lambda.min模型相近的性能。

估计大家看到这里还是有点懵，感觉懂一点了，但是又不太多，所以小罗贴心的准备了下一章节——结合一篇文献，带你加深印象！！

三、LASSO在SCI文章中的应用

在这一部分，小罗挑选了一篇文献，带着大家进一步熟悉LASSO回归在文献中如何运用。

3-1：文献概述

这篇文章是一项关于甲胎蛋白阴性肝细胞癌（AFP-NHCC）患者在接受非手术治疗后的预后评估的研究。

研究团队开发并验证了一种基于LASSO Cox回归的新列线图模型，用于预测这些患者的预后。

背景：AFP-NHCC具有特殊的临床病理特征和预后，传统上依赖于超声成像和AFP进行筛查，但存在局限性，特别是AFP水平正常的患者可能会错失早期诊断和治疗机会。
方法：研究包括410名AFP阴性的HCC患者作为主要队列，148名AFP-NHCC患者作为独立验证队列。使用LASSO Cox回归分析确定独立预后因素，并构建列线图模型（nomogram1），同时使用前向逐步Cox回归构建另一个模型（nomogram2）。通过一致性指数（C-index）、曲线下面积（AUC）、校准曲线和决策曲线分析（DCA）评估模型性能。
结果：nomogram1的C-index为0.708，优于nomogram2和其他传统模型。在验证队列中，nomogram1仍然显示出良好的区分能力。此外，根据nomogram1的总分，患者被分为低风险、中风险和高风险三个不同的预后组。
结论：基于LASSO Cox回归的新型列线图为AFP-NHCC患者提供了更准确和有用的预后预测，有助于这些患者进行个性化的预后评估。

研究还讨论了AFP-NHCC的特殊性，以及现有分期系统在预测AFP-NHCC患者预后方面的局限性。研究结果表明，新开发的列线图模型在预测AFP-NHCC患者预后方面比传统模型更准确。

3-2：基线临床病理特征

表1展示了训练队列和验证队列中AFP-NHCC患者的基线临床病理特征。

以下是对表中数据的分析：

性别和年龄：
- 在两个队列中，大多数患者为男性（训练队列79.8%，验证队列83.8%）。
- 多数患者年龄超过50岁（训练队列77.8%，验证队列83.1%）。
肿瘤大小：
- 训练队列中56.1%的患者肿瘤小于3厘米，而验证队列中这一比例为56.8%。
- 肿瘤在3-5厘米之间的患者分别占19.8%和25.0%。
- 肿瘤大于或等于5厘米的患者分别占24.1%和18.2%。
肿瘤多重度：
- 训练队列中64.9%的患者有单个肿瘤，而验证队列中这一比例为52.0%。
- 有多个肿瘤的患者分别占35.1%和48.0%，在验证队列中比例更高，这在统计上有显著差异（P=0.006）。
肝硬化：
- 绝大多数患者存在肝硬化（训练队列92.0%，验证队列91.2%）。
门静脉癌栓（PVTT）：
- PVTT在训练队列中出现在35.4%的患者中，在验证队列中为27.0%。
腹水：
- 存在腹水的患者分别占39.5%和35.8%。
乙肝表面抗原（HBsAg）：
- 大多数患者HBsAg呈阳性（训练队列78.5%，验证队列81.1%）。
中性粒细胞与淋巴细胞比率（NLR）：
- NLR小于5的患者占大多数（训练队列83.9%，验证队列89.2%）。
血红蛋白（HGB）：
- HGB大于或等于120 g/L的患者占多数（训练队列61.2%，验证队列60.8%）。
血小板（PLT）：
- PLT大于或等于100*10^9/L的患者占46.6%和45.9%。
肌酐（CR）：
- CR小于111μmol/L的患者占绝大多数（训练队列95.4%，验证队列95.3%）。
肝功能指标（ALT, AST, TBIL, ALB等）：
- 多数患者的肝功能指标在正常范围内。
肿瘤分期（BCLC, TNM）和Child-Pugh评分：
- 大多数患者的肿瘤分期为0-B（训练队列78.5%，验证队列80.4%）。
- Child-Pugh评分大多数为A或B（训练队列87.1%，验证队列92.6%）。

表1中的数据提供了两个独立队列中患者的详细基线特征，这些特征对于评估预后模型的适用性和准确性至关重要。不同特征之间的比较有助于识别可能影响预后的变量。

3-3：LASSO Cox回归分析

使用LASSO Cox回归分析在AFP-NHCC患者中识别预后风险因素的过程和结果。

以下是对图中信息的分析：

LASSO系数轮廓（Fig. 1a）：
- 这部分展示了23个潜在风险因素的LASSO系数轮廓。LASSO（最小绝对值收缩和选择算子）是一种变量选择方法，它通过惩罚回归系数的绝对值来实现变量的自动选择和系数的缩减。
- 在LASSO系数轮廓图中，水平线表示惩罚系数，而垂直线表示不同变量的系数大小。当惩罚系数增加时，一些系数会缩减至零，这意味着这些变量从模型中被排除。
五个风险因素的选择（Fig. 1b）：
- 通过LASSO Cox回归分析，最终选择了五个风险因素：PVTT（门静脉癌栓）、腹水、HGB（血红蛋白）、γ-GGT（γ-谷氨酰转肽酶）和CRP（C反应蛋白）。
- 这两个虚线表示的是最小标准和1标准误差（s.e.）标准下的最佳分数。在最小标准下，包括性别、年龄、肿瘤大小、肿瘤数量、肝硬化、PVTT、腹水、HBV（乙型肝炎病毒）、HGB、CR（肌酐）、AST（天门冬氨酸转氨酶）、ALB（白蛋白）、LDH（乳酸脱氢酶）、γ-GGT、CA199（糖类抗原199）和CRP等在内的多个变量被考虑。
- 在1-s.e.标准下，最终确定的五个变量被认为对预后有显著影响，这些变量的系数没有被缩减至零，表明它们在模型中是重要的预测因子。
模型选择和变量重要性：
- LASSO方法通过调整惩罚系数来选择模型，这个过程可以通过交叉验证来确定最佳的惩罚系数，从而选择最佳的模型。
- 在交叉验证过程中，模型的性能通过一致性指数（C-index）来评估，C-index越高，模型的预测能力越好。
- 最终选择的五个变量可能因为它们与总体生存（OS）有显著的相关性，并且在统计上通过了显著性检验。