logistic回归和probit回归预测公司被ST的概率(应用)

来源:转载

1.适合阅读人群:

知道以下知识点:盒状图、假设检验、逻辑回归的理论、probit的理论、看过回归分析,了解AIC和BIC判别准则、能自己跑R语言程序

2.本文目的:用R语言演示一个相对完整的逻辑回归和probit回归建模过程,同时让自己复习一遍在学校时学的知识,记载下来,以后经常翻阅。

3.本文不涉及的部分:

(1)逻辑回归和probit回归参数估计的公式推导,在下一篇写;

(2)由ROC曲线带来的阈值选择,在下下一篇写;

(3)本文用的数据取自王汉生老师《应用商务统计分析》第四章里的数据,直接描述性分析和建模,没有涉及到数据预处理。

4.废话少说,上程序:

#适合人群:知道以下知识点:盒状图、假设检验、逻辑回归的理论、probit的理论、看过回归分析,了解AIC和BIC判别准则、能读R语言程序

1.#########读入数据##############

a=read.csv("C:/Users/Thinkpad/Desktop/ST.csv",header=T)

a1=a[a$year==1999,-1] #训练集

a2=a[a$year==2000,-1] #测试集

a1[c(1:5),]

2.####初步描述性分析######

boxplot(ARA~ST,data=a1,main="ARA") #画出各变量与ST的盒状图,初步查看因变量单独和各个解释性变量的关系

par(mfrow=c(3,2)) #只是初步的描述性分析,没有控制其他因素的影响,没有经过严格的统计检验

boxplot(ASSET~ST,data=a1,main="ASSET")

boxplot(ATO~ST,data=a1,main="ATO")

boxplot(GROWTH~ST,data=a1,main="GROWTH")

boxplot(LEV~ST,data=a1,main="LEV")

boxplot(ROA~ST,data=a1,main="ROA")

boxplot(SHARE~ST,data=a1,main="SHARE")

par(mfrow=c(1,1))

glm0.a=glm(ST~1,family=binomial(link=logit),data=a1) ####逻辑回归时:计算模型的整体显著性水平#####

glm1.a=glm(ST~ARA+ASSET+ATO+GROWTH+LEV+ROA+SHARE, #结果为7.4e-05,说明模型整体高度显著,也就是说所考虑的7个解释性变量中,至少有一个与因变量有关,具体哪一个不知道

family=binomial(link=logit),data=a1)

anova(glm0.a,glm1.a)

1-pchisq(30.565,7)

glm0.b=glm(ST~1,family=binomial(link=probit),data=a1)####probit回归时:计算模型的整体显著性水平#####

glm1.b=glm(ST~ARA+ASSET+ATO+GROWTH+LEV+ROA+SHARE, #和逻辑回归结果一样,显著

family=binomial(link=probit),data=a1)

anova(glm0.b,glm1.b)

1-pchisq(31.702,7) ####看看是哪个自变量对因变量有影响#####

Anova(glm1.a,type="III") #对模型glm1.a做三型方差分析

summary(glm1.a)

Anova(glm1.b,type="III") #对模型glm1.b做三型方差分析

summary(glm1.b)

3.#######模型选择时要解决的问题:(1)选哪个模型;(2)选哪个阈值。

#######其中选6个中的哪个模型用ROC曲线确定(里面涉及到两个指标:TPR,FPR。至于为什么选择用这两个指标来衡量模型的好坏,请往下看,下面会解释,别着急),选择ROC曲线最上面的那条线所对应的模型。

#######模型确定之后,选取阈值可以根据ROC曲线和实际业务确定。(这里还需要查资料,至于什么ROC曲线,别急,继续向下看)

#######6个模型:逻辑回归的全模型,逻辑回归的AIC模型,逻辑回归的BIC模型,probit回归的全模型,probit回归的AIC模型,probit回归的BIC模型,

#我们先随便选两个模型感受一个AIC和BIC值

AIC(glm0.a)#计算逻辑回归方法时,空模型glm0.a的AIC取值

AIC(glm1.a)#计算逻辑回归方法时,全模型glm1.a的AIC取值

AIC(glm0.a,k=log(length(a1[,1])))#计算逻辑回归方法时,空模型glm0.a的BIC取值

AIC(glm1.a,k=log(length(a1[,1])))#计算逻辑回归方法时,全模型glm1.a的BIC取值

#上面只是比较了两个模型的AIC值,BIC值,我们有7个解释变量,一共会有128个不同模型,理论上说需要对这128个模型逐一研究,并选择最有模型,在R中

#我们可以自动的、尽量多的根据AIC搜索最优模型

logit.aic=step(glm1.a,trace=0) #根据AIC准则选择逻辑回归最优模型

summary(logit.aic)

n=length(a1[,1]) #根据BIC准则选择逻辑回归最优模型###

logit.bic=step(glm1.a,k=log(n),trace=0)

summary(logit.bic)

#上面AIC和BIC的结果有点差别,可以理解为AIC三个结果都很重要,而其中的ARA极其重要,BIC选择的模型更简单

#AIC选择的模型的预测精度似乎更好,我们老师当时也说要用AIC准则选模型

probit.aic=step(glm1.b,trace=0) #根据AIC准则选择probit回归最优模型,并赋值给probit.aic

summary(probit.aic)

probit.bic=step(glm1.b,k=log(n),trace=0) #根据bIC准则选择probit回归最优模型,并赋值给probit.bic

