原文链接：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~.,)$coefficientsfor(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 )$coefficientslogistic_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.59619654XFRCAR 0.09217848 0.09217848 0.09217848XINCAR 0.77165707 0.77165707 0.77165707XINSYS 0.69678521 0.69678521 0.69678521XPRDIA -0.29575642 -0.29575642 -0.29575642XPAPUL -0.23921101 -0.23921101 -0.23921101XPVENT -0.33120792 -0.33120792 -0.33120792XREPUL -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）方差减小。

应用glmnetRidge回归

与平常一样，有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)

咱们分明地察看到参数的膨胀，即

利用

让咱们尝试第二组数据

咱们能够尝试各种的值

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))$beta7x1 sparse Matrix of class "dgCMatrix"s0FRCAR . INCAR 0.11005070INSYS 0.03231929PRDIA . PAPUL . PVENT -0.03138089REPUL -0.20962611

没有优化例程，咱们不能冀望有空值

 opt_lasso(.2)FRCAR INCAR INSYS PRDIA0.4810999782 0.0002813658 1.9117847987 -0.3873926427PAPUL PVENT REPUL -0.0863050787 -0.4144139379 -1.3849264055

正交协变量

在学习数学之前，请留神，当协变量是正交的时，有一些十分分明的“变量”抉择过程，

 pca = princomp(X)pca_X = get_pca_ind(pca)$coordplot(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] = beta0betalist[[j+1]] = betaobj[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') &lt; 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

最受欢迎的见解

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指标