應用logistic回歸模型對SoftDrink數據集進行分析

實驗目的：

應用logistic回歸模型對SoftDrink數據集進行分析：

（1）將SoftDrink數據集分成測試集和訓練樣本集

（2）利用訓練樣本集建立logistic模型

（3）對所建立logistic回歸模型進行顯著性檢驗：（1）方差分析（2）計算麥克法登偽R^2,並對結果進行解釋

（4）對所建立的模型進行回歸係數顯著性檢驗，並根據結果挑選顯著性變數進行重新建模

（5）嘗試用兩種不同的方法對數據進行過散布診斷，並對結果進行解釋

（6）嘗試用glm.diag.plots（）函數對數據進行異常值診斷，如存在異常值，請給予相應的處理措施後，然後重新建模

(7) 嘗試用逐步回歸、向前逐步、向後逐步回歸建立logistic方程，並根據測試集的預測效果選取最優模型，寫出logistic方程，同時針對測試集計算該模型在測試集的預測準確率

二、實驗環境：

R語言

三、實驗數據：SoftDrink數據集

一、前言

在日常學習或工作中經常會使用線性回歸模型對某一事物進行預測，例如預測房價、身高、GDP、學生成績等，發現這些被預測的變數都屬於連續型變數。然而有些情況下，被預測變數可能是二元變數，即成功或失敗、流失或不流失、漲或跌等，對於這類問題，線性回歸將束手無策。這個時候就需要另一種回歸方法進行預測，即Logistic回歸——是一種基於最優化演算法的分類方法。

1、優點：計算代價不高，易於理解和實現。

2、缺點：容易欠擬合，分類精度可能不高。

3、適用數據類型：數值型和標稱型數據。

4、實際應用：Logistic模型主要有三大用途：

1）尋找危險因素，找到某些影響因變數的"壞因素"，一般可以通過優勢比發現危險因素；

2）用於預測，可以預測某種情況發生的概率或可能性大小；

3）用於判別，判斷某個新樣本所屬的類別。

5、與普通線性回歸模型的比較

1）Logistic回歸模型的因變數為二分類變數；

2）該模型的因變數和自變數之間不存在線性關係；

3）一般線性回歸模型中需要假設獨立同分布、方差齊性等，而Logistic回歸模型不需要；

4）Logistic回歸沒有關於自變數分布的假設條件，可以是連續變數、離散變數和虛擬變數；

5）由於因變數和自變數之間不存在線性關係，所以參數(偏回歸係數)使用最大似然估計法計算。

二、R語言實施過程

2.1、繪畫logistic分布的相關圖

首先我們對logistic相關圖進行繪製，並了解他與正態分布的差異點。

2.1.1 Logistic分布和標準正態分布的累計概率

2.1.2繪製Logistic分布和標準正態分布的密度曲線

2.2、logistic回歸

2.2.1、數據預處理

讀入數據之後，我們將Choice變數轉化為因子變數（整數型無法進行）。並建立除品牌外的其他所有解釋變數的logistic回歸方程。

setwd("C:/Users/Administrator/Desktop/2018.04.26 logistic分析")SoftDrink<-read.table(file="SoftDrink.txt",header=TRUE)class(SoftDrink$Choice) #顯示Choice變數的類型 SoftDrink$Choice<-as.factor(SoftDrink$Choice) #將choice變數轉換為因子類型Fit1<-glm(Choice~.-Brand,data=SoftDrink,family=binomial(link="logit"))

2.2.2、建立模型

建立除品牌外的其他所有解釋變數的logistic回歸方程

anova(Fit1,test="Chisq") #依據卡方分布進行回歸方程的顯著性檢驗summary(Fit1)

有上表我們得知：常係數p值為0.00987<0.05顯著不為0，其次calories和Fruits的p值都小於0.05，回歸係數較為顯著。於是我們捨棄其他變數，建立包含熱量和果汁含量的二元方程。

建立包含熱量和果汁含量的二元方程

Fit2<-glm(Choice~Calories+Fruits,data=SoftDrink,family=binomial(link="logit"))anova(Fit2,test="Chisq")summary(Fit2)coef(Fit2)exp(coef(Fit2))

由上表我們可以得知：當顯著性水平a為0.1時，回歸係數顯著性均為顯著。結果表明：當果汁含量不變時，熱量每增加一個單位，導致logit P 平均減少0.017；當熱量不變時，果汁含量每增加一個單位，將導致logit P平均增大0.354。

由於logistic回歸方程的係數含義不直觀，於是我們計算優勢比：

exp(coef(Fit2))

