拓端數據tecdat|R語言懲罰logistic邏輯回歸（LASSO,嶺回歸）高維變量選擇的分類模型案例

原文鏈接：http://tecdat.cn/?p=21444

 

邏輯logistic回歸是研究中常用的方法,可以進行影響因素篩選、概率預測、分類等，例如醫學研究中高通里測序技術得到的數據給高維變量選擇問題帶來挑戰，懲罰logisitc回歸可以對高維數據進行變量選擇和系數估計，且其有效的算法保證了計算的可行性。方法本文介紹了常用的懲罰logistic算法如LASSO、嶺回歸。

方法

我們之前已經看到，用於估計參數模型參數的經典估計技術是使用最大似然法。更具體地說，

這里的目標函數只關注擬合優度。但通常，在計量經濟學中，我們相信簡單的理論比更復雜的理論更可取。所以我們想懲罰過於復雜的模型。

這主意不錯。計量經濟學教科書中經常提到這一點，但對於模型的選擇，通常不涉及推理。通常，我們使用最大似然法估計參數，然后使用AIC或BIC來比較兩個模型。Akaike（AIC）標准是基於

我們在左邊有一個擬合優度的度量，而在右邊，該罰則隨着模型的“復雜性”而增加。

這里，復雜性是使用的變量的數量。但是假設我們不做變量選擇，我們考慮所有協變量的回歸。定義

AIC是可以寫為

實際上，這就是我們的目標函數。更具體地說，我們將考慮

在這篇文章中，我想討論解決這種優化問題的數值算法，對於l1（嶺回歸）和l2（LASSO回歸）。

協變量的標准化

這里我們使用從急診室的病人那里觀察到的梗塞數據，我們想知道誰活了下來，得到一個預測模型。第一步是考慮所有協變量x_jxj的線性變換來標准化變量(帶有單位方差)

1.  
    
2.  
   for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])
3.  
    

 

嶺回歸

在運行一些代碼之前，回想一下我們想要解決如下問題

在考慮高斯變量對數似然的情況下，得到殘差的平方和，從而得到顯式解。但不是在邏輯回歸的情況下。
嶺回歸的啟發式方法如下圖所示。在背景中，我們可以可視化logistic回歸的（二維）對數似然，如果我們將優化問題作為約束優化問題重新布線，藍色圓圈就是我們的約束：

 

可以等效地寫（這是一個嚴格的凸問題）

因此，受約束的最大值應該在藍色的圓盤上

1.  
   b0=bbeta[1]
2.  
   beta=bbeta[-1]
3.  
   sum(-y*log(1 + exp(-(b0+X%*%beta))) -
4.  
   (1-y)*log(1 + exp(b0+X%*%beta)))}
5.  
   u = seq(-4,4,length=251)
6.  
   v = outer(u,u,function(x,y) LogLik(c(1,x,y)))
7.  
    
8.  
    
9.  
   lines(u,sqrt(1-u^2),type="l",lwd=2,col="blue")
10.  
    lines(u,-sqrt(1-u^2),type="l",lwd=2,col="blue")

 

 

讓我們考慮一下目標函數，下面的代碼

1.  
    
2.  
   -sum(-y*log(1 + exp(-(b0+X%*%beta))) - (1-y)*
3.  
   log(1 + exp(b0+X%*%beta)))+lambda*sum(beta^2)

為什么不嘗試一個標准的優化程序呢?我們提到過使用優化例程並不明智，因為它們強烈依賴於起點。

1.  
    
2.  
   beta_init = lm(y~.,)$coefficients
3.  
    
4.  
   for(i in 1:1000){
5.  
   vpar[i,] = optim(par = beta_init*rnorm(8,1,2),
6.  
   function(x) LogLik(x,lambda), method = "BFGS", control = list(abstol=1e-9))$par}
7.  
   par(mfrow=c(1,2))
8.  
   plot(density(vpar[,2])

 

顯然，即使我們更改起點，也似乎我們朝着相同的值收斂。可以認為這是最佳的。

然后將用於計算βλ的代碼

1.  
    
2.  
   beta_init = lm(y~.,data )$coefficients
3.  
   logistic_opt = optim(par = beta_init*0, function(x) LogLik(x,lambda),
4.  
   method = "BFGS", control=list(abstol=1e-9))

我們可以將βλ的演化可視化為λ的函數

1.  
    
2.  
   v_lambda = c(exp(seq(-2,5,length=61)))
3.  
    
4.  
   plot(v_lambda,est_ridge[1,],col=colrs[1])
5.  
   for(i in 2:7) lines(v_lambda,est_ridge[i,],

這看起來是有意義的:我們可以觀察到λ增加時的收縮。

Ridge，使用Netwon Raphson算法

我們已經看到，我們也可以使用Newton Raphson解決此問題。沒有懲罰項，算法是

 

其中

 

因此

 

然后是代碼

1.  
    
2.  
   for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])
