拓端tecdat|R語言懲罰logistic邏輯回歸（LASSO,嶺回歸）高維變量選擇的分類模型案例
原文 |?http://tecdat.cn/?p=21444
來源 |?拓端數(shù)據(jù)部落公眾號
邏輯logistic回歸是研究中常用的方法,可以進行影響因素篩選、概率預(yù)測、分類等，例如醫(yī)學(xué)研究中高通里測序技術(shù)得到的數(shù)據(jù)給高維變量選擇問題帶來挑戰(zhàn)，懲罰logisitc回歸可以對高維數(shù)據(jù)進行變量選擇和系數(shù)估計，且其有效的算法保證了計算的可行性。方法本文介紹了常用的懲罰logistic算法如LASSO、嶺回歸。
方法
我們之前已經(jīng)看到，用于估計參數(shù)模型參數(shù)的經(jīng)典估計技術(shù)是使用最大似然法。更具體地說，
這里的目標函數(shù)只關(guān)注擬合優(yōu)度。但通常，在計量經(jīng)濟學(xué)中，我們相信簡單的理論比更復(fù)雜的理論更可取。所以我們想懲罰過于復(fù)雜的模型。
這主意不錯。計量經(jīng)濟學(xué)教科書中經(jīng)常提到這一點，但對于模型的選擇，通常不涉及推理。通常，我們使用最大似然法估計參數(shù)，然后使用AIC或BIC來比較兩個模型。Akaike（AIC）標準是基于
我們在左邊有一個擬合優(yōu)度的度量，而在右邊，該罰則隨著模型的“復(fù)雜性”而增加。
這里，復(fù)雜性是使用的變量的數(shù)量。但是假設(shè)我們不做變量選擇，我們考慮所有協(xié)變量的回歸。定義
AIC是可以寫為
實際上，這就是我們的目標函數(shù)。更具體地說，我們將考慮
在這篇文章中，我想討論解決這種優(yōu)化問題的數(shù)值算法，對于l1（嶺回歸）和l2（LASSO回歸）。
協(xié)變量的標準化
這里我們使用從急診室的病人那里觀察到的梗塞數(shù)據(jù)，我們想知道誰活了下來，得到一個預(yù)測模型。第一步是考慮所有協(xié)變量x_jxj的線性變換來標準化變量(帶有單位方差)
for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])
嶺回歸
在運行一些代碼之前，回想一下我們想要解決如下問題
在考慮高斯變量對數(shù)似然的情況下，得到殘差的平方和，從而得到顯式解。但不是在邏輯回歸的情況下。
嶺回歸的啟發(fā)式方法如下圖所示。在背景中，我們可以可視化logistic回歸的（二維）對數(shù)似然，如果我們將優(yōu)化問題作為約束優(yōu)化問題重新布線，藍色圓圈就是我們的約束：
可以等效地寫（這是一個嚴格的凸問題）
因此，受約束的最大值應(yīng)該在藍色的圓盤上
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")
?
讓我們考慮一下目標函數(shù)，下面的代碼
-sum(-y*log(1 + exp(-(b0+X%*%beta))) - (1-y)*log(1 + exp(b0+X%*%beta)))+lambda*sum(beta^2)
為什么不嘗試一個標準的優(yōu)化程序呢?我們提到過使用優(yōu)化例程并不明智，因為它們強烈依賴于起點。
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])
?
顯然，即使我們更改起點，也似乎我們朝著相同的值收斂?？梢哉J為這是最佳的。
然后將用于計算βλ的代碼
beta_init = lm(y~.,data )$coefficientslogistic_opt = optim(par = beta_init*0, function(x) LogLik(x,lambda),method = "BFGS", control=list(abstol=1e-9))
我們可以將βλ的演化可視化為λ的函數(shù)
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算法
我們已經(jīng)看到，我們也可以使用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
同樣，似乎收斂的速度非?？?。
有趣的是，通過這個算法，我們還可以得到估計量的方差
然后根據(jù) λ函數(shù)計算 βλ的代碼
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)))
我們可以可視化 βλ的演化（作為 λ的函數(shù)）
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個函數(shù)可用于進行嶺回歸。讓我們使用glmnet函數(shù)， α= 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標準范數(shù),
帶正交協(xié)變量的嶺回歸當協(xié)變量是正交的時，得到了一個有趣的例子。這可以通過協(xié)變量的主成分分析得到。
get_pca_ind(pca)$coord
讓我們對這些（正交）協(xié)變量進行嶺回歸
glm_ridge = glmnet(pca_X, y, alpha=0)plot(glm_ridge,xvar="lambda",col=colrs,lwd=2)
我們清楚地觀察到參數(shù)的收縮，即
應(yīng)用
讓我們嘗試第二組數(shù)據(jù)
我們可以嘗試各種λ的值
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標準范數(shù)。
協(xié)變量的標準化
如前所述，第一步是考慮所有協(xié)變量x_jxj的線性變換，使變量標準化(帶有單位方差)
for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])X = as.matrix(X)
嶺回歸
關(guān)于lasso套索回歸的啟發(fā)式方法如下圖所示。在背景中，我們可以可視化logistic回歸的（二維）對數(shù)似然，藍色正方形是我們的約束條件，如果我們將優(yōu)化問題作為一個約束優(yōu)化問題重新考慮，
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")
這里的好處是它可以用作變量選擇工具。
啟發(fā)性地，數(shù)學(xué)解釋如下。考慮一個簡單的回歸方程y_i=xiβ+ε，具有? l1-懲罰和 l2-損失函數(shù)。優(yōu)化問題變成
一階條件可以寫成
則解為
優(yōu)化過程
讓我們從標準（R）優(yōu)化例程開始，比如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)
?
結(jié)果是不穩(wěn)定的。
使用glmnet
為了進行比較，使用專用于lasso的R程序，我們得到以下內(nèi)容
plot(glm_lasso,xvar="lambda",col=colrs,lwd=2)
如果我們仔細觀察輸出中的內(nèi)容，就可以看到存在變量選擇，就某種意義而言，βj，λ= 0，意味著“真的為零”。
,lambda=exp(-4))$beta7x1 sparse Matrix of class "dgCMatrix"s0FRCAR .INCAR 0.11005070INSYS 0.03231929PRDIA .PAPUL .PVENT -0.03138089REPUL -0.20962611
沒有優(yōu)化例程，我們不能期望有空值
opt_lasso(.2)FRCAR INCAR INSYS PRDIA0.4810999782 0.0002813658 1.9117847987 -0.3873926427PAPUL PVENT REPUL-0.0863050787 -0.4144139379 -1.3849264055
正交協(xié)變量
在學(xué)習數(shù)學(xué)之前，請注意，當協(xié)變量是正交的時，有一些非常清楚的“變量”選擇過程，
pca = princomp(X)pca_X = get_pca_ind(pca)$coordplot(glm_lasso,xvar="lambda",col=colrs)plot(glm_lasso,col=colrs)
標準lasso
如果我們回到原來的lasso方法，目標是解決
注意，截距不受懲罰。一階條件是
也就是
我們可以檢查bf0是否至少包含次微分。
對于左邊的項
這樣前面的方程就可以寫出來了
然后我們將它們簡化為一組規(guī)則來檢查我們的解
我們可以將βj分解為正負部分之和，方法是將βj替換為βj+-βj-，其中βj+，βj-≥0。lasso問題就變成了
令αj+，αj?分別表示βj+，βj?的拉格朗日乘數(shù)。
為了滿足平穩(wěn)性條件，我們?nèi)±窭嗜贞P(guān)于βj+的梯度，并將其設(shè)置為零獲得
我們對βj?進行相同操作以獲得
為了方便起見，引入了軟閾值函數(shù)
注意優(yōu)化問題
也可以寫
觀察到
這是一個坐標更新。
現(xiàn)在，如果我們考慮一個（稍微）更一般的問題，在第一部分中有權(quán)重
坐標更新變?yōu)?br>
回到我們最初的問題。
lasso套索邏輯回歸
這里可以將邏輯問題表述為二次規(guī)劃問題?；叵胍幌聦?shù)似然在這里
這是參數(shù)的凹函數(shù)。因此，可以使用對數(shù)似然的二次近似-使用泰勒展開，
其中z_izi是
pi是預(yù)測
這樣，我們得到了一個懲罰的最小二乘問題。我們可以用之前的方法
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') < tol) { break }}return(list(obj=obj[1:j],beta=beta,intercept=beta0)) }
它看起來像是調(diào)用glmnet時得到的結(jié)果，對于一些足夠大的λ，我們確實有空成分。
在第二個數(shù)據(jù)集上的應(yīng)用
現(xiàn)在考慮具有兩個協(xié)變量的第二個數(shù)據(jù)集。獲取lasso估計的代碼是
plot_l = function(lambda){m = apply(df0,2,mean)s = apply(df0,2,sd)for(j in 1:2) df0[,j] ® = 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的一些小值，我們就只有某種程度的參數(shù)縮小
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邏輯回歸 應(yīng)用案例
2.面板平滑轉(zhuǎn)移回歸(PSTR)分析案例實現(xiàn)
3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）
4.R語言泊松Poisson回歸模型分析案例
5.R語言回歸中的Hosmer-Lemeshow擬合優(yōu)度檢驗
6.r語言中對LASSO回歸，Ridge嶺回歸和Elastic Net模型實現(xiàn)
7.在R語言中實現(xiàn)Logistic邏輯回歸
8.python用線性回歸預(yù)測股票價格
9.R語言如何在生存分析與Cox回歸中計算IDI，NRI指標
