拓端数据tecdat|R语言惩罚logistic逻辑回归（LASSO,岭回归）高维变量选择的分类模型案例

原文链接：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)

为什么不尝试一个标准的优化程序呢?我们提到过使用优化例程并不明智，因为它们强烈依赖于起点。

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)))