logistic回歸和probit回歸預(yù)測公司被ST的概率

1.適合閱讀人群：
知道以下知識點(diǎn)：盒狀圖、假設(shè)檢驗(yàn)、邏輯回歸的理論、probit的理論、看過回歸分析，了解AIC和BIC判別準(zhǔn)則、能自己跑R語言程序
2.本文目的：用R語言演示一個(gè)相對完整的邏輯回歸和probit回歸建模過程，同時(shí)讓自己復(fù)習(xí)一遍在學(xué)校時(shí)學(xué)的知識，記載下來，以后經(jīng)常翻閱。
3.本文不涉及的部分：
（1）邏輯回歸和probit回歸參數(shù)估計(jì)的公式推導(dǎo)，在下一篇寫；
（2）由ROC曲線帶來的閾值選擇，在下下一篇寫；
（3）本文用的數(shù)據(jù)取自王漢生老師《應(yīng)用商務(wù)統(tǒng)計(jì)分析》第四章里的數(shù)據(jù)，直接描述性分析和建模，沒有涉及到數(shù)據(jù)預(yù)處理。
4.廢話少說，上程序：
#適合人群：知道以下知識點(diǎn)：盒狀圖、假設(shè)檢驗(yàn)、邏輯回歸的理論、probit的理論、看過回歸分析，了解AIC和BIC判別準(zhǔn)則、能讀R語言程序
1.#########讀入數(shù)據(jù)##############
a=read.csv("C:/Users/Thinkpad/Desktop/ST.csv",header=T)  
a1=a[a$year==1999,-1]                           #訓(xùn)練集   
a2=a[a$year==2000,-1]                           #測試集       
a1[c(1:5),]
2.####初步描述性分析######
boxplot(ARA~ST,data=a1,main="ARA")                   #畫出各變量與ST的盒狀圖，初步查看因變量單獨(dú)和各個(gè)解釋性變量的關(guān)系
par(mfrow=c(3,2))                                #只是初步的描述性分析，沒有控制其他因素的影響，沒有經(jīng)過嚴(yán)格的統(tǒng)計(jì)檢驗(yàn)
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)          ####邏輯回歸時(shí)：計(jì)算模型的整體顯著性水平#####          
glm1.a=glm(ST~ARA+ASSET+ATO+GROWTH+LEV+ROA+SHARE,           #結(jié)果為7.4e-05，說明模型整體高度顯著，也就是說所考慮的7個(gè)解釋性變量中，至少有一個(gè)與因變量有關(guān)，具體哪一個(gè)不知道
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回歸時(shí)：計(jì)算模型的整體顯著性水平#####        
glm1.b=glm(ST~ARA+ASSET+ATO+GROWTH+LEV+ROA+SHARE,           #和邏輯回歸結(jié)果一樣，顯著
family=binomial(link=probit),data=a1)                      
anova(glm0.b,glm1.b)                               
1-pchisq(31.702,7)                             
 
 
                                                  ####看看是哪個(gè)自變量對因變量有影響#####
Anova(glm1.a,type="III")                              #對模型glm1.a做三型方差分析
summary(glm1.a)                                
Anova(glm1.b,type="III")                              #對模型glm1.b做三型方差分析
summary(glm1.b)                                
 
 
3.#######模型選擇時(shí)要解決的問題：（1）選哪個(gè)模型；（2）選哪個(gè)閾值。
#######其中選6個(gè)中的哪個(gè)模型用ROC曲線確定（里面涉及到兩個(gè)指標(biāo)：TPR，F(xiàn)PR。至于為什么選擇用這兩個(gè)指標(biāo)來衡量模型的好壞，請往下看，下面會解釋，別著急），選擇ROC曲線最上面的那條線所對應(yīng)的模型。
#######模型確定之后，選取閾值可以根據(jù)ROC曲線和實(shí)際業(yè)務(wù)確定。（這里還需要查資料，至于什么ROC曲線，別急，繼續(xù)向下看）
#######6個(gè)模型：邏輯回歸的全模型，邏輯回歸的AIC模型，邏輯回歸的BIC模型，probit回歸的全模型，probit回歸的AIC模型，probit回歸的BIC模型，
 
 
                                                            #我們先隨便選兩個(gè)模型感受一個(gè)AIC和BIC值
AIC(glm0.a)                                 #計(jì)算邏輯回歸方法時(shí)，空模型glm0.a的AIC取值
AIC(glm1.a)                                 #計(jì)算邏輯回歸方法時(shí)，全模型glm1.a的AIC取值
AIC(glm0.a,k=log(length(a1[,1])))                   #計(jì)算邏輯回歸方法時(shí)，空模型glm0.a的BIC取值
AIC(glm1.a,k=log(length(a1[,1])))                   #計(jì)算邏輯回歸方法時(shí)，全模型glm1.a的BIC取值
 
 
 
                                                             #上面只是比較了兩個(gè)模型的AIC值，BIC值，我們有7個(gè)解釋變量，一共會有128個(gè)不同模型，理論上說需要對這128個(gè)模型逐一研究，并選擇最有模型，在R中
                                                             #我們可以自動的、盡量多的根據(jù)AIC搜索最優(yōu)模型
logit.aic=step(glm1.a,trace=0)                     #根據(jù)AIC準(zhǔn)則選擇邏輯回歸最優(yōu)模型
summary(logit.aic)
 
                                 
n=length(a1[,1])                                               #根據(jù)BIC準(zhǔn)則選擇邏輯回歸最優(yōu)模型###                           
logit.bic=step(glm1.a,k=log(n),trace=0)                    
summary(logit.bic)
                                                              #上面AIC和BIC的結(jié)果有點(diǎn)差別，可以理解為AIC三個(gè)結(jié)果都很重要，而其中的ARA極其重要，BIC選擇的模型更簡單
                                                              #AIC選擇的模型的預(yù)測精度似乎更好，我們老師當(dāng)時(shí)也說要用AIC準(zhǔn)則選模型
 
 
probit.aic=step(glm1.b,trace=0)                     #根據(jù)AIC準(zhǔn)則選擇probit回歸最優(yōu)模型，并賦值給probit.aic
summary(probit.aic)                            
probit.bic=step(glm1.b,k=log(n),trace=0)                  #根據(jù)bIC準(zhǔn)則選擇probit回歸最優(yōu)模型，并賦值給probit.bic
summary(probit.bic)                            
 
 
##############畫出6個(gè)模型的ROC曲線來確定最終選哪一個(gè)模型################
p=matrix(0,length(a2[,1]),6)                                       #生成矩陣，用于存儲各模型的預(yù)測值
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)]))                        #計(jì)算預(yù)測得到的概率
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)])                                        #計(jì)算預(yù)測得到的概率
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個(gè)模型的預(yù)測概率
    for(i in 1:ngrids){
        p0=i/ngrids                                           #選取閾值
        ST.hat=1*(prob>p0)                                   #根據(jù)閾值生成預(yù)測值
        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個(gè)模型的TPR與FPR的散點(diǎn)圖
}
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.#########預(yù)測與評估,由ROC曲線，我們這里選擇基于AIC準(zhǔn)則的邏輯回歸模型，閾值選擇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)])                             
####對于每個(gè)個(gè)體，最終的預(yù)測結(jié)果為
a2$ST.pred
####TPR=59.57%,FPR=23.89%
TPR=28/(28+19)
FPR=167/(532+167)
 
 
#####################有一定基礎(chǔ)的到這里就可以結(jié)束啦，感興趣的還可以向下看##########################
 