得到：當果汁含量不變時，熱量每增加一個單位，購買私意向的優勢是原來的0.983（優勢比）；當熱量不變時，果汁含量每增加一個單位，購買私意向的優勢是原來的1.425 。可見熱量降低了購買的可能性，而果汁含量增加了購買的可能性。

2.2.3、回歸診斷

Logistic回歸分析的回歸診斷內容與工具與一般線性模型基本相同。再次重點針對異常點的診斷、過散布診斷、擬合優度測度。

Ps：當出現過散布情況並沒有採取相應措施，仍採用極大似然估計法，會導致回歸係數的標準誤被低估，P偏小，進而影響回歸係數的決策。

接下來我們進行一一的分析：

①、異常值診斷

首先我們利用boot包中的glm.diag.plots（）繪圖，通過第四幅圖我們可以直觀的發現異常的觀測點，並通過互動的形式，得出異常點為93、138 。

其次我們採用「探索離群點的統計檢驗方法」利用car包中的outlierTest（）來檢驗。

由此我們可以得知：學生化殘差探測出的離群點為138號觀測，其學生化殘差值為-3.902172。依據近似正態分布計算的概率P值為9.5333e-05。在顯著性為0.05水平下，證明他就是離群點。為了克服異方差的問題，我們將在接下來的步驟對其進行刪除。

#install.packages("boot")library(boot)glm.diag.plots(Fit2,iden=T); 4#install.packages("car")#install.packages("haven")library(car)library(haven)#install.packages("openxlsx")library(openxlsx)#install.packages("readxl")library(readxl)outlierTest(Fit2)

②、過散布診斷

其值為 偏差/自由度。當過散布係數大於1時，成為出現過散布情況。過散布的診斷有兩種常用方法。

方法一：建立二項分布假設的模型，將剩餘的模型偏差除以自由度作為過散布係數。計算散布係數的估計值，小於1表示未出現過散布現象，不用處理；大於1表示出現了過散布情形，需要將參數family=quasibinomial

我們得出值為：0.25>0.05

#方法一：用剩餘的模型偏差除以自由度估計散布係數summary(Fit)$dispersion #顯示默認散布係數Fit$deviance #顯示二元模型的剩餘偏差Fit$df.residual #顯示自由度（樣本量-待估參數）Fit$deviance/Fit$df.residual

方法二：利用軟體包qcc中的qcc.overdispersion.test函數。原假設為散步係數為1 。

我們可以得知，p值=0.316<0.05，拒絕原假設，所以為不存在過散布現象。

#方法二：用軟體包qcc中的qcc.overdispersion.test函數 P265install.packages("qcc")library("qcc")Count.01<-tapply(SoftDrink$Choice,INDEX=SoftDrink$Choice,FUN=length) #計算0和1的觀測個數 qcc.overdispersion.test(c(0,1),type="binomial",size=Count.01)

③、擬合優度測度

Logistic回歸方程的擬合優度測度基於混淆矩陣的預測正確率的測度和麥克法登的偽R^2統計量測度。

方法一：基於混淆矩陣預測正確率

我們可以得知，實際值為0預測值也為0的有130個，預測錯誤的有6個，正確率為95.6%；實際值為1預測值也為1的有147個，預測錯誤的有9個，正確率為94.2%。總正確率為94.9%。

BuyProb<-predict(Fit2,SoftDrink,type="response") #計算預測值為1的概率BuyOrNot<-ifelse(BuyProb>0.5,1,0) #計算預測值，購買概率大於0.5，預測值為1，否則為0 (ConfuseMatrix<-table(SoftDrink$Choice,BuyOrNot))prop.table(ConfuseMatrix,1)*100

方法二：麥克法登的偽R^2

定義為：1-LL/LL0

當R^2越接近1，模型的擬合優度越高，越接近0，模型的擬合優度越低。

我們可以得知，偽R^2=0.82 。可見擬合優度較為理想。

Fit<-glm(Choice~.-Brand-Price-Fat-Age-Vitamin,data=SoftDrink[-138,],family=binomial (link="logit"))anova(Fit)(McR2<-1-anova(Fit)[3,4]/anova(Fit)[1,4])

2.2.4、模型的改進與提升

以上我們建立了logistic回歸模型，再此我們首先通過數據的劃分，更加科學準確的使用訓練樣本進行訓練，而將最終的模型通過測試樣本進行準確性比較。

其次我們可以會發現，這麼多變數我們應該如何選取，總不可能還是和上述一樣，進行人工篩選吧？於是我們使用了幾種篩選變數的方法進行了研究。得出最優的logistic模型。

①、將數據分為訓練樣本與測試樣本，並預測。

