全文下载链接：http://tecdat.cn/?p=21444
在本文中，逻辑logistic回归是研究中常用的方法,可以进行影响因素筛选、概率预测、分类等，例如医学研究中高通里测序技术得到的数据给高维变量选择问题带来挑战，惩罚logisitc回归可以对高维数据进行变量选择和系数估计，且其有效的算法保证了计算的可行性。方法本文介绍了常用的惩罚logistic算法如LASSO、岭回归。
我们之前已经看到，用于估计参数模型参数的经典估计技术是使用最大似然法。更具体地说， 
这里的目标函数只关注拟合优度。但通常，在计量经济学中，我们相信简单的理论比更复杂的理论更可取。所以我们想惩罚过于复杂的模型。
这主意不错。计量经济学教科书中经常提到这一点，但对于模型的选择，通常不涉及推理。通常，我们使用最大似然法估计参数，然后使用AIC或BIC来比较两个模型。Akaike（AIC）标准是基于
我们在左边有一个拟合优度的度量，而在右边，该罚则随着模型的“复杂性”而增加。
这里，复杂性是使用的变量的数量。但是假设我们不做变量选择，我们考虑所有协变量的回归。定义
AIC是可以写为
实际上，这就是我们的目标函数。更具体地说，我们将考虑
在这篇文章中，我想讨论解决这种优化问题的数值算法，对于l1（岭回归）和l2（LASSO回归）。
协变量的标准化
这里我们使用从急诊室的病人那里观察到的梗塞数据，我们想知道谁活了下来，得到一个预测模型。第一步是考虑所有协变量x_jxj的线性变换来标准化变量(带有单位方差)
for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])
岭回归
在运行一些代码之前，回想一下我们想要解决如下问题
在考虑高斯变量对数似然的情况下，得到残差的平方和，从而得到显式解。但不是在逻辑回归的情况下。
岭回归的启发式方法如下图所示。在背景中，我们可以可视化logistic回归的（二维）对数似然，如果我们将优化问题作为约束优化问题重新布线，蓝色圆圈就是我们的约束： 
可以等效地写（这是一个严格的凸问题）
因此，受约束的最大值应该在蓝色的圆盘上
b0=bbeta[1] beta=bbeta[-1] sum(-y*log(1 + exp(-(b0+X%*%beta))) - (1-y)*log(1 + exp(b0+X%*%beta)))} u = seq(-4,4,length=251) v = outer(u,u,function(x,y) LogLik(c(1,x,y))) lines(u,sqrt(1-u^2),type="l",lwd=2,col="blue") lines(u,-sqrt(1-u^2),type="l",lwd=2,col="blue")
让我们考虑一下目标函数，下面的代码
-sum(-y*log(1 + exp(-(b0+X%*%beta))) - (1-y)* log(1 + exp(b0+X%*%beta)))+lambda*sum(beta^2)
点击标题查阅往期内容
R语言自适应LASSO 多项式回归、二元逻辑回归和岭回归应用分析
01
02
03
04
为什么不尝试一个标准的优化程序呢?我们提到过使用优化例程并不明智，因为它们强烈依赖于起点。
beta_init = lm(y~.,)$coefficients for(i in 1:1000){ vpar[i,] = optim(par = beta_init*rnorm(8,1,2), function(x) LogLik(x,lambda), method = "BFGS", control = list(abstol=1e-9))$par} par(mfrow=c(1,2)) plot(density(vpar[,2])
显然，即使我们更改起点，也似乎我们朝着相同的值收敛。可以认为这是最佳的。
然后将用于计算βλ的代码
beta_init = lm(y~.,data )$coefficients logistic_opt = optim(par = beta_init*0, function(x) LogLik(x,lambda), method = "BFGS", control=list(abstol=1e-9))
我们可以将βλ的演化可视化为λ的函数
v_lambda = c(exp(seq(-2,5,length=61))) plot(v_lambda,est_ridge[1,],col=colrs[1]) for(i in 2:7) lines(v_lambda,est_ridge[i,],
这看起来是有意义的:我们可以观察到λ增加时的收缩。
Ridge，使用Netwon Raphson算法
我们已经看到，我们也可以使用Newton Raphson解决此问题。没有惩罚项，算法是
其中
因此
然后是代码
for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j]) for(s in 1:9){ pi = exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) B = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% (t(X)%*%Delta%*%z) beta = cbind(beta,B)} beta[,8:10] [,1] [,2] [,3] XInter 0.59619654 0.59619654 0.59619654 XFRCAR 0.09217848 0.09217848 0.09217848 XINCAR 0.77165707 0.77165707 0.77165707 XINSYS 0.69678521 0.69678521 0.69678521 XPRDIA -0.29575642 -0.29575642 -0.29575642 XPAPUL -0.23921101 -0.23921101 -0.23921101 XPVENT -0.33120792 -0.33120792 -0.33120792 XREPUL -0.84308972 -0.84308972 -0.84308972
同样，似乎收敛的速度非常快。
有趣的是，通过这个算法，我们还可以得到估计量的方差
然后根据 λ函数计算 βλ的代码
for(s in 1:20){ pi = exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) diag(Delta)=(pi*(1-pi)) z = X%*%beta[,s] + solve(Delta)%*%(Y-pi) B = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% (t(X)%*%Delta%*%z) beta = cbind(beta,B)} Varz = solve(Delta) Varb = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% t(X)%*% Delta %*% Varz %*% Delta %*% X %*% solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X)))
我们可以可视化 βλ的演化（作为 λ的函数）
plot(v_lambda,est_ridge[1,],col=colrs[1],type="l") for(i in 2:7) lines(v_lambda,est_ridge[i,],col=colrs[i])
并获得方差的演变
回想一下，当λ=0（在图的左边），β0=βmco（没有惩罚）。因此，当λ增加时（i）偏差增加（估计趋于0）（ii）方差减小。
使用glmnet Ridge回归
与往常一样，有R个函数可用于进行岭回归。让我们使用glmnet函数， α= 0
for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j]) glmnet(X, y, alpha=0) plot(glm_ridge,xvar="lambda")
作为L1标准范数,
带正交协变量的岭回归当协变量是正交的时，得到了一个有趣的例子。这可以通过协变量的主成分分析得到。
get_pca_ind(pca)$coord
让我们对这些（正交）协变量进行岭回归
glm_ridge = glmnet(pca_X, y, alpha=0) plot(glm_ridge,xvar="lambda",col=colrs,lwd=2)
我们清楚地观察到参数的收缩，即
R语言惩罚logistic逻辑回归（LASSO,岭回归）高维变量选择分类心肌梗塞数据模型案例(下)：https://developer.aliyun.com/article/1493442