3.  
    
4.  
   for(s in 1:9){
5.  
   pi = exp(X%*%beta[,s])/(1+exp(X%*%beta[,s]))
6.  
    
7.  
   B = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% (t(X)%*%Delta%*%z)
8.  
   beta = cbind(beta,B)}
9.  
   beta[,8:10]
10.  
    [,1] [,2] [,3]
11.  
    XInter 0.59619654 0.59619654 0.59619654
12.  
    XFRCAR 0.09217848 0.09217848 0.09217848
13.  
    XINCAR 0.77165707 0.77165707 0.77165707
14.  
    XINSYS 0.69678521 0.69678521 0.69678521
15.  
    XPRDIA -0.29575642 -0.29575642 -0.29575642
16.  
    XPAPUL -0.23921101 -0.23921101 -0.23921101
17.  
    XPVENT -0.33120792 -0.33120792 -0.33120792
18.  
    XREPUL -0.84308972 -0.84308972 -0.84308972

同樣，似乎收斂的速度非常快。

有趣的是，通過這個算法，我們還可以得到估計量的方差

然后根據 λ函數計算 βλ的代碼

1.  
    
2.  
   for(s in 1:20){
3.  
   pi = exp(X%*%beta[,s])/(1+exp(X%*%beta[,s]))
4.  
   diag(Delta)=(pi*(1-pi))
5.  
   z = X%*%beta[,s] + solve(Delta)%*%(Y-pi)
6.  
   B = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% (t(X)%*%Delta%*%z)
7.  
   beta = cbind(beta,B)}
8.  
   Varz = solve(Delta)
9.  
   Varb = solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X))) %*% t(X)%*% Delta %*% Varz %*%
10.  
    Delta %*% X %*% solve(t(X)%*%Delta%*%X+2*lambda*diag(ncol(X)))

我們可以可視化 βλ的演化（作為 λ的函數）

1.  
    
2.  
   plot(v_lambda,est_ridge[1,],col=colrs[1],type="l")
3.  
   for(i in 2:7) lines(v_lambda,est_ridge[i,],col=colrs[i])

 

並獲得方差的演變

 

回想一下，當\λ=0（在圖的左邊），β0=βmco（沒有懲罰）。因此，當λ增加時（i）偏差增加（估計趨於0）（ii）方差減小。

使用glmnetRidge回歸

與往常一樣，有R個函數可用於進行嶺回歸。讓我們使用glmnet函數， α= 0

1.  
    
2.  
   for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])
3.  
   glmnet(X, y, alpha=0)
4.  
   plot(glm_ridge,xvar="lambda")

 

作為L1標准范數,

帶正交協變量的嶺回歸當協變量是正交的時，得到了一個有趣的例子。這可以通過協變量的主成分分析得到。

1.  
    
2.  
   get_pca_ind(pca)$coord

讓我們對這些（正交）協變量進行嶺回歸

1.  
    
2.  
   glm_ridge = glmnet(pca_X, y, alpha=0)
3.  
   plot(glm_ridge,xvar="lambda",col=colrs,lwd=2)

 

我們清楚地觀察到參數的收縮，即

應用

讓我們嘗試第二組數據

我們可以嘗試各種λ的值

1.  
   glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=0)
2.  
   plot(reg,xvar="lambda",col=c("blue","red"),lwd=2)
3.  
   abline(v=log(.2))

或者

1.  
    
2.  
   abline(v=log(1.2))
3.  
   plot_lambda(1.2)

下一步是改變懲罰的標准，使用l1標准范數。

協變量的標准化

如前所述，第一步是考慮所有協變量x_jxj的線性變換，使變量標准化(帶有單位方差)

1.  
    
2.  
   for(j in 1:7) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])
3.  
   X = as.matrix(X)

嶺回歸

關於lasso套索回歸的啟發式方法如下圖所示。在背景中，我們可以可視化logistic回歸的（二維）對數似然，藍色正方形是我們的約束條件，如果我們將優化問題作為一個約束優化問題重新考慮，

1.  
   LogLik = function(bbeta){
2.  
    
3.  
   sum(-y*log(1 + exp(-(b0+X%*%beta))) -
4.  
   (1-y)*log(1 + exp(b0+X%*%beta)))}
