基于R语言实现Lasso回归分析

Python032

基于R语言实现Lasso回归分析,第1张

基于R语言实现Lasso回归分析

主要步骤:

将数据存成csv格式,逗号分隔

在R中,读取数据,然后将数据转成矩阵形式

加载lars包,先安装

调用lars函数

确定Cp值最小的步数

确定筛选出的变量,并计算回归系数

具体代码如下:

需要注意的地方:

1、数据读取的方法,这里用的file.choose( ),这样做的好处是,会弹出窗口让你选择你要加载进来的文件,免去了输入路径的苦恼。

2、数据要转为矩阵形式

3、(la) 可以看到R方,这里为0.66,略低

4、图如何看? summary的结果里,第1步是Cp最小的,在图里,看到第1步与横轴0.0的交界处,只有变量1是非0的。所以筛选出的是nongyangungun

Ps: R语言只学习了数据输入,及一些简单的处理,图形可视化部分尚未学习,等论文写完了,再把这部分认真学习一下~~在这里立个flag

原文链接:http://tecdat.cn/?p=20882

1导言

这篇文章探讨了为什么使用广义相加模型 是一个不错的选择。为此,我们首先需要看一下线性回归,看看为什么在某些情况下它可能不是最佳选择。

2回归模型

假设我们有一些带有两个属性Y和X的数据。如果它们是线性相关的,则它们可能看起来像这样:

a<-ggplot(my_data, aes(x=X,y=Y))+geom_point()+

为了检查这种关系,我们可以使用回归模型。线性回归是一种使用X来预测变量Y的方法。将其应用于我们的数据将预测成红线的一组值:

a+geom_smooth(col="red", method="lm")+

这就是“直线方程式”。根据此等式,我们可以从直线在y轴上开始的位置(“截距”或α)开始描述,并且每个单位的x都增加了多少y(“斜率”),我们将它称为x的系数,或称为β)。还有一点自然的波动,如果没有的话,所有的点都将是完美的。我们将此称为“残差”(ϵ)。数学上是:

或者,如果我们用实际数字代替,则会得到以下结果:

这篇文章通过考虑每个数据点和线之间的差异(“残差)然后最小化这种差异来估算模型。我们在线的上方和下方都有正误差和负误差,因此,通过对它们进行平方并最小化“平方和”,使它们对于估计都为正。这称为“普通最小二乘法”或OLS。

3非线性关系如何?

因此,如果我们的数据看起来像这样,我们该怎么办:

我们刚刚看到的模型的关键假设之一是y和x线性相关。如果我们的y不是正态分布的,则使用广义线性模型 (Nelder&Wedderburn,1972),其中y通过链接函数进行变换,但再次假设f(y)和x线性相关。如果不是这种情况,并且关系在x的范围内变化,则可能不是最合适的。我们在这里有一些选择:

我们可以使用线性拟合,但是如果这样做的话,我们会在数据的某些部分上面或者下面。

我们可以分为几类。我在下面的图中使用了三个,这是一个合理的选择。同样,我们可能处于数据某些部分之下或之上,而在类别之间的边界附近似乎是准确的。例如,如果x = 49时,与x = 50相比,y是否有很大不同?

我们可以使用多项式之类的变换。下面,我使用三次多项式,因此模型适合:。这些的组合使函数可以光滑地近似变化。这是一个很好的选择,但可能会极端波动,并可能在数据中引起相关性,从而降低拟合度。

请点击输入图片描述

请点击输入图片描述

4样条曲线

多项式的进一步细化是拟合“分段”多项式,我们在数据范围内将多项式链在一起以描述形状。“样条线”是分段多项式,以绘图员用来绘制曲线的工具命名。物理样条曲线是一种柔性条,可以弯曲成形,并由砝码固定。在构造数学样条曲线时,我们有多项式函数,二阶导数连续,固定在“结”点上。

下面是一个ggplot2 对象,该 对象的 geom_smooth 的公式包含ns 函数中的“自然三次样条”  。这种样条曲线为“三次”,并且使用10个结

请点击输入图片描述

请点击输入图片描述

5光滑函数

