Glmnet是一个通过惩罚最大似然关系拟合广义线性模型的软件包。正则化路径是针对正则化参数λ的值网格处的lasso或Elastic Net(弹性网络)惩罚值计算的。该算法非常快,并且可以利用输入矩阵中的稀疏性 x
。它适合线性,逻辑和多项式,泊松和Cox回归模型。可以从拟合模型中做出各种预测。它也可以拟合多元线性回归。
glmnet
解决以下问题
在覆盖整个范围的λ值网格上。这里l(y,η)是观察i的负对数似然贡献;例如对于高斯分布是 。 _弹性网络_惩罚由α控制,LASSO(α= 1,默认),Ridge(α= 0)。调整参数λ控制惩罚的总强度。
众所周知,岭惩罚使相关预测因子的系数彼此缩小,而套索倾向于选择其中一个而丢弃其他预测因子。_弹性网络_则将这两者混合在一起。
glmnet
算法使用循环坐标下降法,该方法在每个参数固定不变的情况下连续优化目标函数,并反复循环直到收敛,我们的算法可以非常快速地计算求解路径。
代码可以处理稀疏的输入矩阵格式,以及系数的范围约束,还包括用于预测和绘图的方法,以及执行K折交叉验证的功能。
快速开始
首先,我们加载 glmnet
包:
library(glmnet)
包中使用的默认模型是高斯线性模型或“最小二乘”。我们加载一组预先创建的数据以进行说明。用户可以加载自己的数据,也可以使用工作空间中保存的数据。
该命令 从此保存的R数据中加载输入矩阵 x
和因向量 y
。
我们拟合模型 glmnet
。
fit = glmnet(x, y)
可以通过执行plot
函数来可视化系数 :
plot(fit)
每条曲线对应一个变量。它显示了当λ变化时,其系数相对于整个系数向量的ℓ1范数的路径。上方的轴表示当前λ处非零系数的数量,这是套索的有效自由度(_df_)。用户可能还希望对曲线进行注释。这可以通过label = TRUE
在plot命令中进行设置来完成 。
如果我们只是输入对象名称或使用
glmnetprint
函数,则会显示每个步骤的路径 摘要 :
print(fit)
## ## Call: glmnet(x = x, y = y) ## ## Df %Dev Lambda ## \[1,\] 0 0.0000 1.63000 ## \[2,\] 2 0.0553 1.49000 ## \[3,\] 2 0.1460 1.35000 ## \[4,\] 2 0.2210 1.23000 ## \[5,\] 2 0.2840 1.12000 ## \[6,\] 2 0.3350 1.02000 ## \[7,\] 4 0.3900 0.93300 ## \[8,\] 5 0.4560 0.85000 ## \[9,\] 5 0.5150 0.77500 ## \[10,\] 6 0.5740 0.70600 ## \[11,\] 6 0.6260 0.64300 ## \[12,\] 6 0.6690 0.58600 ## \[13,\] 6 0.7050 0.53400 ## \[14,\] 6 0.7340 0.48700 ## \[15,\] 7 0.7620 0.44300 ## \[16,\] 7 0.7860 0.40400 ## \[17,\] 7 0.8050 0.36800 ## \[18,\] 7 0.8220 0.33500 ## \[19,\] 7 0.8350 0.30600 ## \[20,\] 7 0.8460 0.27800
它从左到右显示了非零系数的数量(Df
),解释的(零)偏差百分比(%dev
)和λ(Lambda
)的值。
我们可以在序列范围内获得一个或多个λ处的实际系数:
coef(fit,s=0.1)
## 21 x 1 sparse Matrix of class "dgCMatrix" ## 1 ## (Intercept) 0.150928 ## V1 1.320597 ## V2 . ## V3 0.675110 ## V4 . ## V5 -0.817412 ## V6 0.521437 ## V7 0.004829 ## V8 0.319416 ## V9 . ## V10 . ## V11 0.142499 ## V12 . ## V13 . ## V14 -1.059979 ## V15 . ## V16 . ## V17 . ## V18 . ## V19 . ## V20 -1.021874
还可以使用新的输入数据在特定的λ处进行预测:
predict(fit,newx=nx,s=c(0.1,0.05))
## 1 2 ## \[1,\] 4.4641 4.7001 ## \[2,\] 1.7509 1.8513 ## \[3,\] 4.5207 4.6512 ## \[4,\] -0.6184 -0.6764 ## \[5,\] 1.7302 1.8451 ## \[6,\] 0.3565 0.3512 ## \[7,\] 0.2881 0.2662 ## \[8,\] 2.7776 2.8209 ## \[9,\] -3.7016 -3.7773 ## \[10,\] 1.1546 1.1067
该函数 glmnet
返回一系列模型供用户选择。交叉验证可能是该任务最简单,使用最广泛的方法。
cv.glmnet
是交叉验证的主要函数。
cv.glmnet
返回一个 cv.glmnet
对象,此处为“ cvfit”,其中包含交叉验证拟合的所有成分的列表。
我们可以绘制对象。
它包括交叉验证曲线(红色虚线)和沿λ序列的上下标准偏差曲线(误差线)。垂直虚线表示两个选定的λ。
我们可以查看所选的λ和相应的系数。例如,
cvfit$lambda.min
## \[1\] 0.08307
lambda.min
是给出最小平均交叉验证误差的λ值。保存的另一个λ是 lambda.1se
,它给出了的模型,使得误差在最小值的一个标准误差以内。我们只需要更换 lambda.min
到lambda.1se
以上。
coef(cvfit, s = "lambda.min")
## 21 x 1 sparse Matrix of class "dgCMatrix" ## 1 ## (Intercept) 0.14936 ## V1 1.32975 ## V2 . ## V3 0.69096 ## V4 . ## V5 -0.83123 ## V6 0.53670 ## V7 0.02005 ## V8 0.33194 ## V9 . ## V10 . ## V11 0.16239 ## V12 . ## V13 . ## V14 -1.07081 ## V15 . ## V16 . ## V17 . ## V18 . ## V19 . ## V20 -1.04341
注意,系数以稀疏矩阵格式表示。原因是沿着正则化路径的解通常是稀疏的,因此使用稀疏格式在时间和空间上更为有效。
可以根据拟合的cv.glmnet
对象进行预测 。让我们看一个示例。
## 1 ## \[1,\] -1.3647 ## \[2,\] 2.5686 ## \[3,\] 0.5706 ## \[4,\] 1.9682 ## \[5,\] 1.4964
newx
与新的输入矩阵 s
相同,如前所述,是预测的λ值。
线性回归
这里的线性回归是指两个模型系列。一个是 gaussian
正态_分布_,另一个是 mgaussian
多元正态_分布_。
正态_分布_
假设我们有观测值xi∈Rp并且yi∈R,i = 1,...,N。目标函数是
其中λ≥0是复杂度参数,0≤α≤1在岭回归(α=0)和套索LASSO(α=1)之间。
应用坐标下降法解决该问题。具体地说,通过计算βj=β〜j处的梯度和简单的演算,更新为
其中 。
当x
变量标准化为具有单位方差(默认值)时,以上公式适用 。
glmnet
提供各种选项供用户自定义。我们在这里介绍一些常用的选项,它们可以在glmnet
函数中指定 。
alpha
表示弹性网混合参数α,范围α∈[0,1]。α=1是套索(默认),α=0是Ridge。weights
用于观察权重。每个观察值的默认值为1。nlambda
是序列中λ值的数量。默认值为100。lambda
可以提供,但通常不提供,程序会构建一个序列。自动生成时,λ序列由lambda.max
和 确定lambda.min.ratio
。standardize
是x
在拟合模型序列之前进行变量标准化的逻辑标志 。
例如,我们设置α=0.2,并对后半部分的观测值赋予两倍的权重。为了避免在此处显示太长时间,我们将其设置 nlambda
为20。但是,实际上,建议将λ的数量设置为100(默认值)或更多。
然后我们可以输出glmnet
对象。
print(fit)
## ## Call: glmnet(x = x, y = y, weights = c(rep(1, 50), rep(2, 50)), alpha = 0.2, nlambda = 20) ## ## Df %Dev Lambda ## \[1,\] 0 0.000 7.94000 ## \[2,\] 4 0.179 4.89000 ## \[3,\] 7 0.444 3.01000 ## \[4,\] 7 0.657 1.85000 ## \[5,\] 8 0.785 1.14000 ## \[6,\] 9 0.854 0.70300 ## \[7,\] 10 0.887 0.43300 ## \[8,\] 11 0.902 0.26700 ## \[9,\] 14 0.910 0.16400 ## \[10,\] 17 0.914 0.10100 ## \[11,\] 17 0.915 0.06230 ## \[12,\] 17 0.916 0.03840 ## \[13,\] 19 0.916 0.02360 ## \[14,\] 20 0.916 0.01460 ## \[15,\] 20 0.916 0.00896 ## \[16,\] 20 0.916 0.00552 ## \[17,\] 20 0.916 0.00340
这将显示生成对象的调用 fit
以及带有列Df
(非零系数的数量), %dev
(解释的偏差百分比)和Lambda
(对应的λ值) 的三列矩阵 。
我们可以绘制拟合的对象。
让我们针对log-lambda值标记每个曲线来绘制“拟合”。
这是训练数据中的偏差百分比。我们在这里看到的是,在路径末端时,该值变化不大,但是系数有点“膨胀”。这使我们可以将注意力集中在重要的拟合部分上。
我们可以提取系数并在某些特定值的情况下进行预测。两种常用的选项是:
s
指定进行提取的λ值。exact
指示是否需要系数的精确值。
一个简单的例子是:
## 21 x 2 sparse Matrix of class "dgCMatrix" ## 1 1 ## (Intercept) 0.19657 0.199099 ## V1 1.17496 1.174650 ## V2 . . ## V3 0.52934 0.531935 ## V4 . . ## V5 -0.76126 -0.760959 ## V6 0.46627 0.468209 ## V7 0.06148 0.061927 ## V8 0.38049 0.380301 ## V9 . . ## V10 . . ## V11 0.14214 0.143261 ## V12 . . ## V13 . . ## V14 -0.91090 -0.911207 ## V15 . . ## V16 . . ## V17 . . ## V18 . 0.009197 ## V19 . . ## V20 -0.86099 -0.863117
左列是,exact = TRUE
右列是 FALSE
。从上面我们可以看到,0.01不在序列中,因此尽管没有太大差异,但还是有一些差异。如果没有特殊要求,则线性插补就足够了。
用户可以根据拟合的对象进行预测。除中的选项外 coef
,主要参数是 newx
的新值矩阵 x
。type
选项允许用户选择预测类型:*“链接”给出拟合值
- 因变量与正态分布的“链接”相同。
- “系数”计算值为的系数
s
例如,
## 1 ## \[1,\] -0.9803 ## \[2,\] 2.2992 ## \[3,\] 0.6011 ## \[4,\] 2.3573 ## \[5,\] 1.7520
给出在λ=0.05时前5个观测值的拟合值。如果提供的多个值, s
则会生成预测矩阵。
用户可以自定义K折交叉验证。除所有 glmnet
参数外, cv.glmnet
还有特殊的参数,包括 nfolds
(次数), foldid
(用户提供的次数), type.measure
(用于交叉验证的损失):*“ deviance”或“ mse”
- “ mae”使用平均绝对误差
举个例子,
cvfit = cv.glmnet(x, y, type.measure = "mse", nfolds = 20)
根据均方误差标准进行20折交叉验证。
并行计算也受 cv.glmnet
。为我们在这里给出一个简单的比较示例。
system.time(cv.glmnet(X, Y))
## user system elapsed ## 3.591 0.103 3.724
system.time(cv.glmnet(X, Y, parallel = TRUE))
## user system elapsed ## 4.318 0.391 2.700
从上面的建议可以看出,并行计算可以大大加快计算过程。
- “ lambda.min”:达到最小MSE的λ。
cvfit$lambda.min
## \[1\] 0.08307
## 21 x 1 sparse Matrix of class "dgCMatrix" ## 1 ## (Intercept) 0.14936 ## V1 1.32975 ## V2 . ## V3 0.69096 ## V4 . ## V5 -0.83123 ## V6 0.53670 ## V7 0.02005 ## V8 0.33194 ## V9 . ## V10 . ## V11 0.16239 ## V12 . ## V13 . ## V14 -1.07081 ## V15 . ## V16 . ## V17 . ## V18 . ## V19 . ## V20 -1.04341
在这里,我们使用相同的k折,为α选择一个值。
将它们全部放置在同一绘图上:
我们看到lasso(alpha=1
)在这里表现最好。
系数上下限
假设我们要拟合我们的模型,但将系数限制为大于-0.7且小于0.5。这可以通过upper.limits
和 lower.limits
参数实现 :
通常,我们希望系数为正,因此我们只能lower.limit
将其设置 为0。
r语言中对LASSO回归,Ridge岭回归和弹性网络Elastic Net模型实现2:https://developer.aliyun.com/article/1485109