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

一、Logistic模型簡介

Logistic迴歸模型公式如下：

xn的情況下，標籤變數y=1時的概率。顯然，該模型是一個非線性模型，具有S型分佈，可見下圖：

#繪製Logistic曲線
x <- seq(from = -10, to = 10, by = 0.01)
y = exp(x)/(1+exp(x))
library(ggplot2)
p <- ggplot(data = , mapping = aes(x = x,y = y))
p + geom_line(colour = 'blue')+ annotate('text', x = 1, y = 0.3, label ='y==e^x / 1-e^x', parse = TRUE)+ ggtitle('Logistic曲線')
 

由於是非線性模型，從而就少了像線性模型那樣的約束，如自變數與因變數具有線性關係、隨機誤差滿足方差齊性等。

以上的模型公式其實是可以變換成線性形式的，只需要一個簡單的logit變換即可，即：

二、模型用途

Logistic模型主要有三大用途：

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

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

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

三、模型應用

下文使用Logistic模型對電信行業的客戶流失資料進行建模，資料來源為R中C50包自帶資料集churnTrain和churnTest。

#使用C50包中自帶的電信行業客戶流失資料
library(C50)
data(churn)
train <- churnTrain
test <- churnTest
str(train)

資料集中包含了20個變數，其中變數洲(state)、國際長途計劃(international_plan)、信箱語音計劃(voice_mail_plan)和是否流失(churn)為因子變數，其餘變數均為數值變數，而且這裡的區域編碼變數(area_code)沒有任何實際意義，故考慮排除該變數。

#剔除無意義的區域編碼變數
train <- churnTrain[,-3]
test <- churnTest[,-3]
#由於模型中，更關心的是流失這個結果(churn=yes)，所以對該因子進行排序
train$churn <- factor(train$churn,levels = c('no','yes'), order = TRUE)
test$churn <- factor(test$churn, ,levels = c('no','yes'), order = TRUE)
#構建Logistic模型
model <- glm(formula = churn ~ ., data =train, family = 'binomial')
summary(model)
 

發現有很多變數並不顯著，故考慮剔除這些不顯著的變數，這裡使用逐步迴歸法進行變數的選擇(需要注意的是，Logistic為非線性模型，迴歸係數是通過極大似然估計方法計算所得)。

#step函式實現逐步迴歸法
model2 <- step(object = model, trace = 0)
summary(model2)

從結果中發現，所有變數的P值均小於0.05，通過顯著性檢驗，保留了相對重要的變數。模型各變數通過顯著性檢驗的同時還需確保整個模型是顯著的，只有這樣才能保證模型是正確的、有意義的，下面對模型進行卡方檢驗。、

#模型的顯著性檢驗
anova(object = model2, test = 'Chisq')

從上圖中可知，隨著變數從第一個到最後一個逐個加入模型，模型最終通過顯著性檢驗，說明由上述這些變數組成的模型是有意義的，並且是正確的。

雖然模型的偏回歸係數和模型均通過顯著性檢驗，但不代表模型能夠非常準確的擬合實際值，這就需要對模型進行擬合優度檢驗，即通過比較模型的預測值與實際值之間的差異情況來進行檢驗。

Logistic迴歸模型的擬合優度檢驗一般使用偏差卡方檢驗、皮爾遜卡方檢驗和HL統計量檢驗三種方法，其中前兩種檢驗適合模型中只有離散的自變數，而後一種適合模型中包含連續的自變數。擬合優度檢驗的原假設為“模型的預測值與實際值不存在差異”。

#模型的擬合優度檢驗
library(sjmisc)
HL_test <- hoslem_gof(x = model)
HL_test

從模型的擬合優度檢驗結果可知，P.value>0.05，該模型接受擬合優度檢驗的原假設，即可以認為實際值與模型的預測值之間比較接近，不存在顯著差異。

以上各項指標均表示模型對電信行業客戶流失資料擬合的比較理想，接下來就用該模型對測試集進行預測，預測一個未知的客戶是否可能流失，從而起到流失預警的作用。

#模型對樣本外資料(測試集)的預測精度
prob <- predict(object = model2, newdata= test, type = 'response')
pred <- ifelse(prob>= 0.5, 'yes','no')
pred <- factor(pred, levels =c('no','yes'), order = TRUE)
f <- table(test$churn, pred)
f

從上圖中我們發現：

1).模型對非流失客戶(no)的預測還是非常準確的（1408/(1408+35)=97.6%）;

2).模型對流失客戶(yes)的預測非常不理想（42/(182+42)=18.8%）

3).模型的整體預測準確率為87.0%((1408+42)/(1408+35+182+42))，還算說得過去。