样条曲线可以是光滑的或“摇摆的”,这可以通过改变节点数(k)或使用光滑惩罚γ来控制。如果我们增加结的数目,它将更“摇摆”。这可能会更接近数据,而且误差也会更小,但我们开始“过度拟合”关系,并拟合我们数据中的噪声。当我们结合光滑惩罚时,我们会惩罚模型中的复杂度,这有助于减少过度拟合。

请点击输入图片描述

6广义相加模型(GAM)

广义加性模型(GAM)(Hastie,1984)使用光滑函数(如样条曲线)作为回归模型中的预测因子。这些模型是严格可加的,这意味着我们不能像正常回归那样使用交互项,但是我们可以通过重新参数化作为一个更光滑的模型来实现同样的效果。事实并非如此,但本质上,我们正转向一种模型,如:

请点击输入图片描述

摘自Wood (2017)的GAM的更正式示例 是:

请点击输入图片描述

其中:

μi≡E(Yi),Y的期望

Yi〜EF(μi,ϕi),Yi是一个响应变量,根据均值μi和形状参数ϕ的指数族分布。

Ai是任何严格参数化模型分量的模型矩阵的一行,其中θ为对应的参数向量。

fi是协变量xk的光滑函数,其中k是每个函数的基础。

如果您要建立回归模型,但怀疑光滑拟合会做得更好,那么GAM是一个不错的选择。它们适合于非线性或有噪声的数据。

7 gam拟合

那么,如何 为上述S型数据建立 GAM模型?在这里,我将使用三次样条回归 :

