欢迎关注”生信修炼手册”!

回归我们并不陌生，线性回归和最小二乘法，逻辑回归和最大似然法，这些都是我们耳熟能详的事物，在生物信息学中的应用也比较广泛, 回归中经常出现两类问题，欠拟合和过拟合。

对于欠拟合，简单而言就是我们考虑的少了，一般通过在回归模型中增加自变量或者扩大样本数量来解决；对于过拟合，简单而言就是考虑的太多了，模型过于复杂了，这时候可以对已有的自变量进行筛选，在代价函数中增加惩罚项来限制模型的复杂度，增加的惩罚项我们称之为正则化，正则化常用的有L1正则化和L2正则化，

所谓正则化Regularization, 指的是在回归模型代价函数后面添加一个约束项， 在线性回归模型中，有两种不同的正则化项

1. 所有参数绝对值之和，即L1范数，对应的回归方法叫做Lasso回归

2. 所有参数的平方和，即L2范数，对应的回归方法叫做Ridge回归，岭回归

lasso回归对应的代价函数如下

岭回归对应的代价函数如下

红框标记的就是正则项，需要注意的是，正则项中的回归系数为每个自变量对应的回归系数，不包含回归常数项。

在预后建模的文章中，我们需要针对多个marker基因的表达量汇总形成一个指标，使用该指标来作为最终的maker, 而这个指标在文章中被称之为各种risk score, 比如NAD+基因的预后模型，构建的maker就叫做NPRS, 全称的解释如下

The NAD+ metabolism-related prognostic risk score (NPRS) of each sample was calculated using the formula: NPRS = ΣExp (mRNAί) × Coefficient (mRNAί)

所以各种的预后建模，其实都是lasso回归技术在生物信息学领域的应用。注意观察上述的Lasso回归代价函数，，可以看到有一个未知数λ， 这个参数是一个惩罚项的系数，数值越大，惩罚项对应的影响就越大，我们求解的目标是代价函数值最小，λ = 0时，惩罚项失去意义，代价函数变成了普通的线性回归，而λ过大，惩罚项的影响被放的过大，过小时，惩罚项又失去了原本的意义，所以使用lasso回归，第一个问题是设置合理的λ 值。

这个λ 值 如何设置呢？最简单的办法是找到两个队列，训练集和验证集，适应一系列的λ值对训练集进行建模，观察模型在验证集上的表现，然后选择在验证集上表现最佳模型的λ值，当没有额外的验证集时，就只能通过交叉验证的方式将数据集人工划分为训练集和验证集，然后进行分析。在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/

正则项本身只是一个代价函数中的添加项，所以其应用范围不仅局限于线性回归，逻辑回归，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
# 自变量数据，每一行表示一个患者，每一列表示一个自变量
> 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
# 因变量数据，生存数据的因变量为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. 交叉验证

通过交叉验证，选择最佳的λ值。在选择λ值时，我们需要指定评价指标，就是根据评价指标的值来选择最佳模型和最佳λ值，对应的是typpe.measure参数，对于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 λ值，纵坐标为每个λ值对应的评价指标，用error bar的形式展现了多个模型评价指标的均值+标准误，可以看到在图中有两条垂直的虚线，左边的虚线对应评价指标最佳的λ值，即lambda.min, c-index值越大越好，deviance值越小越好；右边的虚线表示评价指标在最佳值1个标准误范围的模型的λ值，即lambda.1se, 通过以下方式可以提取对应的值

> cvfit$lambda.min
[1] 0.01749823
> cvfit$lambda.1se
[1] 0.04868986

通过print函数可以看到交叉验证的关键信息

> print(cvfit)Call:  cv.glmnet(x = x, y = y, type.measure = "deviance", nfolds = 10, family = "cox")Measure: Partial Likelihood DevianceLambda Index Measure      SE Nonzero
min 0.01750    29   13.08 0.06221      15
1se 0.04869    18   13.14 0.05369      10

通过coef函数可以显示自变量的回归系数，可以看到很多自变量的回归系数都是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提取回归系数，我们就得到了最终的回归模型。

·end·

—如果喜欢，快分享给你的朋友们吧—

原创不易，欢迎收藏，点赞，转发！生信知识浩瀚如海，在生信学习的道路上，让我们一起并肩作战！

本公众号深耕耘生信领域多年，具有丰富的数据分析经验，致力于提供真正有价值的数据分析服务，擅长个性化分析，欢迎有需要的老师和同学前来咨询。

更多精彩

• KEGG数据库，除了pathway你还知道哪些

• 全网最完整的circos中文教程

• DNA甲基化数据分析专题

• 突变检测数据分析专题

• mRNA数据分析专题

• lncRNA数据分析专题

• circRNA数据分析专题

• miRNA数据分析专题

• 单细胞转录组数据分析专题

• chip_seq数据分析专题

• Hi-C数据分析专题

• HLA数据分析专题

• TCGA肿瘤数据分析专题

• 基因组组装数据分析专题

• CNV数据分析专题

• GWAS数据分析专题

• 机器学习专题

• 2018年推文合集

• 2019年推文合集

• 2020推文合集

写在最后

转发本文至朋友圈，后台私信截图即可加入生信交流群，和小伙伴一起学习交流。

扫描下方二维码，关注我们，解锁更多精彩内容！

一个只分享干货的

生信公众号