

线性回归（6）—— Lasso的推导与实战

小明的数据分析笔记本

2021-10-18

导读：终于说到Lasso了，Lasso在生信领域应该是一个老生常谈的算法，无论懂不懂，一波单因素，Lasso，多因

终于说到Lasso了，Lasso在生信领域应该是一个老生常谈的算法，无论懂不懂，一波单因素，Lasso，多因素就往上怼。这是生信套路市场化的结果。这一期，我将从矩阵推导，与岭回归的区别，偏差与方差，Lasso的意义，还有代码实战五个方面来聊这个事儿。看来这应该是篇长文，而且比较晦涩，浏览量应该不会多，我想能看到最后的都是好朋友。

1.Lasso的矩阵推导

和岭回归一样，Lasso是被创造来作用于多重共线性问题的算法(注意用词，不是解决)，不过Lasso使用的是系数的L1范式（L1范式则是系数的绝对值）乘以正则化系数，所以Lasso的损失函数表达式为：

对损失函数进行求导：

所以，这里的问题又回到了OLS中XTX的逆必须存在，回顾下岭回归最后的结果：

在岭回归中，我们通过正则化系数α能够向方阵XTX加上一个单位矩阵，以此来防止方阵XTX的行列式为0，而现在L2范式所带的正则项α在求导后，并不带有w这个项，所以它无法对XTX造成任何影响。也就是说，Lasso无法解决特征之间“精确相关”的问题。当我们使用OLS求解线性回归的时候，如果线性回归无解或者报除0错误，换Lasso不能解决任何问题。

岭回归 vs Lasso

岭回归可以解决特征间的精确相关关系导致的最小二乘法无法使用的问题，而Lasso不行。

幸运的是，在现实中我们其实会比较少遇到“精确相关”的多重共线性问题，大部分多重共线性问题应该是“高度关“，而如果我们假设方阵XTX的逆是一定存在的，那我们可以有：

通过增大α，我们可以为w的计算增加一个负项，从而限制参数估计中w的大小，而防止多重共线性引起的参数被估计过大导致模型失准的问题。Lasso不是从根本上解决多重共线性问题，而是限制多重共线性带来的影响。何况，这还是在我们假设所有的系数都为正的情况下，假设系数w无法为正，则很有可能我们需要将我们的正则项参数α设定为负，因此α可以取负数，并且负数越大，对共线性的限制也越大。

推导完成了，但是问题也来了，Lasso是如何把系数惩罚到0了，下一部分会说到。

2.Lasso vs Ridge

假设w是二维的，那么：

上面分别是Lasso和岭回归的几何描述。蓝色区域是对w的约束区域，比如|β1|+|β2|≤s是Lasso的限制条件，β1^2+β2^2≤s是岭回归的约束区域。LASSO由于是绝对值形式，其约束条件更为"尖锐"。红色的圆圈是RSS的等高线，我们期望找到既满足在限制范围内又满足RSS最小的β1和β2的取值。

在岭回归中，是跟圆的交点，等高线很难交点到（0，β2）或（β1，0）的点，则系数基本不会等于0。但是对于Lasso就不太一样了，等高线很容易就与四个角相交，这样另外一个系数必然等于0。所以，Lasso就起到了降维的作用。

岭回归 vs Lasso

Lasso可以做特征选择，岭回归不行。

3.偏差与方差

写文章的时候，会写到，为了防止模型的过拟合，我们进一步使用了Lasso回归.......

那么，为什么会过拟合呢？因为一个数据集中不可能是完美的数据的，比如肯定会有噪音存在，而且，我们取得数据集会不会集中在总体数据的冰山一角。所以，我们对良好的模型有个要求，就是其泛化能力是ok的，即泛化误差较小。那么什么是泛化误差呢，前面其实也说过，我们一般用均方误差MSE表示，这里则细细地剖析下泛化误差的构成和对建模过程中的启发。

“偏差-方差分解” (bias􀁭variance decomposition)也是解释学习算法泛化性能的一种重要工具

以上是周志华老师在《机器学习》一书中的描述，我觉得比我讲得好，所以直接截取书中片段来阐明这个问题。

泛化误差被分解为偏差、方差与噪声之和。我们来分别看看这几个名词：

偏差(bias): 偏差度量了模型的预测均值与真实结果的偏离程度，即刻画了学习算法本身的拟合能力。
方差(variance): 方差度量了同样大小的训练集的变动所导致的学习性能的变化，即刻画了数据扰动所造成的影响。
噪音(noise): 噪音表达了当前任务上任何学习算法所能达到期望泛化误差MSE的下界，即刻画了学习问题本身的难度。

所以，偏差-方差分解说明，泛化性能是由学习算法的能力，数据的充分性以及学习任务本身的难度共同决定的。

当给定学习任务(noise固定)，为了取得好的泛化性能，则需要偏差较小，这样能够充分拟合数据，方差较小，这样使得数据扰动产生的影响小。

我们从下面的靶心图来对偏差和方差有个直观的感受：