set.seed(1)train=sample(1:nrow(SoftDrink),nrow(SoftDrink)*3/4)test=(-train)train_sample=SoftDrink[train,]head(train_sample)test_sample=SoftDrink[test,]dim(train_sample)dim(test_sample)Fit3=glm(Choice~.-Brand,data=train_sample,family=binomial(link="logit"))anova(Fit3,test="Chisq")summary(Fit3)coef(Fit3)exp(coef(Fit3))#計算模型內預測準確率BuyProb_train<-predict(Fit3,train_sample,type="response") #計算預測值為1的概率BuyOrNot<-ifelse(BuyProb_train>0.5,1,0) #計算預測值，購買概率大於0.5，預測值為1，否則為0 (ConfuseMatrix<-table(train_sample$Choice,BuyOrNot))prop.table(ConfuseMatrix,1)*100#計算模型外即測試樣本準確率BuyProb_test<-predict(Fit3,test_sample,type="response") #計算預測值為1的概率BuyOrNot<-ifelse(BuyProb_test>0.5,1,0) #計算預測值，購買概率大於0.5，預測值為1，否則為0 (ConfuseMatrix<-table(test_sample$Choice,BuyOrNot))prop.table(ConfuseMatrix,1)*100

可知預測結果比上述要高。說明模型較為合適。

②、利用逐步回歸篩選變數

Fit_both=step(Fit3,direction="both")anova(Fit_both,test="Chisq")summary(Fit_both)coef(Fit_both)exp(coef(Fit_both))#計算模型內預測準確率BuyProb_train<-predict(Fit_both,train_sample,type="response") #計算預測值為1的概率BuyOrNot<-ifelse(BuyProb_train>0.5,1,0) #計算預測值，購買概率大於0.5，預測值為1，否則為0 (ConfuseMatrix<-table(train_sample$Choice,BuyOrNot))prop.table(ConfuseMatrix,1)*100#計算模型外即測試樣本準確率BuyProb_test<-predict(Fit_both,test_sample,type="response") #計算預測值為1的概率BuyOrNot<-ifelse(BuyProb_test>0.5,1,0) #計算預測值，購買概率大於0.5，預測值為1，否則為0 (ConfuseMatrix<-table(test_sample$Choice,BuyOrNot))prop.table(ConfuseMatrix,1)*100

③、利用向前逐步回歸篩選變數

Fit_forward=step(Fit3,direction="forward")anova(Fit_forward,test="Chisq")summary(Fit_forward)coef(Fit_forward)exp(coef(Fit_forward))#計算模型內預測準確率BuyProb_train<-predict(Fit_forward,train_sample,type="response") #計算預測值為1的概率BuyOrNot<-ifelse(BuyProb_train>0.5,1,0) #計算預測值，購買概率大於0.5，預測值為1，否則為0 (ConfuseMatrix<-table(train_sample$Choice,BuyOrNot))prop.table(ConfuseMatrix,1)*100#計算模型外即測試樣本準確率BuyProb_test<-predict(Fit_forward,test_sample,type="response") #計算預測值為1的概率BuyOrNot<-ifelse(BuyProb_test>0.5,1,0) #計算預測值，購買概率大於0.5，預測值為1，否則為0 (ConfuseMatrix<-table(test_sample$Choice,BuyOrNot))prop.table(ConfuseMatrix,1)*100

④、利用向後逐步回歸篩選變數

Fit_backward=step(Fit3,direction="backward")anova(Fit_backward,test="Chisq")summary(Fit_backward)coef(Fit_backward)exp(coef(Fit_backward))#計算模型內預測準確率BuyProb_train<-predict(Fit_backward,train_sample,type="response") #計算預測值為1的概率BuyOrNot<-ifelse(BuyProb_train>0.5,1,0) #計算預測值，購買概率大於0.5，預測值為1，否則為0 (ConfuseMatrix<-table(train_sample$Choice,BuyOrNot))prop.table(ConfuseMatrix,1)*100#計算模型外即測試樣本準確率BuyProb_test<-predict(Fit_backward,test_sample,type="response") #計算預測值為1的概率BuyOrNot<-ifelse(BuyProb_test>0.5,1,0) #計算預測值，購買概率大於0.5，預測值為1，否則為0 (ConfuseMatrix<-table(test_sample$Choice,BuyOrNot))prop.table(ConfuseMatrix,1)*100

通過上述三種變數選取的方法的混淆矩陣我們得知，三種篩選變數的出的結論大致相同。由於我們的變數之間的獨立性很強，並且變數的數量也不是很多。所以變數篩選的方法，再此顯得有些多餘了。但如果推廣到其他問題或者其他數據集，可能發揮巨大的的作用。

三、總結

再次我們對logistic回歸模型進行了詳細的闡述以及通過R語言將其實現。並通過假設檢驗、回歸診斷、模型的提升等方法，對模型進行檢驗與提升。