预后建模绕不开的lasso cox回归
回归我们并不陌生，线性回归和最小二乘法，逻辑回归和最大似然法，这些都是我们耳熟能详的事物，在生物信息学中的应用也比较广泛, 回归中经常出现两类问题，欠拟合和过拟合。
对于欠拟合，简单而言就是我们考虑的少了，一般通过在回归模型中增加自变量或者扩大样本数量来解决；对于过拟合，简单而言就是考虑的太多了，模型过于复杂了，这时候可以对已有的自变量进行筛选，在代价函数中增加惩罚项来限制模型的复杂度，增加的惩罚项我们称之为正则化，正则化常用的有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 02. 交叉验证
通过交叉验证，选择最佳的λ值。在选择λ值时，我们需要指定评价指标，就是根据评价指标的值来选择最佳模型和最佳λ值，对应的是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 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函数可以显示自变量的回归系数，可以看到很多自变量的回归系数都是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·