共计 6433 个字符，预计需要花费 17 分钟才能阅读完成。
原文链接：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))) 咱们能够可视化 βλ 的演变（作为 λ 的函数）
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))$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
最受欢迎的见解
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 指标