关于数据挖掘:R语言用普通最小二乘OLS广义相加模型GAM-样条函数进行逻辑回归LOGISTIC分类附代码数据

29次阅读

共计 5660 个字符，预计需要花费 15 分钟才能阅读完成。

最近咱们被客户要求撰写对于回归的钻研报告，包含一些图形和统计输入。

本文咱们对 [逻辑回归和样条曲线进行]() 介绍。

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))

输入在这里

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

当然，能够应用 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(&gt;|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)

应用 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(&gt;|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(&gt;|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）–以及右侧的局部。当然，对于每个局部，只有一个自由度，因为咱们有一个抛物线函数（三个系数），然而有两个束缚（连续性和一阶导数的连续性）。

在图上，咱们失去以下内容

当然，咱们能够应用 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(&gt;|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 指标

正文完

数据挖掘

发表至：数据挖掘

2023-09-25

0

关于数据挖掘:R语言分布滞后线性和非线性模型DLMs和DLNMs分析时间序列数据附代码数据

关于数据挖掘:R语言广义加性混合模型GAMM分析长沙气象因子空气污染PM25浓度显著性检验逐日变化可视化

关于数据挖掘:Matlab中的偏最小二乘法PLS回归模型离群点检测和变量选择附代码数据

关于数据挖掘:ARIMA模型ARIMAX模型预测冰淇淋消费时间序列数据附代码数据

关于数据挖掘:R语言文本挖掘kmeans聚类分析上海玛雅水公园景区五一假期评论词云可视化附代码数据

关于数据挖掘:R语言用普通最小二乘OLS广义相加模型GAM-样条函数进行逻辑回归LOGISTIC分类附代码数据

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

原文出处：拓端数据部落公众号

（对数）似然函数

牛顿（或费舍尔）算法

规范最小二乘

规范逻辑回归 glm 函数：

可视化

分段线性样条函数

应用 bs（）线性样条曲线

分段二次样条

应用 bs（）二次样条

三次样条

结的地位

相加模型

Just My Socks（注册教程内含优惠码）

关于数据挖掘:R语言用普通最小二乘OLS广义相加模型GAM-样条函数进行逻辑回归LOGISTIC分类附代码数据

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

原文出处：拓端数据部落公众号

（对数）似然函数

牛顿（或费舍尔）算法

规范最小二乘

规范逻辑回归 glm 函数：

可视化

分段线性样条函数

应用 bs（）线性样条曲线

分段二次样条

应用 bs（）二次样条

三次样条

结的地位

相加模型

Just My Socks（注册教程 内含优惠码）

Just My Socks（注册教程内含优惠码）