假设红色的靶心区域是学习算法完美的正确预测值，蓝色点为训练数据集所训练出的模型对样本的预测值，当我们从靶心逐渐往外移动时，预测效果逐渐变差。

从上面的图片中很容易可以看到，左边一列的蓝色点比较集中，右边一列的蓝色点比较分散，它们描述的是方差的两种情况。比较集中的属于方差比较小，比较分散的属于方差比较大的情况。

我们再从蓝色点与红色靶心区域的位置关系来看，靠近红色靶心的属于偏差较小的情况，远离靶心的属于偏差较大的情况。

Question：从上面的图中可以看出，模型不稳定时会出现偏差小、方差大的情况，那么偏差和方差作为两种度量方式有什么区别呢？

Answer：Bias的对象是单个模型，是期望输出与真实标记的差别。它描述了模型对本训练集的拟合程度。Variance的对象是多个模型，是相同分布的不同数据集训练出模型的输出值之间的差异。它刻画的是数据扰动对模型的影响。

但是，偏差和方差是有冲突的，这称为偏差-方差窘境(bias-variance dilemma). 下图给出了一个示意图。给定学习任务，假定我们能控制学习算法的训练程度，则在训练不足时，学习器的拟合能力不够强，训练数据的扰动不足以使学习器产生显著变化，此时偏差主导了泛化错误率；随着训练程度的加深，学习器的拟合能力逐渐增强，训练数据发生的扰动渐渐能被学习器学到，方差逐渐主导了泛化错误率；在训练程度充足后，学习器的拟合能力已经非常强，训练数据发生的轻微扰动都会导致学习器发生显著变化，若训练数据自身的、非全局的特性被学习器学到了，则将发生过拟合。

4.Lasso的意义

先起个序，在文明高度发达的21世纪，互联网成为了我们获取信息的主要途径，我们拿起手机，想获得什么样的信息，简直不能太方便。那么虽然获取信息的总量多了太多，但每个时代的大师好像都是那些寥寥数人。这说明有效的信息量并没有怎么变过。

这也是Lasso的世界观——稀疏性假设。

简单来说，我们认为，尽管世界如此复杂，但有用的信息却非常有限。套用到常见的统计学模型中，假如，我们考虑一个线性回归模型，有一个因变量Y，但有成百上千的自变量X。我们假设，只有有限个X的回归系数不为0，但其余的都是0。也就是说他们跟Y并没有啥子特别显著的关系。找到其中重要的X，对我们理解数据有重要的意义。对应前文的例子，生物学家想要研究基因对于某类疾病的影响，面对上万个可能的基因，生物学家们倾向于认为只有其中的一小部分对于该类疾病有着显著的影响；而为了预测消费者对于电影和书籍的喜好，线上电影和书店也倾向于认为一个消费者的喜好可以从少量的评分数据中得到。而“稀疏性假设”就是对于人们这种倾向的具体体现，举个例子，在分析基于和疾病的关系时，对于我们放入的10000个可能的基因，我们认为这10000个基因在回归模型中的回归系数只有少量的不为0。

说到这里，可能大家想到了WGCNA，是的，Lasso的世界观和WGCNA的世界观是一致的。WGCNA认为，在诺大的生物网络中，只有少数的节点起到了枢纽作用，这便是无尺度网络。我想，Lasso和WGCNA现实中见面应该是特别不错的朋友。

运用到管理学中，有种理论叫做“奥卡姆剃刀”，意思是简约之法则，是由14世纪逻辑学家奥卡姆的威廉提出的一个解决问题的法则，他在《箴言书注》2卷15题说“切勿浪费较多东西，去做‘用较少的东西，同样可以做好的事情’。”换一种说法，如果关于同一个问题有许多种理论，每一种都能作出同样准确的预言，那么应该挑选其中使用假定最少的。尽管越复杂的方法通常能做出越好的预言，但是在不考虑预言能力（即结果大致相同）的情况下，假设越少越好。

道德经中也说：万物之始，大道至简，衍化至繁。

到这里奇怪的发现，Lasso属于数学领域，WGCNA属于生物学领域，奥卡姆剃刀属于管理学领域，道德经属于哲学领域。这些不同领域的前辈们竟然得出了类似的结论，真是殊途同归，条条大路通罗马啊。

5.Lasso的实战

这里本来想详解下glmnet函数的，但内容着实不少，所以决定这期glmnet不作为重点，后面我们再聊聊glmnet的细节。

应用Lasso也比较简单，相比较于岭回归，只需要把glmnet的alpha=0改为alpha=1即可。

library(ElemStatLearn) #contains the datalibrary(car) #package to calculate Variance Inflation Factorlibrary(glmnet) #allows ridge regression, LASSO and elastic nedata(prostate)str(prostate)##载入包和数据

train代表训练集和测试集；lpsa在这里作为响应y，1：8列数据在这里作为预测变量x

prostate$gleason <- ifelse(prostate$gleason == 6, 0, 1) ##divide gleason score into <=6 and >6 groupstable(prostate$gleason)##将gleason评分变为二分类变量 <=6 和 >6
train <- subset(prostate, train == TRUE)[,1:9]test = subset(prostate, train==FALSE)[,1:9]##将数据集分为训练集和测试集