gam(Y ~ s(X, bs="cr")

上面的设置意味着:

s()指定光滑器。还有其他选项,但是s是一个很好的默认选项

bs=“cr”告诉它使用三次回归样条('basis')。

s函数计算出要使用的默认结数,但是您可以将其更改为k=10,例如10个结。

8模型输出:

查看模型摘要:

#### Family: gaussian## Link function: identity## Parametric coefficients:##             Estimate Std. Error t value Pr(>|t|)## (Intercept)  43.9659     0.8305   52.94   <2e-16 ***## ---## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### Approximate significance of smooth terms:##        edf Ref.df     F p-value## s(X) 6.087  7.143 296.3  <2e-16 ***## ---## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### R-sq.(adj) =  0.876   Deviance explained = 87.9%## GCV = 211.94  Scale est. = 206.93    n = 300

显示了我们截距的模型系数,所有非光滑参数将在此处显示

每个光滑项的总体含义如下。

这是基于“有效自由度”(edf)的,因为我们使用的样条函数可以扩展为许多参数,但我们也在惩罚它们并减少它们的影响。

9检查模型:

该 gam.check() 函数可用于查看残差图,但它也可以测试光滑器以查看是否有足够的结来描述数据。但是如果p值很低,则需要更多的结。

请点击输入图片描述

#### Method: GCV   Optimizer: magic## Smoothing parameter selection converged after 4 iterations.## The RMS GCV score gradient at convergence was 1.107369e-05 .## The Hessian was positive definite.## Model rank =  10 / 10#### Basis dimension (k) checking results. Low p-value (k-index<1) may## indicate that k is too low, especially if edf is close to k'.####        k'  edf k-index p-value## s(X) 9.00 6.09     1.1    0.97

10它比线性模型好吗?

让我们对比具有相同数据的普通线性回归模型:

anova(my_lm, my_gam)

## Analysis of Variance Table#### Model 1: Y ~ X## Model 2: Y ~ s(X, bs = "cr")##   Res.Df   RSS     Df Sum of Sq      F    Pr(>F)## 1 298.00 88154## 2 292.91 60613 5.0873     27540 26.161 <2.2e-16 ***## ---## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

我们的方差分析函数在这里执行了f检验,我们的GAM模型明显优于线性回归。

11小结

所以,我们看了什么是回归模型,我们是如何解释一个变量y和另一个变量x的。其中一个基本假设是线性关系,但情况并非总是这样。当关系在x的范围内变化时,我们可以使用函数来改变这个形状。一个很好的方法是在“结”点处将光滑曲线链接在一起,我们称之为“样条曲线”

我们可以在常规回归中使用这些样条曲线,但是如果我们在GAM的背景中使用它们,我们同时估计了回归模型以及如何使我们的模型更光滑。

上面的示例显示了基于样条的GAM,其拟合度比线性回归模型好得多。

12参考:

NELDER, J. A. &WEDDERBURN, R. W. M. 1972. Generalized Linear Models. Journal of the Royal Statistical Society. Series A (General), 135, 370-384.

HARRELL, F. E., JR. 2001. Regression Modeling Strategies, New York, Springer-Verlag New York.

请点击输入图片描述

最受欢迎的见解

1.R语言多元Logistic逻辑回归 应用案例

2.面板平滑转移回归(PSTR)分析案例实现

3.matlab中的偏最小二乘回归(PLSR)和主成分回归(PCR)

4.R语言泊松Poisson回归模型分析案例

5.R语言回归中的Hosmer-Lemeshow拟合优度检验

6.r语言中对LASSO回归,Ridge岭回归和Elastic Net模型实现

7.在R语言中实现Logistic逻辑回归

8.python用线性回归预测股票价格

9.R语言如何在生存分析与Cox回归中计算IDI,NRI指标

假设有k个自变量的多元线性回归模型:

其中误差项是一个期望值为0且服从正态分布的随机变量:

则利用最小二乘法可得参数的估计值为:

该求解公式唯一的条件是矩阵X是列满秩的,不然会有无穷多解:

当各变量之间存在共线性问题,即各变量之间存在部分线性相关时,例如:

易知此时X近乎是不满秩的(实际情况很难完全共线性),X^TX近乎是奇异的,X的最小奇异值会非常小,那它的影响到底有多大呢?我们先从 矩阵计算 的角度来看。

对于一个方程或者系统而言,当输入有一个非常微小的扰动时,我们希望方程或系统的输出变化也非常微小,如果输出的变化非常大,且不能被控制,那这个系统的预测就无效了,蝴蝶效应讲的就是这个。在矩阵计算中,这叫做 扰动分析

可以看到矩阵的条件数越大,扰动就越大,即x的求解值会变得非常不准确。回到上面讲的线性回归问题,容易证明最小二乘法的解满足下面的正定方程:

此时

当方程有共线性问题时,X的最小特征值非常小,相应的,上述的条件数会非常大。也就是说机器学习中的共线性问题实际上就是矩阵计算中的条件数问题。 从实际应用的角度,一般若K<100,则认为多重共线性的程度很小,若是100<=K<=1000,则认为存在一般程度上的多重共线性,若是K>1000,则就认为存在严重的多重共线性。

再从统计学的角度来看共线性。可以证明参数$\theta$的协方差矩阵为

又对任意的常数矩阵A和随机变量x有

代入上式即可得

具体到每个参数,有:

其中$R i 2 $是将第i个变量$x_i$作为因变量,其他k-1个变量作为自变量进行线性回归获得的$R 2 $,且令

方差膨胀因子 ( variance inflation factor ,VIF)。当

时,即当第i个变量和其他变量之间存在线性关系时,VIF趋于无穷大。所以 VIF 的大小反应了变量的共线性程度。一般地,当VIF大于5或10时,认为模型存在严重的共线性问题。

同时考虑参数显著性检验的 t 统计量

当存在共线性时,参数的标准差偏大,相应的 t 统计量 会偏小,这样容易淘汰一些不应淘汰的解释变量,使统计检验的结果失去可靠性。

另外考虑线性回归的残差

其中M是一个投影矩阵,且满足

易证明

而矩阵M的范数与X的条件数毫无关系,于是可以得出 共线性并不影响模型的训练精度 。但是对于泛化精度,由于参数的估计已经不准确啦,所以泛化误差肯定要差些,具体差多少,我还很难用公式表示出来。

总结一下,共线性问题对线性回归模型有如下影响:

根据上一节的描述,共线性问题有如下几种检验方法:

当变量数不多,样本数不是很大时,上述的方法是没问题的,检验某个变量有共线性问题时,可以结合实际业务考虑直接剔除该变量。但是有的时候变量数大到有上千个,VIF的计算需要建立上千个回归模型(条件数仅能判定是否存在共线性,但不能找到对应的变量),这将耗费很长时间。

事实上我们可以从模型角度来直接规避共线性问题。

主成分分析法作为多元统计分析的一种常用方法在处理多变量问题时具有其一定的优越性,其降维的优势是明显的,主成分回归方法对于一般的多重共线性问题还是适用的,尤其是对共线性较强的变量之间。当采取主成分提取了新的变量后,往往这些变量间的组内差异小而组间差异大,起到了消除共线性的问题。

逐步回归(Stepwise Regression)是一种常用的消除多重共线性、选取“最优”回归方程的方法。其做法是将逐个引入自变量,引入的条件是该自变量经F检验是显著的,每引入一个自变量后,对已选入的变量进行逐个检验,如果原来引入的变量由于后面变量的引入而变得不再显著,那么就将其剔除。引入一个变量或从回归方程中剔除一个变量,为逐步回归的一步,每一步都要进行F 检验,以确保每次引入新变量之前回归方程中只包含显著的变量。这个过程反复进行,直到既没有不显著的自变量选入回归方程,也没有显著自变量从回归方程中剔除为止。

岭回归是一种可用于共线性数据分析的有偏估计回归方法,它是一种改良的最小二乘估计法,通过放弃最小二乘法的无偏性,以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法,对条件数很大(病态数据)的拟合要强于最小二乘法。

在线性回归问题中,最小二乘法实际上是最小化问题:

而岭回归则是加入了L2惩罚项:

这样参数的方差不会过大,且随着惩罚项系数C的增大,共线性的影响将越来也小。在这个过程中,可以记录$\theta(k)$(岭迹)的变化情况,通过对岭迹的波动来判断我们是否要剔除该变量。

那为什么说岭回归能解决共线性问题呢?从矩阵计算的角度来看,L2正则化下方程的解为:

在上一节我们讲到共线性代表正定矩阵X T X的条件数很大:

而当条件数很大时,矩阵的逆的数值计算也是非常不准确的,但是当我们给矩阵加上一个单位矩阵时,奇异性(不可逆)问题就完全没有啦。

进一步考虑对惩罚项对奇异值的影响,假设X的奇异值(SVD)分解为:

则容易证明

其中D是对角矩阵,且满足

其反应了惩罚项是如何影响到条件数的。

LASSO回归和岭回归类似,只不过将惩罚项由L2范数改为了L1范数

L1范数没有L2范数那么圆润,毕竟存在不可导点,而且在L1范数下LASSO回归也给不出解析解啦,但是相对于岭回归,LASSO估计的参数能更容易收敛到0

ElasticNet回归同时兼顾了L1和L2惩罚项:

当许多变量是相关的时候,Elastic-net是有用的。Lasso一般会随机选择其中一个,而Elastic-net则会选在两个。

除此之外,还有L0范数(非零元的个数)、L1/2范数等。

首先捏造一份好的数据,样本量为100,特征数为8,且满足方程:

其中误差项是期望为0,标准差为1.5的正态分布随机变量。

此时平均准确率为0.934955,拟合的系数MSE为0.203657

然后我们基于这份数据另外构造出两份数据,第二份数据增加两个随机的特征用作对比,第一份数据则增加两个共线性特征:

先来看下它们的条件数

可以看到X2的条件数很搭,最小奇异值为0.213,此时还不至于完全共线性。

拿这两份数据重新用线性回归拟合模型。

对于第二份共线性构造数据X2,有平均测试集准确率为0.932070,拟合的参数MSE为7.697837。可以看到MSE增加了很多,准确率也下降了0.2%,测试拟合的系数为:

在来看对比用的数据X3,其平均测试集准确率为0.934952,参数MSE为0.171651,与X1无异。

以上是直接的结果,我们再来看VIF

可以看到第0、1、2、3、8、9个特征的VIF都过高。且可以看出第1个特征相对第0、2、3个特征的VIF较高。

最后我们试着用模型的方法来检测共线性问题

其中当alpha取0.1时,岭回归估计的系数分别为

可以看到第0、1、2、3、8、9个变量都出现了波动,代表它们之间存在一定的共线性。观察岭迹,我们可以考虑剔除其中波动比较大的第1、8、9个变量。

另外Lasso回归类似,可以用sklearn中的linear_model.Lasso来学习,这里就不展示了。最后对于逻辑回归任务,sklearn函数内部提供了L1或L2正则化方案,通过它们也可以去检测共线性问题。

[1]. variance inflation factor

[2]. 多重共线性的解决方法之——岭回归与LASSO

[3].ridge regression