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

原文出處：拓端數(shù)據(jù)部落公眾號

本文我們對邏輯回歸和樣條曲線進(jìn)行介紹。

logistic回歸基于以下假設(shè)：給定協(xié)變量x，Y具有伯努利分布，

?

目的是估計(jì)參數(shù)β。

回想一下，針對該概率使用該函數(shù)是

?

（對數(shù)）似然函數(shù)

對數(shù)似然

?

其中

。數(shù)值方法基于（數(shù)值）下降梯度來計(jì)算似然函數(shù)的?最大值。對數(shù)似然（負(fù)）是以下函數(shù)

1. 


2. negLogLik = function(beta){

3. -sum(-y*log(1 + exp(-(X%*%beta))) - (1-y)*log(1 + exp(X%*%beta)))

4. }

現(xiàn)在，我們需要一個起始點(diǎn)來啟動算法

optim(par = beta_init, negLogLik, hessian=TRUE, method = "BFGS", control=list(abstol=1e-9))

在這里，我們得到

1. logistic_opt$par

2. (Intercept) ? ? ? ?FRCAR ? ? ? ?INCAR ? ? ? ?INSYS

3. 1.656926397 ?0.045234029 -2.119441743 ?0.204023835

4. PRDIA ? ? ? ?PAPUL ? ? ? ?PVENT ? ? ? ?REPUL

5. -0.102420095 ?0.165823647 -0.081047525 -0.005992238

讓我們在這里驗(yàn)證該輸出是否有效。例如，如果我們（隨機(jī)）更改起點(diǎn)的值會怎么樣

1. 


2. plot(v_beta)

3. par(mfrow=c(1,2))

4. hist(v_beta[,1],xlab=names( )[ ])

5. hist(v_beta[,2],xlab=names( )[2])

這里有個問題。注意，我們不能在這里進(jìn)行數(shù)值優(yōu)化。我們可以考慮使用其他優(yōu)化方法

1. 


2. logLikelihoodLogitStable = function(vBeta, mX, vY) {

3. -sum(vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta) +

4. (1-vY)*(-log(1 + exp(mX %*% vBeta))

5. 


6. 


7. optimLogitLBFGS = optimx(beta_init, logLikelihoodLogitStable,

最優(yōu)點(diǎn)

結(jié)果不理想。

我們使用的技術(shù)基于以下思想，

?

問題是我的計(jì)算機(jī)不知道一階和二階導(dǎo)數(shù)。

可以使用這種計(jì)算的函數(shù)

1. 


2. logit = function(x){1/(1+exp(-x))}

3. 


4. for(i in 1:num_iter){

5. grad = (t(X)%*%(logit(X%*%beta) - y))

6. beta = beta - ginv(H)%*%grad

7. LL[i] = logLik(beta, X, y)

以我們的OLS起點(diǎn)，我們獲得

如果我們嘗試另一個起點(diǎn)

一些系數(shù)非常接近。然后我們嘗試其他方法。

牛頓（或費(fèi)舍爾）算法

在計(jì)量經(jīng)濟(jì)學(xué)教科書里，您可以看到：

?

?

1. 


2. beta=as.matrix(lm(Y~0+X)$coefficients

3. for(s in 1:9){

4. pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s]))

5. gradient=t(X)%*%(Y-pi)

6. omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi))

在這里觀察到，我僅使用該算法的十次迭代。

事實(shí)是，收斂似乎非?？?。而且它相當(dāng)魯棒，看看我們改變起點(diǎn)會得到什么

1. beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1)*runif(8)

2. for(s in 1:9){

3. pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s]))

4. gradient=t(X)%*%(Y-pi)

5. omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi))

6. Hessian=-t(X)%*%omega%*%X

7. beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)}

8. beta[,8:10]

效果提高了，并且可以使用矩陣的逆獲得標(biāo)準(zhǔn)偏差。

標(biāo)準(zhǔn)最小二乘

我們更進(jìn)一步。我們已經(jīng)看到想要計(jì)算類似

?

?但是實(shí)際，這是一個標(biāo)準(zhǔn)的最小二乘問題

?

這里唯一的問題是權(quán)重Δold是未知β的函數(shù)。但是實(shí)際上，如果我們繼續(xù)迭代，我們應(yīng)該能夠解決它：給定β，我們得到了權(quán)重，并且有了權(quán)重，我們可以使用加權(quán)的OLS來獲取更新的β。這就是迭代最小二乘的想法。

該算法

1. 


2. beta_init = lm(PRONO~.,data=df)$coefficients

3. 