模型對非流失客戶預測精準，而對流失客戶預測非常差，我認為的可能原因是模型對非平衡資料非常敏感。即構建模型的訓練集中流失客戶為483例，而非流失客戶為2850例，兩者相差非常大。

上文對模型偏回歸係數、模型整體和模型擬合優度進行了顯著性檢驗，結果均表明模型比較理想，同時也對模型的預測精度進行驗證，也說明了模型的整體預測能力比較理想。接下來我們通過另一種視覺化的方法衡量模型的優劣，即ROC曲線，該曲線的橫座標和縱座標各表示1-反例的覆蓋率（正確預測到的正例數/實際正例總數）和正例的覆蓋率（正確預測到的負例數/實際負例總數）。

#繪製ROC曲線
library(pROC)
roc_curve <- roc(test$churn,prob)
names(roc_curve)
x <- 1-roc_curve$specificities
y <- roc_curve$sensitivities
library(ggplot2)
p <- ggplot(data = , mapping = aes(x= x, y = y))
p + geom_line(colour = 'red') +geom_abline(intercept = 0, slope = 1)+ annotate('text', x = 0.4, y = 0.5, label =paste('AUC=',round(roc_curve$auc,2)))+ labs(x = '1-specificities',y = 'sensitivities', title = 'ROC Curve')

通過比較ROC曲線與45°直線可以直觀的反映模型的好壞，但並不能從定量的角度反饋模型好是好到什麼程度或模型差是差到什麼程度。那麼就引申出了AUC的概念，這裡的AUC為ROC曲線和y=x直線之間的面積。在實際應用中，多個模型的比較可以通過面積大小來選擇更佳的模型，選擇標準是AUC越大越好。對於一個模型而言，一般AUC大於0.8就能夠說明模型是比較合理的了。

統計學中，對於單樣本的K-S檢驗就是利用樣本資料來推斷其是否服從某種分佈，對於兩樣本的K-S檢驗主要推測的是兩個樣本是否具有相同的分佈，對於模型的評估，希望正例的累積概率分佈與負例的累積概率分佈存在顯著差異。

所以我們使用K-S統計量刻畫模型的優劣，即使正例與負例的累積概率差達到最大。這是一個定量的判斷規則，通常要求模型KS值在0.4以上；如下圖所示，為傳統的評價準則

#計算累計分佈函式值
#Kolmogorov-Smirnov檢驗（柯爾莫哥洛夫-斯摩洛夫），用來檢驗資料的分佈是不是符合一個理論的已知分佈
KS_Value <- function(x, y)
{
  gaps_x <- seq(min(x), max(x), length=1000)
  cauculate_x <- numeric()
  for(i in 1:length(gaps_x)){
    cauculate_x[i] <- sum(x<=gaps_x[i])/length(x)
  }
  gaps_x <- sort((gaps_x-min(gaps_x))/(max(gaps_x)-min(gaps_x)))
  gaps_y <- seq(min(y), max(y), length=1000)
  cauculate_y <- numeric()
  for(i in 1:length(gaps_y)){
    cauculate_y[i] <- sum(y<=gaps_y[i])/length(y)
  }
  gaps_y <- sort((gaps_y-min(gaps_y))/(max(gaps_y)-min(gaps_y)))
  return(list(df = data.frame(rbind(data.frame(Gaps = gaps_x,Cauculate = cauculate_x,Type = 'Positive'),data.frame(Gaps = gaps_y,Cauculate = cauculate_y,Type = 'Negtive'))), 
   KS = max(abs(cauculate_y-cauculate_x)), 
   x = gaps_y[which.max(abs(cauculate_y-cauculate_x))],
   y = abs(cauculate_x[which.max(abs(cauculate_y-cauculate_x))]-cauculate_y[which.max(abs(cauculate_y+cauculate_x))])/2))
}
#準備K-S資料
testks <- ifelse(test$churn == "yes", 1, 0)
ks_data <- as.data.frame(cbind(good_bad=testks, prob ))
good_ks <- ks_data[which(ks_data$good_bad==0),'prob']
bad_ks <- ks_data[which(ks_data$good_bad==1),'prob']
#生成K-S值
ks <- KS_Value(bad_ks,good_ks) #KS_Value屬於上面自定義的KS值函式
#繪製K-S曲線
ggplot(data = ks$df, mapping = aes(x = Gaps, y = Cauculate, colour = Type)) + geom_line() + theme(legend.position='none') + annotate(geom = 'text', x = ks$x, y = ks$y, label = paste('K-S Value: ', round(ks$KS,2))) + labs(x = 'Probility', y = 'CDF')  

上圖結果顯示，K-S統計量的值為0.52，根據傳統的評價準則，也說明該模型很好。