R语言惩罚logistic逻辑回归（LASSO,岭回归）高维变量选择分类心肌梗塞数据模型案例(上)：/article/1493440
应用
让我们尝试第二组数据
我们可以尝试各种λ的值
glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=0) plot(reg,xvar="lambda",col=c("blue","red"),lwd=2) abline(v=log(.2))
或者
abline(v=log(1.2)) plot_lambda(1.2)
下一步是改变惩罚的标准，使用l1标准范数。
协变量的标准化
如前所述，第一步是考虑所有协变量x_jxj的线性变换，使变量标准化(带有单位方差)
for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j]) X = as.matrix(X)
岭回归
关于lasso套索回归的启发式方法如下图所示。在背景中，我们可以可视化logistic回归的（二维）对数似然，蓝色正方形是我们的约束条件，如果我们将优化问题作为一个约束优化问题重新考虑，
LogLik = function(bbeta){ sum(-y*log(1 + exp(-(b0+X%*%beta))) - (1-y)*log(1 + exp(b0+X%*%beta)))} contour(u,u,v,add=TRUE) polygon(c(-1,0,1,0),c(0,1,0,-1),border="blue")
这里的好处是它可以用作变量选择工具。
启发性地，数学解释如下。考虑一个简单的回归方程y_i=xiβ+ε，具有 l1-惩罚和 l2-损失函数。优化问题变成
一阶条件可以写成
则解为
优化过程
让我们从标准（R）优化例程开始，比如BFGS
logistic_opt = optim(par = beta_init*0, function(x) PennegLogLik(x,lambda), hessian=TRUE, method = "BFGS", control=list(abstol=1e-9)) plot(v_lambda,est_lasso[1,],col=colrs[1],type="l") for(i in 2:7) lines(v_lambda,est_lasso[i,],col=colrs[i],lwd=2)
结果是不稳定的。
使用glmnet
为了进行比较，使用专用于lasso的R程序，我们得到以下内容
plot(glm_lasso,xvar="lambda",col=colrs,lwd=2)
如果我们仔细观察输出中的内容，就可以看到存在变量选择，就某种意义而言，βj，λ= 0，意味着“真的为零”。
,lambda=exp(-4))$beta 7x1 sparse Matrix of class "dgCMatrix" s0 FRCAR . INCAR 0.11005070 INSYS 0.03231929 PRDIA . PAPUL . PVENT -0.03138089 REPUL -0.20962611
没有优化例程，我们不能期望有空值
opt_lasso(.2) FRCAR INCAR INSYS PRDIA 0.4810999782 0.0002813658 1.9117847987 -0.3873926427 PAPUL PVENT REPUL -0.0863050787 -0.4144139379 -1.3849264055
正交协变量
在学习数学之前，请注意，当协变量是正交的时，有一些非常清楚的“变量”选择过程，
pca = princomp(X) pca_X = get_pca_ind(pca)$coord plot(glm_lasso,xvar="lambda",col=colrs) plot(glm_lasso,col=colrs)
-
标准lasso
如果我们回到原来的lasso方法，目标是解决
注意，截距不受惩罚。一阶条件是
也就是
我们可以检查bf0是否至少包含次微分。
对于左边的项
这样前面的方程就可以写出来了
然后我们将它们简化为一组规则来检查我们的解
我们可以将βj分解为正负部分之和，方法是将βj替换为βj+-βj-，其中βj+，βj-≥0。lasso问题就变成了
令αj+，αj?分别表示βj+，βj?的拉格朗日乘数。
为了满足平稳性条件，我们取拉格朗日关于βj+的梯度，并将其设置为零获得
我们对βj?进行相同操作以获得
为了方便起见，引入了软阈值函数
注意优化问题
也可以写
观察到
这是一个坐标更新。
现在，如果我们考虑一个（稍微）更一般的问题，在第一部分中有权重
坐标更新变为
回到我们最初的问题。
lasso套索逻辑回归
这里可以将逻辑问题表述为二次规划问题。回想一下对数似然在这里
这是参数的凹函数。因此，可以使用对数似然的二次近似-使用泰勒展开，
其中z_izi是
pi是预测
这样，我们得到了一个惩罚的最小二乘问题。我们可以用之前的方法
beta0 = sum(y-X%*%beta)/(length(y)) beta0list[j+1] = beta0 betalist[[j+1]] = beta obj[j] = (1/2)*(1/length(y))*norm(omega*(z - X%*%beta - beta0*rep(1,length(y))),'F')^2 + lambda*sum(abs(beta)) if (norm(rbind(beta0list[j],betalist[[j]]) - rbind(beta0,beta),'F') < tol) { break } } return(list(obj=obj[1:j],beta=beta,intercept=beta0)) }
它看起来像是调用glmnet时得到的结果，对于一些足够大的λ，我们确实有空成分。
在第二个数据集上的应用
现在考虑具有两个协变量的第二个数据集。获取lasso估计的代码是
plot_l = function(lambda){ m = apply(df0,2,mean) s = apply(df0,2,sd) for(j in 1:2) df0[,j] & reg = glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=1,lambda=lambda) u = seq(0,1,length=101) p = function(x,y){ predict(reg,newx=cbind(x1=xt,x2=yt),type="response")} image(u,u,v,col=clr10,breaks=(0:10)/10) points(df$x1,df$x2,pch=c(1,19)[1+z],cex=1.5) contour(u,u,v,levels = .5,add=TRUE)
考虑 lambda的一些小值，我们就只有某种程度的参数缩小
plot(reg,xvar="lambda",col=c("blue","red"),lwd=2) abline(v=exp(-2.8)) plot_l(exp(-2.8))
但是使用较大的λ，则存在变量选择：β1，λ= 0