4. for(s in 1:1000){

5. omega = diag(nrow(df))

6. diag(omega) = (p*(1-p))

7. 


輸出在這里

結(jié)果很好，我們在這里也有估計(jì)量的標(biāo)準(zhǔn)差

標(biāo)準(zhǔn)邏輯回歸glm函數(shù)：

當(dāng)然，可以使用R內(nèi)置函數(shù)

可視化

讓我們在第二個數(shù)據(jù)集上可視化從邏輯回歸獲得的預(yù)測

1. 


2. image(u,u,v ,breaks=(0:10)/10)

3. points(x,y,pch=19 )

4. points(x,y,pch=c(1,19)

5. contour(u,u,v,levels = .5

這里的水平曲線-或等概率-是線性的，因此該空間被一條直線（或更高維的超平面）一分為二（0和1，生存和死亡，白色和黑色）此外，由于我們是線性模型，因此，如果更改截距（為創(chuàng)建兩個類別的閾值），我們將獲得平行的另一條直線（或超平面）。

接下來，我們將約會樣條曲線以平滑那些連續(xù)的協(xié)變量。

分段線性樣條函數(shù)

我們從“簡單”回歸開始（只有一個解釋變量），我們可以想到的最簡單的模型來擴(kuò)展我們上面的線性模型，?是考慮一個分段線性函數(shù)，它分為兩部分。最方便的方法是使用正部函數(shù)

（如果該差為正，則為x和s之間的差，否則為0）。如

是以下連續(xù)的分段線性函數(shù)，在s處劃分。

對于較小的x值，線性增加，斜率β1；對于較大的x值，線性減少。因此，β2被解釋為斜率的變化。

當(dāng)然，可以考慮多個結(jié)。獲得正值的函數(shù)如下

pos = function(x,s) (x-s)*(x<=s)

然后我們可以在回歸模型中直接使用它

回歸的輸出在這里

1. 


2. 


3. Coefficients:

4. Estimate Std. Error z value Pr(>|z|)

5. (Intercept) ? ? -0.1109 ? ? 3.2783 ?-0.034 ? 0.9730

6. INSYS ? ? ? ? ? -0.1751 ? ? 0.2526 ?-0.693 ? 0.4883

7. pos(INSYS, 15) ? 0.7900 ? ? 0.3745 ? 2.109 ? 0.0349 *

8. pos(INSYS, 25) ?-0.5797 ? ? 0.2903 ?-1.997 ? 0.0458 *

因此，對于很小的值，原始斜率并不重要，但是在15以上時，它會變得明顯為正。而在25以上，又發(fā)生了重大變化。我們可以對其進(jìn)行繪圖以查看發(fā)生了什么

1. 


2. plot(u,v,type="l")

3. points(INSYS,PRONO)

4. abline(v=c(5,15,25,55)

使用bs（）線性樣條曲線

使用GAM模型，情況略有不同。我們將在這里使用所謂的?b樣條曲線，

我們可以用邊界結(jié)點(diǎn)（5,55）和結(jié) {15,25}定義樣條函數(shù)

1. 


2. B = bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=1)

3. matplot(x,B,type="l",lty=1,lwd=2,col=clr6)

如我們所見，此處定義的函數(shù)與之前的函數(shù)不同，但是在每個段（5,15）（15,25）和（25,55）。但是這些函數(shù)（兩組函數(shù)）的線性組合將生成相同的空間。換個角度說，對輸出的解釋會不同，預(yù)測應(yīng)該是一樣的。

1. 


2. 


3. Coefficients:

4. Estimate Std. Error z value Pr(>|z|)

5. (Intercept) ? ?-0.9863 ? ? 2.0555 ?-0.480 ? 0.6314

6. bs(INSYS,..)1 ?-1.7507 ? ? 2.5262 ?-0.693 ? 0.4883

7. bs(INSYS,..)2 ? 4.3989 ? ? 2.0619 ? 2.133 ? 0.0329 *

8. bs(INSYS,..)3 ? 5.4572 ? ? 5.4146 ? 1.008 ? 0.3135

觀察到像以前一樣存在三個系數(shù)，但是這里的解釋更加復(fù)雜了

但是，預(yù)測結(jié)果很好。

分段二次樣條

讓我們再往前走一步...我們是否也可以具有導(dǎo)數(shù)的連續(xù)性？考慮拋物線函數(shù)，不要對

和

進(jìn)行分解，考慮對

和

進(jìn)行分解。

1. Coefficients:

2. Estimate Std. Error z value Pr(>|z|)

3. (Intercept) ? ? ?29.9842 ? ?15.2368 ? 1.968 ? 0.0491 *

4. poly(INSYS, 2)1 408.7851 ? 202.4194 ? 2.019 ? 0.0434 *

5. poly(INSYS, 2)2 199.1628 ? 101.5892 ? 1.960 ? 0.0499 *

6. pos2(INSYS, 15) ?-0.2281 ? ? 0.1264 ?-1.805 ? 0.0712 .

7. pos2(INSYS, 25) ? 0.0439 ? ? 0.0805 ? 0.545 ? 0.5855

不出所料，這里有五個系數(shù)：截距和拋物線函數(shù)的三個參數(shù)，然后是中間兩個附加項(xiàng)–此處（15,25）–以及右側(cè)的部分。當(dāng)然，對于每個部分，只有一個自由度，因?yàn)槲覀冇幸粋€拋物線函數(shù)（三個系數(shù)），但是有兩個約束（連續(xù)性和一階導(dǎo)數(shù)的連續(xù)性）。

在圖上，我們得到以下內(nèi)容

使用bs（）二次樣條

當(dāng)然，我們可以使用R函數(shù)執(zhí)行相同的操作。但是和以前一樣，這里的函數(shù)有所不同

1. 


2. matplot(x,B,type="l",col=clr6)

如果我們運(yùn)行R代碼，得到

1. glm(y~bs(INSYS knots=c(15,25),

2. Boundary.knots=c(5,55),degre=2)

3. 


4. Coefficients:

5. Estimate Std. Error z value Pr(>|z|)

6. (Intercept) ? ? ? 7.186 ? ? ?5.261 ? 1.366 ? 0.1720

7. bs(INSYS, ..)1 ?-14.656 ? ? ?7.923 ?-1.850 ? 0.0643 .

8. bs(INSYS, ..)2 ? -5.692 ? ? ?4.638 ?-1.227 ? 0.2198

9. bs(INSYS, ..)3 ? -2.454 ? ? ?8.780 ?-0.279 ? 0.7799

10. bs(INSYS, ..)4 ? ?6.429 ? ? 41.675 ? 0.154 ? 0.8774

預(yù)測是完全相同的

1. 


2. plot(u,v,ylim=0:1,type="l",col="red")

三次樣條

我們可以使用三次樣條曲線。我們將考慮對

進(jìn)行分解，得到時間連續(xù)性，以及前兩個導(dǎo)數(shù)的連續(xù)性。如果我們使用bs函數(shù)，則如下

1. matplot(x,B,type="l",lwd=2,col=clr6,lty=1

2. abline(v=c(5,15,25,55),lty=2)

現(xiàn)在的預(yù)測將是

1. bs(x,knots=c(15,25),

2. Boundary.knots=c(5,55),degre=3

?

結(jié)的位置

在許多應(yīng)用程序中，我們不想指定結(jié)的位置。我們只想說（三個）中間結(jié)?？梢允褂?/p>

bs(x,degree=1,df=4)

可以查看

1. 


2. bs(x, degree = 1L, knots = c(15.8, 21.4, 27.15),

3. Boundary.knots = c(8.7, 54), intercept = FALSE)

它為我們提供了邊界結(jié)的位置（樣本中的最小值和最大值），也為我們提供了三個中間結(jié)。觀察到實(shí)際上，這五個值只是（經(jīng)驗(yàn)）分位數(shù)

1. quantile( ,(0:4)/4)

2. 0% ? 25% ? 50% ? 75% ?100%

3. 8.70 15.80 21.40 27.15 54.00

如果我們繪制預(yù)測，我們得到

plot(u,v,ylim=0:1,type="l",col="red",lwd=2)

如果我們回到logit變換之前的計(jì)算，我們清楚地看到斷點(diǎn)是不同的分位數(shù)?

1. 


2. plot(x,y,type="l",col="red",lwd=2)

3. abline(v=quantile(my ,(0:4)/4),lty=2)

如果我們沒有指定，則不會得到任何結(jié)…

1. 


2. bs(x, degree = 2L, knots = numeric(0),

3. Boundary.knots = c(8.7,54), intercept = FALSE)

如果我們看一下預(yù)測?

predict(reg,newdata=data.frame(u),type="response")

實(shí)際上，這和二次多項(xiàng)式回歸是一樣的（如預(yù)期的那樣）

相加模型?

現(xiàn)在考慮第二個數(shù)據(jù)集，包含兩個變量。這里考慮一個模型

然后我們用glm函數(shù)來實(shí)現(xiàn)相加模型的思想。

1. glm(y~bs(x1,degree=1,df=3)+bs(x2,degree=1,df=3), family=binomial(link =

2. v = outer(u,u,p)

3. image(u,u,v, ",col=clr10,breaks=(0:10)/10)

現(xiàn)在，我們能夠得到一個“完美”的模型，所以，結(jié)果似乎不再連續(xù)

persp(u,u,v,theta=20,phi=40,col="green"

當(dāng)然，它是分段線性的，有超平面，有些幾乎是垂直的。
我們也可以考慮分段二次函數(shù)??

1. 


2. 


3. contour(u,u,v,levels = .5,add=TRUE)

有趣的是，我們現(xiàn)在有兩個“完美”的模型，白點(diǎn)和黑點(diǎn)的區(qū)域不同。?

在R中，可以使用mgcv包來運(yùn)行g(shù)am回歸。它用于廣義相加模型，但這里只有一個變量，所以實(shí)際上很難看到“可加”部分，可以參考其他GAM文章。

最受歡迎的見解

1.R語言多元Logistic邏輯回歸 應(yīng)用案例

2.面板平滑轉(zhuǎn)移回歸(PSTR)分析案例實(shí)現(xiàn)

3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）

4.R語言泊松Poisson回歸模型分析案例

5.R語言回歸中的Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn)

6.r語言中對LASSO回歸，Ridge嶺回歸和Elastic Net模型實(shí)現(xiàn)

7.在R語言中實(shí)現(xiàn)Logistic邏輯回歸

8.python用線性回歸預(yù)測股票價格

9.R語言如何在生存分析與Cox回歸中計(jì)算IDI，NRI指標(biāo)