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 预测模型的精度