##########下面我們隨便選幾個(gè)模型，來解釋下為什么要使用TPR和FPR這兩個(gè)指標(biāo)衡量模型的精度，然后畫出ROC曲線，提供邏輯回歸全模型時(shí)，在眾多不同的FPR下的TPR取值######################################
summary(glm1.a)
p=predict(glm1.a,a2)                                #利用邏輯回歸的全模型glm1.a對數(shù)據(jù)a2進(jìn)行預(yù)測
p=exp(p)/(1+exp(p))                             #計(jì)算預(yù)測得到的概率
a2$ST.pred=1*(p>0.5)                             #以0.5為閾值生成預(yù)測值
table(a2[,c(8,9)])                             
 
                                                                   ###從結(jié)果看來，預(yù)測精度699/746=93.7%，沒有正確預(yù)測一家ST公司
                                                                   #####說明不能用總體精度來衡量預(yù)測的好壞，我們有可能犯兩類錯(cuò)誤：（1）把真實(shí)的ST公司預(yù)測為0；（2）把真實(shí)的非ST公司預(yù)測為1。由于我們關(guān)心的是找出那些ST公司
                                                                   #####，可以通過下面兩個(gè)指標(biāo)來度量上面兩種錯(cuò)誤
                                                                   #####TPR：把真實(shí)的ST公司正確地預(yù)測為ST=1的概率；
                                                                   #####FPR：把真實(shí)的非ST公司錯(cuò)誤地預(yù)測為ST=1的概率
                                                                   #####上面預(yù)測TPR=0（很糟糕），F(xiàn)PR=0（非常好），下面我們把閾值改為0試試結(jié)果
 
a2$ST.pred=1*(p>0)                               #以0為閾值生成預(yù)測值，TPR=100%（非常好），F(xiàn)PR=100%（很糟糕）
table(a2[,c(8,9)])                                                                                                
                                                                  ######由結(jié)果可知這兩個(gè)指標(biāo)的取值是魚和熊掌不可兼得
 
 
                                                                    
a2$ST.pred=1*(p>0.05)                                #以0.05為閾值生成預(yù)測值
table(a2[,c(8,9)])                              #計(jì)算預(yù)測值與真實(shí)值的2維頻數(shù)表
 
 
######上面一直說了ROC曲線，這里開始解釋ROC曲線是何方神圣，上面說了FPR和TPR是魚和熊掌不可兼得，那么現(xiàn)在我們便以FPR為橫坐標(biāo)，TPR為縱坐標(biāo)，畫出他們的曲線，看看他們究竟是什么關(guān)系，而這個(gè)曲線的名字就是ROC曲線
#########下面為了得到全面的分析，我們寫了循環(huán)，以邏輯回歸的全模型為例，提供在眾多不同的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的散點(diǎn)圖，即ROC曲線
points(c(0,1),c(0,1),type="l",lty=2)                        #添加對角線

5.結(jié)果:

圖1 箱形圖：用來觀察哪個(gè)變量對因變量有影響

圖2 ROC曲線：為了確定選擇哪個(gè)模型以及作為閾值選擇的初步參考

圖3 預(yù)測結(jié)果：1：ST了，0：未被ST

圖4 預(yù)測模型的精度