x <- as.matrix(train[, 1:8])y <- train[, 9]fit <- glmnet(x, y, family = "gaussian", alpha = 1)print(fit)

列 Df 是自由度，代表了非零的线性模型拟合系数的个数。列 %Dev 代表了由模型解释的残差的比例，对于线性模型来说就是模型拟合的R2（R-squred）。它在 0 和 1 之间，越接近 1 说明模型的表现越好，如果是 0，说明模型的预测结果还不如直接把因变量的均值作为预测值来的有效。列 Lambda 当然就是每个模型对应的 λ 值。我们可以看到，随着 λ 的变小，越来越多的自变量被模型接纳进来，%Dev 也越来越大。第 65 行时，模型包含了所有 8 个自变量，%Dev 也在 0.7 以上。其实我们本应该得到 100个不同的模型，但是连续几个 %Dev 变化很小时 glmnet() 会自动停止。分析模型输出我们可以看到当 Df 大于 6 的时候，%Dev 就达到了 0.66，而且继续缩小 λ，即增加更多的自变量到模型中，也不能显著提高 %Dev。所以我们可以认为当 λλ 接近 0.1 时，得到的包含 6 个自变量的模型，可以相当不错的描述这组数据。

plot(fit)

## 我们也可以通过指定λ值，抓取出某一个模型的系数：predict(fit,s = 0.1,type = 'coefficient')## coef(fit,s = 0.1)

计算一下在测试集的MSE

pred.lasso = predict(fit,newx = as.matrix(test[,-9]),s = 0.1,type = 'response')mean((test$lpsa - pred.lasso)^2)

上期计算了几种方法的MSE，我们来比较一下：

OLS：0.4925

ridge（lambda=0.1）：0.4784

ridge（lambda=optimal-λ）：0.4777

lasso（lambda=0.1）：0.4447

所以目前为止，lasso的泛化能力要比上面的结果都优秀

毕竟这个lambda是通过我们的观察选的，不是很科学，我们接下来用交叉验证来选择一个最优的lambda。

cv_fit <- cv.glmnet(x, y,                     family = 'gaussian',                    alpha = 1,                    # nfolds = 10,                      type.measure="mse")plot(cv_fit)

参数family规定了回归模型的类型:
-family="gaussian"适用于一维连续因变量
-family=mgaussian"适用于多维连续因变量
-family="poisson"适用于非负次数因变量(count)
-family="binomial"适用于二元离散因变量(binary)
-family="multinomial"适用于多元离散因变量(category)

cv_fit$lambda.min  #最佳lambda值cv_fit$lambda.1se #指在lambda.min一个标准差范围内得到的最简单模型的那一个lambda值。因为lambda值达到一定大小之后，继续增加模型自变量个数及缩小lambda值，并不能显著提高模型性能，lambda.lse给出的就是一个具备优良性能但是自变量个数最少的模型coefficients<-coef(fit,s=cv_fit$lambda.min)> 最佳lambda值Active.Index<-which(coefficients!=0)    > 系数不为0的特征索引Active.coefficients<-coefficients[Active.Index]   > 系数不为0的特征系数值

上图中，左边这条线是x=cv_fit$lambda.min，右边这条线是x=cv_fit$lambda.1se

当lambda为最优时，模型此时保留了所有变量

最佳lambda值：cv_fit$lambda.min

输出最佳lambda值时，lasso回归的系数

predict(fit,s = cv_fit$lambda.min,type = 'coefficient')

利用最佳lambda值拟合的lasso模型输出测试集的结果，并计算MSE

pred.lasso2 = predict(fit,newx = as.matrix(test[,-9]),s = cv_fit$lambda.min,type = 'response')mean((test$lpsa - pred.lasso2)^2)

我们发现一个现象，用最优的lambda建模发现泛化误差比以上几种方法的都大，这其实也能理解，毕竟模型复杂度改变不是很多，那么我们选择lambda.lse来建模，lambda.lse给出的就是一个具备优良性能但是自变量个数最少的模型，也就是更简单的模型。

pred.lasso3 = predict(fit,newx = as.matrix(test[,-9]),s = cv_fit$lambda.1se,type = 'response')mean((test$lpsa - pred.lasso3)^2)

这时候发现得到的MSE是这些方法中最小的，这也说明了机器学习的建模不是一种方法就可以走到底的，需要以结果为导向，不断的调试才能得到满意的模型。

参考资料：

周志华-机器学习
知乎-狗熊会-统计学习：变量选择之Lasso
统计之都：Lasso回归
Understanding the Bias-Variance Tradeoff
菜菜的sklearn

【声明】内容源于网络

小明的数据分析笔记本

分享R语言和python在生物信息领域做数据分析和数据可视化的简单小例子；偶尔会分享一些组学数据处理相关的内容

内容 971

粉丝 0

小明的数据分析笔记本分享R语言和python在生物信息领域做数据分析和数据可视化的简单小例子；偶尔会分享一些组学数据处理相关的内容

总阅读218

粉丝0

内容971