5.  
    
6.  
   contour(u,u,v,add=TRUE)
7.  
   polygon(c(-1,0,1,0),c(0,1,0,-1),border="blue")

這里的好處是它可以用作變量選擇工具。

啟發性地，數學解釋如下。考慮一個簡單的回歸方程y_i=xiβ+ε，具有  l1-懲罰和 l2-損失函數。優化問題變成

一階條件可以寫成

則解為

 

優化過程

讓我們從標准（R）優化例程開始，比如BFGS

1.  
    
2.  
   logistic_opt = optim(par = beta_init*0, function(x) PennegLogLik(x,lambda),
3.  
   hessian=TRUE, method = "BFGS", control=list(abstol=1e-9))
4.  
    
5.  
    
6.  
   plot(v_lambda,est_lasso[1,],col=colrs[1],type="l")
7.  
   for(i in 2:7) lines(v_lambda,est_lasso[i,],col=colrs[i],lwd=2)

 

結果是不穩定的。

使用glmnet

為了進行比較，使用專用於lasso的R程序，我們得到以下內容

1.  
    
2.  
   plot(glm_lasso,xvar="lambda",col=colrs,lwd=2)

如果我們仔細觀察輸出中的內容，就可以看到存在變量選擇，就某種意義而言，βj，λ= 0，意味着“真的為零”。

1.  
   ,lambda=exp(-4))$beta
2.  
   7x1 sparse Matrix of class "dgCMatrix"
3.  
   s0
4.  
   FRCAR .
5.  
   INCAR 0.11005070
6.  
   INSYS 0.03231929
7.  
   PRDIA .
8.  
   PAPUL .
9.  
   PVENT -0.03138089
10.  
    REPUL -0.20962611

沒有優化例程，我們不能期望有空值

1.  
    
2.  
   opt_lasso(.2)
3.  
   FRCAR INCAR INSYS PRDIA
4.  
   0.4810999782 0.0002813658 1.9117847987 -0.3873926427
5.  
   PAPUL PVENT REPUL
6.  
   -0.0863050787 -0.4144139379 -1.3849264055

 

正交協變量

在學習數學之前，請注意，當協變量是正交的時，有一些非常清楚的“變量”選擇過程，

1.  
    
2.  
   pca = princomp(X)
3.  
   pca_X = get_pca_ind(pca)$coord
4.  
    
5.  
   plot(glm_lasso,xvar="lambda",col=colrs)
6.  
   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是預測

這樣，我們得到了一個懲罰的最小二乘問題。我們可以用之前的方法

1.  
    
2.  
   beta0 = sum(y-X%*%beta)/(length(y))
3.  
   beta0list[j+1] = beta0
4.  
   betalist[[j+1]] = beta
5.  
   obj[j] = (1/2)*(1/length(y))*norm(omega*(z - X%*%beta -
6.  
   beta0*rep(1,length(y))),'F')^2 + lambda*sum(abs(beta))
7.  
    
8.  
   if (norm(rbind(beta0list[j],betalist[[j]]) -
9.  
   rbind(beta0,beta),'F') < tol) { break }
10.  
    }
11.  
    return(list(obj=obj[1:j],beta=beta,intercept=beta0)) }

它看起來像是調用glmnet時得到的結果，對於一些足夠大的λ，我們確實有空成分。

在第二個數據集上的應用

現在考慮具有兩個協變量的第二個數據集。獲取lasso估計的代碼是

1.  
    
2.  
   plot_l = function(lambda){
3.  
   m = apply(df0,2,mean)
4.  
   s = apply(df0,2,sd)
5.  
   for(j in 1:2) df0[,j] &
6.  
   reg = glmnet(cbind(df0$x1,df0$x2), df0$y==1, alpha=1,lambda=lambda)
7.  
   u = seq(0,1,length=101)
8.  
   p = function(x,y){
9.  
    
10.  
    predict(reg,newx=cbind(x1=xt,x2=yt),type="response")}
11.  
     
12.  
    image(u,u,v,col=clr10,breaks=(0:10)/10)
13.  
    points(df$x1,df$x2,pch=c(1,19)[1+z],cex=1.5)
14.  
    contour(u,u,v,levels = .5,add=TRUE)

考慮 lambda的一些小值，我們就只有某種程度的參數縮小

1.  
    
2.  
   plot(reg,xvar="lambda",col=c("blue","red"),lwd=2)
3.  
   abline(v=exp(-2.8))
4.  
   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指標