如何進行邏輯回歸分析
邏輯回歸是當y=f(x)，而y為分類變量的時候的邏輯曲線擬合的方法。這種模型通常的用法就是通過給定的一個x的預測值來預測y。這些預測值可以說連續的、分類的，或者是混合的。
通常來說，分類變量y有多種不同的假設值。其中，最簡單的一個例子就是y為一個二元變量，這意味着我們可以假設每個y值為1或0。在機器學習里，一個非常常用的案例就是郵件分類：有一個針對每封郵件的給定屬性集，比如字數、鏈接和圖片。這時，我們要用的算法就是根據它判斷這封郵件是垃圾郵件（1）或者不是非垃圾郵件（0）。在這篇博客當中，基於這個預測變量是二元的，我把它稱之為“二元邏輯回歸”。當然，邏輯回歸同樣也可以預測一個超過2個值的獨立變量。在這篇博客中，我們稱它為“多元邏輯回歸”。在這里，一個典型的例子就是分類一部電影，它是“有趣”、“一般般”還是“無聊”。
R運行邏輯回歸
R可以讓邏輯回歸建模過程變得很簡單。我們可以使用glm()函數進行邏輯回歸建模，而且，它的擬合過程和我們之前所做的線性回歸沒有很大的區別。在這篇博客中，我們將介紹如何一步一步的進行邏輯回歸建模。
數據集
我們將以Titanic數據集作為實例進行分析。關於這個數據集，在網上流傳了好幾個免費的版本，而我則建議大家使用在Kaggle上的那個版本，因為它已經被用過很多次了（如果你要下載這個數據集，你需要在Kaggle注冊一個賬號）。
這個（訓練）數據集是關於一些乘客的信息（精確的數值為889），而這個比賽的目標就是基於諸如服務等級、性別、年齡等一些因素來預測一個顧客是否幸存（1表示乘客幸存，0則表示死亡）。在這里，你已經看到了，這個實例即有分類變量，也有連續變量。
數據清洗過程
當我們拿到一個真實的數據集的時候，我們需要確認這個數據集有哪些缺失值或者異常值。因此，我們需要做這方面的准備工作以便我們后面對此進行數據分析。第一步則是我們使用read.csv()來讀取數據集，可以在這里下載數據集。
確認參數na.strings相當於c(“”)，因此，每個缺失值可以用NA表示。這可以幫助我們進行下一步的工作。
-
training.data.raw <- read.csv('train.csv',header=T,na.strings=c(""))
現在，我們要查閱一下缺失值，並觀察一下每個變量存在多少異常值，使用sapply()函數來展示數據框的每一列的值。
-
sapply(training.data.raw,function(x) sum(is.na(x))) -
PassengerId Survived Pclass Name Sex -
0 0 0 0 0 -
Age SibSp Parch Ticket Fare -
177 0 0 0 0 -
Cabin Embarked -
687 2 -
sapply(training.data.raw, function(x) length(unique(x))) -
PassengerId Survived Pclass Name Sex -
891 2 3 891 2 -
Age SibSp Parch Ticket Fare -
89 7 7 681 248 -
Cabin Embarked -
148 4
在這里對缺失值進行可視化會有幫助：Amelia包有一個特別的作圖函數missmap()來針對我們的數據集進行作圖，並且對缺失值做標記。
-
library(Amelia) -
missmap(training.data.raw, main = "Missing values vs observed")
這個變量有非常多的缺失值，我們都不會使用他們。我們也會把Passengerld去掉，因為它只有索引和票。
使用subset()函數，我們可以從原始數據集中只選擇與之相關的列。
-
data <- subset(training.data.raw,select=c(2,3,5,6,7,8,10,12))
觀察一下缺失值
現在，我們需要觀察一下別的缺失值。R可以通過在擬合函數內設定一個參數集來產生擬合一個廣義線性模型。不過，我個人認為，我們應當盡可能的手動去除這些缺失值。這里有很多種方法來說實現，而一種常用的方法就是以均值，或中位數，或者是一個眾數來代替缺失值。這里，我們會使用均值。
data$Age[is.na(data$Age)] <- mean(data$Age,na.rm=T)
當我們把它看成是分類變量的時候，使用read.table()函數或read.csv()函數，而默認情況下是分類變量被轉成因子類型。一個因子就是R處理分類變量的方式。
我們可以使用下面的代碼來查閱它：
-
is.factor(data$Sex) -
TRUE -
-
is.factor(data$Embarked) -
TRUE
為了能更容易的理解R是怎么處理分類變量的，我們可以使用contrasts()函數。這個函數向我們展示了R是如何優化變量的，並把它轉譯成模型。
-
contrasts(data$Sex) -
male -
female 0 -
male 1 -
contrasts(data$Embarked) -
Q S -
C 0 0 -
Q 1 0 -
S 0 1
比如，你可以看見在變量sex，female可被作為參考。在乘船上的缺失值，由於這里只有兩行，我們可以把這兩行刪掉(我們也可以用眾數取代缺失值，然后保留這些數據點）。
-
data <- data[!is.na(data$Embarked),] -
rownames(data) <- NULL
在繼續擬合過程之前，讓我們回顧一下數據清洗和重構的重要性。這樣的預處理通常對於獲得一個好而且預測能力強的模型是非常重要的。
模型擬合
我們把數據分成兩個部分，訓練和測試集。訓練數據集可用於模型的擬合，而測試數據集則對此進行測試。
-
train <- data[1:800,] -
test <- data[801:889,]
現在，讓我們擬合一下這個模型。確認我們在glm()函數使用參數family=binomial。
-
model <- glm(Survived ~.,family=binomial(link='logit'),data=train)
現在，使用summary()函數來獲得這個模型的結果。
-
summary(model) -
Call: -
glm(formula = Survived ~ ., family = binomial(link = "logit"), -
data = train) -
-
Deviance Residuals: -
Min 1Q Median 3Q Max -
-2.6064 -0.5954 -0.4254 0.6220 2.4165 -
-
Coefficients: -
Estimate Std. Error z value Pr(>|z|) -
(Intercept) 5.137627 0.594998 8.635 < 2e-16 *** -
Pclass -1.087156 0.151168 -7.192 6.40e-13 *** -
Sexmale -2.756819 0.212026 -13.002 < 2e-16 *** -
Age -0.037267 0.008195 -4.547 5.43e-06 *** -
SibSp -0.292920 0.114642 -2.555 0.0106 * -
Parch -0.116576 0.128127 -0.910 0.3629 -
Fare 0.001528 0.002353 0.649 0.5160 -
EmbarkedQ -0.002656 0.400882 -0.007 0.9947 -
EmbarkedS -0.318786 0.252960 -1.260 0.2076 -
--- -
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -
-
(Dispersion parameter for binomial family taken to be 1) -
-
Null deviance: 1065.39 on 799 degrees of freedom -
Residual deviance: 709.39 on 791 degrees of freedom -
AIC: 727.39 -
-
Number of Fisher Scoring iterations: 5
解讀邏輯回歸模型的結果
現在，我們可以分析這個擬合的模型，閱讀一下它的結果。首先，我們可以看到SibSp，Fare和Embarked並不是重要的變量。對於那些比較重要的變量，sex的p值最低，這說明sex與乘客的生存幾率的關系是最強的。這個負系數預測變量表明其它變量都一樣，男乘客生存的機率更低。記住，邏輯模型的因變量是對數機率：ln(odds) = ln(p/(1-p)) = ax1 + bx2 + … + z*xn。male是一個優化變量，男性的生還機率下載2.75個對數機率，而年齡的增大則下降0.037個對數機率。
現在，我們可以使用anova()函數來分析這個數據表的偏差了：
-
anova(model, test="Chisq") -
Analysis of Deviance Table -
Model: binomial, link: logit -
Response: Survived -
Terms added sequentially (first to last) -
-
Df Deviance Resid. Df Resid. Dev Pr(>Chi) -
NULL 799 1065.39 -
Pclass 1 83.607 798 981.79 < 2.2e-16 *** -
Sex 1 240.014 797 741.77 < 2.2e-16 *** -
Age 1 17.495 796 724.28 2.881e-05 *** -
SibSp 1 10.842 795 713.43 0.000992 *** -
Parch 1 0.863 794 712.57 0.352873 -
Fare 1 0.994 793 711.58 0.318717 -
Embarked 2 2.187 791 709.39 0.334990 -
--- -
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
空偏差和殘差的不同展示了我們的模型是空模型（只有截距）。這個溝閱讀，效果越好。分析這個圖表，我們可以看到每增加一個變量，同時也會降低其中的偏差。再一次，添加Pclass、Sex和Age能顯著的降低殘差。其它變量對於模型的提高沒什么用，盡管SibSp的p值很低。一個p值較大的值預示着一個沒有變量的模型多少解釋了模型的變化。最終，你想要看到的就是偏差的降低和AIC。
然而，這里並不存在線性回歸的殘差R^2,McFaddenR^2指數可用於模型的擬合。
-
library(pscl) -
pR2(model) -
llh llhNull G2 McFadden r2ML r2CU -
-354.6950111 -532.6961008 356.0021794 0.3341513 0.35917
評估模型的預測能力
在上述的步驟中，我們簡單的評估一下模型的擬合效果。現在，我們看一下使用心得數據集，它的預測結果y是什么樣的。設定參數type=’response’，R可以輸出形如P(Y=1|X）的概率。我們的預測邊界就是0.5。如果P(Y=1|X)>0.5，y=1或y=0。記住，其它的決策邊界可能是更好的選擇。
-
fitted.results <- predict(model,newdata=subset(test,select=c(2,3,4,5,6,7,8)),type='response') -
fitted.results <- ifelse(fitted.results > 0.5,1,0) -
-
misClasificError <- mean(fitted.results != test$Survived) -
print(paste('Accuracy',1-misClasificError)) -
"Accuracy 0.842696629213483"
精度為0.84，說明這個模型擬合效果不錯。然而，你要記住，結果多少都依賴於手動改變數據集的數據。因此，如果你想要得到更精確的精度，你最好執行一些交叉檢驗，例如k次交叉檢驗。
作為最后一步，我們會做ROC曲線並計算AUC（曲線下的面積），它常用於預測二元分類器的模型表現。
ROC曲線是一種曲線，它可以通過設定各種極值來讓正例律（TPR)來抵消反正例律(FPR)，它就在ROC曲線之下。通常來說，一個預測能力強的模型應當能讓ROC接近1（1是理想的）而不是0.5。
-
library(ROCR) -
p <- predict(model, newdata=subset(test,select=c(2,3,4,5,6,7,8)), type="response") -
pr <- prediction(p, test$Survived) -
prf <- performance(pr, measure = "tpr", x.measure = "fpr") -
plot(prf) -
-
auc <- performance(pr, measure = "auc") -
auc <- auc@y.values[[1]] -
auc -
0.8647186
下面是ROC圖：
文章來源：數據人網