共计 5660 个字符,预计需要花费 15 分钟才能阅读完成。
原文链接: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")
三次样条
咱们能够应用三次样条曲线。咱们将思考对进行合成,失去工夫连续性,以及前两个导数的连续性。如果咱们应用 bs 函数,则如下
matplot(x,B,type="l",lwd=2,col=clr6,lty=1
abline(v=c(5,15,25,55),lty=2)
当初的预测将是
bs(x,knots=c(15,25),
Boundary.knots=c(5,55),degre=3
\
结的地位
在许多应用程序中,咱们不想指定结的地位。咱们只想说(三个)两头结。能够应用
bs(x,degree=1,df=4)
能够查看
bs(x, degree = 1L, knots = c(15.8, 21.4, 27.15),
Boundary.knots = c(8.7, 54), intercept = FALSE)
它为咱们提供了边界结的地位(样本中的最小值和最大值),也为咱们提供了三个两头结。察看到实际上,这五个值只是(教训)分位数
quantile(,(0:4)/4)
0% 25% 50% 75% 100%
8.70 15.80 21.40 27.15 54.00
如果咱们绘制预测,咱们失去
plot(u,v,ylim=0:1,type="l",col="red",lwd=2)
如果咱们回到 logit 变换之前的计算,咱们分明地看到断点是不同的分位数
plot(x,y,type="l",col="red",lwd=2)
abline(v=quantile(my ,(0:4)/4),lty=2)
如果咱们没有指定,则不会失去任何结…
bs(x, degree = 2L, knots = numeric(0),
Boundary.knots = c(8.7,54), intercept = FALSE)
如果咱们看一下预测
predict(reg,newdata=data.frame(u),type="response")
实际上,这和二次多项式回归是一样的(如预期的那样)
相加模型
当初思考第二个数据集,蕴含两个变量。这里思考一个模型
而后咱们用 glm 函数来实现相加模型的思维。
glm(y~bs(x1,degree=1,df=3)+bs(x2,degree=1,df=3), family=binomial(link =
v = outer(u,u,p)
image(u,u,v, ",col=clr10,breaks=(0:10)/10)
当初,咱们可能失去一个“完满”的模型,所以,后果仿佛不再间断
persp(u,u,v,theta=20,phi=40,col="green"
当然,它是分段线性的,有超平面,有些简直是垂直的。\
咱们也能够思考分段二次函数
contour(u,u,v,levels = .5,add=TRUE)
乏味的是,咱们当初有两个“完满”的模型,白点和黑点的区域不同。
在 R 中,能够应用 mgcv 包来运行 gam 回归。它用于狭义相加模型,但这里只有一个变量,所以实际上很难看到“可加”局部,能够参考其余 GAM 文章。
最受欢迎的见解
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 指标