summary(probit.bic)

##############画出6个模型的ROC曲线来确定最终选哪一个模型################

p=matrix(0,length(a2[,1]),6) #生成矩阵,用于存储各模型的预测值

p[,1]=predict(glm1.a,a2)

p[,2]=predict(logit.aic,a2)

p[,3]=predict(logit.bic,a2)

p[,c(1:3)]=exp(p[,c(1:3)])/(1+exp(p[,c(1:3)])) #计算预测得到的概率

p[,4]=predict(glm1.b,a2)

p[,5]=predict(probit.aic,a2)

p[,6]=predict(probit.bic,a2)

p[,c(4:6)]=pnorm(p[,c(4:6)]) #计算预测得到的概率

plot(c(0,1),c(0,1),type="l",main="FPR vs. TPR",xlab="FPR",ylab="TPR")#画图,生成基本框架

FPR=rep(0,ngrids)

TPR=rep(0,ngrids)

for(k in 1:6){prob=p[,k] #取出p中第K列的值,即第K个模型的预测概率for(i in 1:ngrids){p0=i/ngrids #选取阈值ST.hat=1*(prob>p0) #根据阈值生成预测值FPR[i]=sum((1-ST.true)*ST.hat)/sum(1-ST.true)TPR[i]=sum(ST.true*ST.hat)/sum(ST.true)}points(FPR,TPR,type="b",col=k,lty=k,pch=k) #向图上添加第k个模型的TPR与FPR的散点图

}

legend(0.6,0.5,c("LOGIT FULL MODEL","LOGIT AIC MODEL",

"LOGIT BIC MODEL","PROBIT FULL MODEL","PROBIT AIC MODEL",

"PROBIT BIC MODEL"),lty=c(1:6),col=c(1:6),pch=c(1:6))

4.#########预测与评估,由ROC曲线,我们这里选择基于AIC准则的逻辑回归模型,阈值选择0.05,这块的选择还需要再查阅资料确定###########

p=predict(logit.aic,a2)

p=exp(p)/(1+exp(p))

a2$ST.pred=1*(p>0.05)

table(a2[,c(8,9)])

####对于每个个体,最终的预测结果为

a2$ST.pred

####TPR=59.57%,FPR=23.89%

TPR=28/(28+19)

FPR=167/(532+167)

#####################有一定基础的到这里就可以结束啦,感兴趣的还可以向下看##########################

##########下面我们随便选几个模型,来解释下为什么要使用TPR和FPR这两个指标衡量模型的精度,然后画出ROC曲线,提供逻辑回归全模型时,在众多不同的FPR下的TPR取值######################################

summary(glm1.a)

p=predict(glm1.a,a2)#利用逻辑回归的全模型glm1.a对数据a2进行预测

p=exp(p)/(1+exp(p))#计算预测得到的概率

a2$ST.pred=1*(p>0.5)#以0.5为阈值生成预测值

table(a2[,c(8,9)])

###从结果看来,预测精度699/746=93.7%,没有正确预测一家ST公司

#####说明不能用总体精度来衡量预测的好坏,我们有可能犯两类错误:(1)把真实的ST公司预测为0;(2)把真实的非ST公司预测为1。由于我们关心的是找出那些ST公司

#####,可以通过下面两个指标来度量上面两种错误

#####TPR:把真实的ST公司正确地预测为ST=1的概率;

#####FPR:把真实的非ST公司错误地预测为ST=1的概率

#####上面预测TPR=0(很糟糕),FPR=0(非常好),下面我们把阈值改为0试试结果

a2$ST.pred=1*(p>0)#以0为阈值生成预测值,TPR=100%(非常好),FPR=100%(很糟糕)

table(a2[,c(8,9)])

######由结果可知这两个指标的取值是鱼和熊掌不可兼得

a2$ST.pred=1*(p>0.05)#以0.05为阈值生成预测值

table(a2[,c(8,9)])#计算预测值与真实值的2维频数表

######上面一直说了ROC曲线,这里开始解释ROC曲线是何方神圣,上面说了FPR和TPR是鱼和熊掌不可兼得,那么现在我们便以FPR为横坐标,TPR为纵坐标,画出他们的曲线,看看他们究竟是什么关系,而这个曲线的名字就是ROC曲线

#########下面为了得到全面的分析,我们写了循环,以逻辑回归的全模型为例,提供在众多不同的FPR下的TPR取值

ngrids=100

TPR=rep(0,ngrids)

FPR=rep(0,ngrids)

p0=rep(0,ngrids)

for(i in 1:ngrids){

p0[i]=i/ngrids;#选取阈值p0

ST.true=a2$ST

ST.pred=1*(p>p0[i])

TPR[i]=sum(ST.pred*ST.true)/sum(ST.true)

FPR[i]=sum(ST.pred*(1-ST.true))/sum(1-ST.true)

}

plot(FPR,TPR,type="l",col=2)#画出FPR与TPR的散点图,即ROC曲线

points(c(0,1),c(0,1),type="l",lty=2)#添加对角线

5.结果:

图1 箱形图:用来观察哪个变量对因变量有影响

图2 ROC曲线:为了确定选择哪个模型以及作为阈值选择的初步参考

图3 预测结果:1:ST了,0:未被ST

图4 预测模型的精度

6.参考资料:

王汉生《应用商务统计分析》第四章。

 

 

end!

分享给朋友:
您可能感兴趣的文章:
随机阅读: