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

本文我们对逻辑回归和样条曲线进行介绍。

logistic回归基于以下假设：给定协变量x，Y具有伯努利分布，

目的是估计参数β。

回想一下，针对该概率使用该函数是

（对数）似然函数

对数似然

其中

。数值方法基于（数值）下降梯度来计算似然函数的最大值。对数似然（负）是以下函数 negLogLik = function(beta){ -sum(-y*log(1 + exp(-(X%*%beta))) - (1-y)*log(1 + exp(X%*%beta))) }

现在，我们需要一个起始点来启动算法

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

在这里，我们得到

logistic_opt$par (Intercept) FRCAR INCAR INSYS 1.656926397 0.045234029 -2.119441743 0.204023835 PRDIA PAPUL PVENT REPUL -0.102420095 0.165823647 -0.081047525 -0.005992238

让我们在这里验证该输出是否有效。例如，如果我们（随机）更改起点的值会怎么样

plot(v_beta) par(mfrow=c(1,2)) hist(v_beta[,1],xlab=names( )[ ]) hist(v_beta[,2],xlab=names( )[2])

这里有个问题。注意，我们不能在这里进行数值优化。我们可以考虑使用其他优化方法

logLikelihoodLogitStable = function(vBeta, mX, vY) { -sum(vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta) + (1-vY)*(-log(1 + exp(mX %*% vBeta))

optimLogitLBFGS = optimx(beta_init, logLikelihoodLogitStable,

最优点

结果不理想。

我们使用的技术基于以下思想，

问题是我的计算机不知道一阶和二阶导数。

可以使用这种计算的函数

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

for(i in 1:num_iter){ grad = (t(X)%*%(logit(X%*%beta) - y)) beta = beta - ginv(H)%*%grad LL[i] = logLik(beta, X, y)

以我们的OLS起点，我们获得

如果我们尝试另一个起点

一些系数非常接近。然后我们尝试其他方法。

牛顿（或费舍尔）算法

在计量经济学教科书里，您可以看到：

beta=as.matrix(lm(Y~0+X)$coefficients for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi))

在这里观察到，我仅使用该算法的十次迭代。

事实是，收敛似乎非常快。而且它相当鲁棒，看看我们改变起点会得到什么

beta=as.matrix(lm(Y~0+X)$coefficients,ncol=1)*runif(8) for(s in 1:9){ pi=exp(X%*%beta[,s])/(1+exp(X%*%beta[,s])) gradient=t(X)%*%(Y-pi) omega=matrix(0,nrow(X),nrow(X));diag(omega)=(pi*(1-pi)) Hessian=-t(X)%*%omega%*%X beta=cbind(beta,beta[,s]-solve(Hessian)%*%gradient)} beta[,8:10]

效果提高了，并且可以使用矩阵的逆获得标准偏差。

标准最小二乘

我们更进一步。我们已经看到想要计算类似

但是实际，这是一个标准的最小二乘问题

这里唯一的问题是权重Δold是未知β的函数。但是实际上，如果我们继续迭代，我们应该能够解决它：给定β，我们得到了权重，并且有了权重，我们可以使用加权的OLS来获取更新的β。这就是迭代最小二乘的想法。

该算法

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

for(s in 1:1000){ omega = diag(nrow(df)) diag(omega) = (p*(1-p))

输出在这里

结果很好，我们在这里也有估计量的标准差

标准逻辑回归glm函数：

当然，可以使用R内置函数

可视化

让我们在第二个数据集上可视化从逻辑回归获得的预测

image(u,u,v ,breaks=(0:10)/10) points(x,y,pch=19 ) points(x,y,pch=c(1,19) contour(u,u,v,levels = .5

这里的水平曲线-或等概率-是线性的，因此该空间被一条直线（或更高维的超平面）一分为二（0和1，生存和死亡，白色和黑色）此外，由于我们是线性模型，因此，如果更改截距（为创建两个类别的阈值），我们将获得平行的另一条直线（或超平面）。

接下来，我们将约会样条曲线以平滑那些连续的协变量。

分段线性样条函数

我们从“简单”回归开始（只有一个解释变量），我们可以想到的最简单的模型来扩展我们上面的线性模型，是考虑一个分段线性函数，它分为两部分。最方便的方法是使用正部函数

（如果该差为正，则为x和s之间的差，否则为0）。如

是以下连续的分段线性函数，在s处划分。

对于较小的x值，线性增加，斜率β1；对于较大的x值，线性减少。因此，β2被解释为斜率的变化。

当然，可以考虑多个结。获得正值的函数如下

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

然后我们可以在回归模型中直接使用它

回归的输出在这里

Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.1109 3.2783 -0.034 0.9730 INSYS -0.1751 0.2526 -0.693 0.4883 pos(INSYS, 15) 0.7900 0.3745 2.109 0.0349 * pos(INSYS, 25) -0.5797 0.2903 -1.997 0.0458 *

因此，对于很小的值，原始斜率并不重要，但是在15以上时，它会变得明显为正。而在25以上，又发生了重大变化。我们可以对其进行绘图以查看发生了什么

plot(u,v,type="l") points(INSYS,PRONO) abline(v=c(5,15,25,55)

使用bs（）线性样条曲线

使用GAM模型，情况略有不同。我们将在这里使用所谓的 b样条曲线，

我们可以用边界结点（5,55）和结 {15,25}定义样条函数

B = bs(x,knots=c(15,25),Boundary.knots=c(5,55),degre=1) matplot(x,B,type="l",lty=1,lwd=2,col=clr6)

如我们所见，此处定义的函数与之前的函数不同，但是在每个段（5,15）（15,25）和（25,55）。但是这些函数（两组函数）的线性组合将生成相同的空间。换个角度说，对输出的解释会不同，预测应该是一样的。

Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.9863 2.0555 -0.480 0.6314 bs(INSYS,..)1 -1.7507 2.5262 -0.693 0.4883 bs(INSYS,..)2 4.3989 2.0619 2.133 0.0329 * bs(INSYS,..)3 5.4572 5.4146 1.008 0.3135

观察到像以前一样存在三个系数，但是这里的解释更加复杂了

但是，预测结果很好。

分段二次样条

让我们再往前走一步...我们是否也可以具有导数的连续性？考虑抛物线函数，不要对

和

进行分解，考虑对

和

进行分解。Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 29.9842 15.2368 1.968 0.0491 * poly(INSYS, 2)1 408.7851 202.4194 2.019 0.0434 * poly(INSYS, 2)2 199.1628 101.5892 1.960 0.0499 * pos2(INSYS, 15) -0.2281 0.1264 -1.805 0.0712 . pos2(INSYS, 25) 0.0439 0.0805 0.545 0.5855

不出所料，这里有五个系数：截距和抛物线函数的三个参数，然后是中间两个附加项–此处（15,25）–以及右侧的部分。当然，对于每个部分，只有一个自由度，因为我们有一个抛物线函数（三个系数），但是有两个约束（连续性和一阶导数的连续性）。

在图上，我们得到以下内容

使用bs（）二次样条

当然，我们可以使用R函数执行相同的操作。但是和以前一样，这里的函数有所不同

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

如果我们运行R代码，得到

glm(y~bs(INSYS knots=c(15,25), Boundary.knots=c(5,55),degre=2) Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 7.186 5.261 1.366 0.1720 bs(INSYS, ..)1 -14.656 7.923 -1.850 0.0643 . bs(INSYS, ..)2 -5.692 4.638 -1.227 0.2198 bs(INSYS, ..)3 -2.454 8.780 -0.279 0.7799 bs(INSYS, ..)4 6.429 41.675 0.154 0.8774

预测是完全相同的

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