隨機森林——二元分類的利器之Kaggle初體驗Titanic: Machine Learning from Disaster

01-24

題圖來自新浪公開課：泰坦尼克號全紀錄_新浪公開課_考試頻道_新浪網

本文主要參考以下文章：

Exploring Survival on the Titanic

【Kaggle實例分析】Titanic Machine Learning from Disaster

機器學習（二）如何做到Kaggle排名前2%

從2月份開始，學習R語言和數據分析已經一段時間了，前面也通讀了《R語言實戰》並且盡量用自己的語言寫了十來篇學習筆記，並且緊跟大數據分析社群的學習進度陸續了學習了以下內容：

第一講：零基礎入門方法論
第二講：數據結構入門
第三講：簡單數據處理和分析
第四講：複雜數據處理和分析
第五講：SQL從入門到精通

本文是第五講的課後作業之一，旨在梳理一遍之前所學知識，特別是《R語言實戰》當中的系統性的知識，以Kaggle網站入門問題：泰坦尼克生存率預測為基礎，完整地走一遍從數據預處理、特徵工程、建模、預測、驗證的整個過程。

1 載入數據

在載入數據之前，先安裝和導入需要用到的相關包。

#管理安裝包#數據輸入，輸出packages <- c("readr")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }#圖形包packages <- c("ggplot2")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }#圖形包packages <- c("ggthemes")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("scales")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("plyr")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("stringr")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("InformationValue")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("MLmetrics")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }#決策樹模型包packages <- c("rpart")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("randomForest")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("dplyr")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("e1071")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("Amelia")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("party")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("gbm")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }packages <- c("class")if (length(setdiff(packages, rownames(installed.packages()))) > 0) { install.packages(setdiff(packages, rownames(installed.packages()))) }library(readr) # File read / writelibrary(ggplot2) # Data visualizationlibrary(ggthemes) # Data visualizationlibrary(scales) # Data visualizationlibrary(plyr)library(stringr) # String manipulationlibrary(InformationValue) # IV / WOE calculationlibrary(MLmetrics) # Mache learning metrics.e.g. Recall, Precision, Accuracy, AUClibrary(rpart) # Decision tree utilslibrary(randomForest) # Random Forestlibrary(dplyr) # Data manipulationlibrary(e1071) # SVMlibrary(Amelia) # Missing value utilslibrary(party) # Conditional inference treeslibrary(gbm) # AdaBoostlibrary(class) # KNN

導入相關數據：

#數據導入> train <- read.csv("train.csv", stringsAsFactors = F)> test <- read.csv("test.csv", stringsAsFactors = F)

初步觀察數據：

> head(train) PassengerId Survived Pclass1 1 0 32 2 1 13 3 1 34 4 1 15 5 0 36 6 0 3 Name Sex Age SibSp Parch1 Braund, Mr. Owen Harris male 22 1 02 Cumings, Mrs. John Bradley (Florence Briggs Thayer) female 38 1 03 Heikkinen, Miss. Laina female 26 0 04 Futrelle, Mrs. Jacques Heath (Lily May Peel) female 35 1 05 Allen, Mr. William Henry male 35 0 06 Moran, Mr. James male NA 0 0 Ticket Fare Cabin Embarked1 A/5 21171 7.2500 S2 PC 17599 71.2833 C85 C3 STON/O2. 3101282 7.9250 S4 113803 53.1000 C123 S5 373450 8.0500 S6 330877 8.4583 Q> head(test) PassengerId Pclass Name Sex Age1 892 3 Kelly, Mr. James male 34.52 893 3 Wilkes, Mrs. James (Ellen Needs) female 47.03 894 2 Myles, Mr. Thomas Francis male 62.04 895 3 Wirz, Mr. Albert male 27.05 896 3 Hirvonen, Mrs. Alexander (Helga E Lindqvist) female 22.06 897 3 Svensson, Mr. Johan Cervin male 14.0 SibSp Parch Ticket Fare Cabin Embarked1 0 0 330911 7.8292 Q2 1 0 363272 7.0000 S3 0 0 240276 9.6875 Q4 0 0 315154 8.6625 S5 1 1 3101298 12.2875 S6 0 0 7538 9.2250 S>

不難發現，兩個數據集除了Suvived欄位不同以外，其他欄位均相同，由於後續要對訓練數據和測試數據做相同的轉換，為了避免重複操作和出現不一致的情況，更為了避免可能碰到的Categorical類型新level的問題，這裡建議將兩個數據集合同，統一操作。

合併數據：

#數據合併> data <- bind_rows(train,test)> train.row <- 1:nrow(train)> test.row <- (1+nrow(train)):(nrow(train)+nrow(test))

2 數據預覽

觀察合併後的數據：

> str(data)"data.frame": 1309 obs. of 12 variables: $ PassengerId: int 1 2 3 4 5 6 7 8 9 10 ... $ Survived : int 0 1 1 1 0 0 0 0 1 1 ... $ Pclass : int 3 1 3 1 3 3 1 3 3 2 ... $ Name : chr "Braund, Mr. Owen Harris" "Cumings, Mrs. John Bradley (Florence Briggs Thayer)" "Heikkinen, Miss. Laina" "Futrelle, Mrs. Jacques Heath (Lily May Peel)" ... $ Sex : chr "male" "female" "female" "female" ... $ Age : num 22 38 26 35 35 NA 54 2 27 14 ... $ SibSp : int 1 1 0 1 0 0 0 3 0 1 ... $ Parch : int 0 0 0 0 0 0 0 1 2 0 ... $ Ticket : chr "A/5 21171" "PC 17599" "STON/O2. 3101282" "113803" ... $ Fare : num 7.25 71.28 7.92 53.1 8.05 ... $ Cabin : chr "" "C85" "" "C123" ... $ Embarked : chr "S" "C" "S" "S" ...> summary(data) PassengerId Survived Pclass Name Min. : 1 Min. :0.0000 Min. :1.000 Length:1309 1st Qu.: 328 1st Qu.:0.0000 1st Qu.:2.000 Class :character Median : 655 Median :0.0000 Median :3.000 Mode :character Mean : 655 Mean :0.3838 Mean :2.295 3rd Qu.: 982 3rd Qu.:1.0000 3rd Qu.:3.000 Max. :1309 Max. :1.0000 Max. :3.000 NA"s :418 Sex Age SibSp Parch Length:1309 Min. : 0.17 Min. :0.0000 Min. :0.000 Class :character 1st Qu.:21.00 1st Qu.:0.0000 1st Qu.:0.000 Mode :character Median :28.00 Median :0.0000 Median :0.000 Mean :29.88 Mean :0.4989 Mean :0.385 3rd Qu.:39.00 3rd Qu.:1.0000 3rd Qu.:0.000 Max. :80.00 Max. :8.0000 Max. :9.000 NA"s :263 Ticket Fare Cabin Embarked Length:1309 Min. : 0.000 Length:1309 Length:1309 Class :character 1st Qu.: 7.896 Class :character Class :character Mode :character Median : 14.454 Mode :character Mode :character Mean : 33.295 3rd Qu.: 31.275 Max. :512.329 NA"s :1

從上可知，合併後的數據集包含12個變數，1309條數據，其中891條為訓練數據，418條為測試數據。

其中生存情況（Survived）中缺失值NA有418個（需要預測的），年齡（Age）中缺失值有263個，船票費用（Fare）中缺失值有1個。變數解釋如下：

PassengerId 整型變數，標識乘客的ID，遞增變數，對預測無幫助
Survived 整型變數，標識該乘客是否倖存。0表示遇難，1表示倖存。將其轉換為factor變數比較方便處理
Pclass 整型變數，標識乘客的社會-經濟狀態，1代表Upper，2代表Middle，3代表Lower
Name 字元型變數，除包含姓和名以外，還包含Mr. Mrs. Dr.這樣的具有西方文化特點的信息
Sex 字元型變數，標識乘客性別，適合轉換為factor類型變數
Age 整型變數，標識乘客年齡，有缺失值
SibSp 整型變數，代表兄弟姐妹及配偶的個數。其中Sib代表Sibling也即兄弟姐妹，Sp代表Spouse也即配偶
Parch 整型變數，代表父母或子女的個數。其中Par代表Parent也即父母，Ch代表Child也即子女
Ticket 字元型變數，代表乘客的船票號
Fare 數值型，代表乘客的船票價
Cabin 字元型，代表乘客所在的艙位，有缺失值
Embarked 字元型，代表乘客登船口岸，適合轉換為factor型變數

我們需要根據這些數據對生存情況（Survived）——因變數進行預測，可供使用的自變數一共有11個。

3 探索式數據分析

3.1 乘客社會等級越高，辛存率越高

將Survived變數轉換為factor類型變數。

> data$Survived <- factor(data$Survived)

通過以下代碼統計並繪製不同Pclass的乘客倖存和遇難的人數。

> ggplot(data = data[1:nrow(train),], mapping = aes(x = Pclass, y = ..count.., fill=Survived)) + geom_bar(stat = "count", position="dodge") + xlab("Pclass") + ylab("Count") + ggtitle("How Pclass impact survivor") + scale_fill_manual(values=c("blue", "yellow")) + geom_text(stat = "count", aes(label = ..count..), position=position_dodge(width=1), , vjust=-0.5) + theme(plot.title = element_text(hjust = 0.5), legend.position="bottom")>

從上圖可見，隨著Pclass的增加，即社會地位的降低，倖存率不斷下降，最低等級的乘客只有不到25%倖存。

可以通過計算Pclass的WOE和IV值來更為定量地計算出該變數的預測價值，下面為相關代碼：

>WOETable(X=factor(data$Pclass[1:nrow(train)]),Y=data$Survived[1:nrow(train)]) CAT GOODS BADS TOTAL PCT_G PCT_B WOE IV1 1 136 80 216 0.3976608 0.1457195 1.0039160 0.252927922 2 87 97 184 0.2543860 0.1766849 0.3644848 0.028320873 3 119 372 491 0.3479532 0.6775956 -0.6664827 0.21970095> IV(X=factor(data$Pclass[1:nrow(train)]), Y=data$Survived[1:nrow(train)])[1] 0.5009497attr(,"howgood")[1] "Highly Predictive"

從結果可以看出，Pclass的IV為0.5，且「Highly Predictive」。因此有充分的理由暫將Pclass作為預測模型的特徵變數之一。

3.2 不同頭銜的乘客倖存率不同

由於乘客姓名的重複度太低，不適合直接使用。然後我們注意到在乘客名字（Name）中，有一個非常顯著的特點：包含Mr. Mrs. Dr.等具有文化特徵的信息，可將之抽取出來並作為一個獨立的新變數——Title。

以下代碼從姓名中抽取乘客的Title並進行factor類型轉換。

> data$Title <- sapply(data$Name, FUN=function(x) {strsplit(x, split="[,.]")[[1]][2]})> data$Title <- sub(" ", "", data$Title)> data$Title[data$Title %in% c("Mme", "Mlle")] <- "Mlle"> data$Title[data$Title %in% c("Capt", "Don", "Major", "Sir")] <- "Sir"> data$Title[data$Title %in% c("Dona", "Lady", "the Countess", "Jonkheer")] <- "Lady"> data$Title <- factor(data$Title)>

通過以下代碼統計並繪製不同Title的乘客倖存和遇難的人數。

> ggplot(data = data[1:nrow(train),], mapping = aes(x = Title, y = ..count.., fill=Survived)) + geom_bar(stat = "count", position="stack") + xlab("Title") + ylab("Count") + ggtitle("How Title impact survivor") + scale_fill_discrete(name="Survived", breaks=c(0, 1), labels=c("Perish", "Survived")) + geom_text(stat = "count", aes(label = ..count..), position=position_stack(vjust = 0.5)) + theme(plot.title = element_text(hjust = 0.5), legend.position="bottom")

從上圖可看出，Title為Mr的乘客倖存比例非常小，而Title為Mrs和Miss的乘客倖存比例非常大。

為了定量地評估Title的預測價值，這裡再次使用WOE和IV來進行定量計算，以下為相關代碼：

> WOETable(X=data$Title[1:nrow(train)], Y=data$Survived[1:nrow(train)]) CAT GOODS BADS TOTAL PCT_G PCT_B WOE IV1 Col 1 1 2 0.002873563 0.001808318 0.46315552 4.933741e-042 Dr 3 4 7 0.008620690 0.007233273 0.17547345 2.434548e-043 Lady 2 1 3 0.005747126 0.001808318 1.15630270 4.554455e-034 Master 23 17 40 0.066091954 0.030741410 0.76543639 2.705859e-025 Miss 127 55 182 0.364942529 0.099457505 1.30000942 3.451330e-016 Mlle 3 3 3 0.008620690 0.005424955 0.46315552 1.480122e-037 Mr 81 436 517 0.232758621 0.788426763 -1.22003757 6.779360e-018 Mrs 99 26 125 0.284482759 0.047016275 1.80017883 4.274821e-019 Ms 1 1 1 0.002873563 0.001808318 0.46315552 4.933741e-0410 Rev 6 6 6 0.017241379 0.010849910 0.46315552 2.960244e-0311 Sir 2 3 5 0.005747126 0.005424955 0.05769041 1.858622e-05> > IV(X=data$Title[1:nrow(train)], Y=data$Survived[1:nrow(train)])[1] 1.487853attr(,"howgood")[1] "Highly Predictive">

從計算結果可見，IV為1.520702，且」Highly Predictive」。因此，有理由暫將Title作為預測模型中的一個特徵變數。

3.3 女性倖存率遠高於男性

由Titanic號沉沒的背景可知，逃生時遵循「婦女與小孩先走」的規則，由此猜想，Sex變數應該對預測乘客倖存有所幫助。

下列數據驗證了這一猜想：大部分女性（233/(233+81)=74.20%）得以倖存，而男性中只有很小部分（109/(109+468)=22.85%）倖存。

#數據因子化> data$Sex <- as.factor(data$Sex)#數據可視化> ggplot(data = data[1:nrow(train),], mapping = aes(x = Sex, y = ..count.., fill=Survived)) + geom_bar(stat = "count", position="dodge") + xlab("Sex") + ylab("Count") + ggtitle("How Sex impact survivo") + geom_text(stat = "count", aes(label = ..count..), position=position_dodge(width=1), , vjust=-0.5) + theme(plot.title = element_text(hjust = 0.5), legend.position="bottom")>

WOE和IV分析可知，Sex變數的IV為1.34並且指示」Highly Predictive」，因此也可以暫時將Sex作為後續模型的特徵變數之一。

> WOETable(X=data$Sex[1:nrow(train)], Y=data$Survived[1:nrow(train)]) CAT GOODS BADS TOTAL PCT_G PCT_B WOE IV1 female 233 81 314 0.6812865 0.147541 1.5298770 0.81656512 male 109 468 577 0.3187135 0.852459 -0.9838327 0.5251163> > IV(X=data$Sex[1:nrow(train)], Y=data$Survived[1:nrow(train)])[1] 1.341681attr(,"howgood")[1] "Highly Predictive">

3.4 未成年人倖存率高於成年人

按照「婦女與小孩先走」的規則，未成年人應該有更大可能倖存。如下圖所示，Age < 18的乘客中，倖存人數確實高於遇難人數。同時青壯年乘客中，遇難人數遠高於倖存人數。因此也可以暫時將Age作為後續模型的特徵變數之一。

> ggplot(data = data[(!is.na(data$Age)) & row(as.matrix(data[, "Age"])) <= 891, ], aes(x = Age, color=Survived)) + geom_line(aes(label=..count..), stat = "bin", binwidth=5) + labs(title = "How Age impact survivor", x = "Age", y = "Count", fill = "Survived")Warning: Ignoring unknown aesthetics: label

3.5 配偶及兄弟姐妹數十種的乘客更易倖存

對於SibSp變數，分別統計繪製出倖存與遇難人數。

> ggplot(data = data[1:nrow(train),], mapping = aes(x = SibSp, y = ..count.., fill=Survived)) + geom_bar(stat = "count", position="dodge") + labs(title = "How SibSp impact survivor", x = "Sibsp", y = "Count", fill = "Survived") + geom_text(stat = "count", aes(label = ..count..), position=position_dodge(width=1), , vjust=-0.5) + theme(plot.title = element_text(hjust = 0.5), legend.position="bottom")>

從上圖可見，SibSp為0的乘客，倖存率低於1/3；SibSp為1或2的乘客，倖存率高於50%；SibSp大於等於3的乘客，倖存率非常低。

同樣可以通過計算WOE和IV定量計算SibSp的預測價值，其中IV為0.1448994，且」Highly Predictive」。因此也可以暫時將SibSp作為後續模型的特徵變數之一。

> WOETable(X=as.factor(data$SibSp[1:nrow(train)]), Y=data$Survived[1:nrow(train)]) CAT GOODS BADS TOTAL PCT_G PCT_B WOE IV1 0 210 398 608 0.593220339 0.724954463 -0.2005429 0.0264183492 1 112 97 209 0.316384181 0.176684882 0.5825894 0.0813873343 2 13 15 28 0.036723164 0.027322404 0.2957007 0.0027798114 3 4 12 16 0.011299435 0.021857923 -0.6598108 0.0069666045 4 3 15 18 0.008474576 0.027322404 -1.1706364 0.0220639536 5 5 5 5 0.014124294 0.009107468 0.4388015 0.0022013917 8 7 7 7 0.019774011 0.012750455 0.4388015 0.003081947> > IV(X=as.factor(data$SibSp[1:nrow(train)]), Y=data$Survived[1:nrow(train)])[1] 0.1448994attr(,"howgood")[1] "Highly Predictive">

3.6 父母與子女數為1到3的乘客更可能倖存

對於Parch變數，分別統計繪製出倖存與遇難人數。

> ggplot(data = data[1:nrow(train),], mapping = aes(x = Parch, y = ..count.., fill=Survived)) + geom_bar(stat = "count", position="dodge") + labs(title = "How Parch impact survivor", x = "Parch", y = "Count", fill = "Survived") + geom_text(stat = "count", aes(label = ..count..), position=position_dodge(width=1), , vjust=-0.5) + theme(plot.title = element_text(hjust = 0.5), legend.position="bottom")>

從上圖可知，Parch為0的乘客，倖存率低於1/3；Parch為1到3的乘客，倖存率高於50%；Parch大於等於4的乘客，幾乎沒有倖存的。

同樣可以通過計算WOE和IV定量計算Parch的預測價值，其中IV為0.1166611，且」Highly Predictive」。因此也可以暫時將Parch作為後續模型的特徵變數之一。

> WOETable(X=as.factor(data$Parch[1:nrow(train)]), Y=data$Survived[1:nrow(train)]) CAT GOODS BADS TOTAL PCT_G PCT_B WOE IV1 0 233 445 678 0.671469741 0.810564663 -0.1882622 0.0261863122 1 65 53 118 0.187319885 0.096539162 0.6628690 0.0601757283 2 40 40 80 0.115273775 0.072859745 0.4587737 0.0194584404 3 3 2 5 0.008645533 0.003642987 0.8642388 0.0043233945 4 4 4 4 0.011527378 0.007285974 0.4587737 0.0019458446 5 1 4 5 0.002881844 0.007285974 -0.9275207 0.0040849227 6 1 1 1 0.002881844 0.001821494 0.4587737 0.000486461> IV(X=as.factor(data$Parch[1:nrow(train)]), Y=data$Survived[1:nrow(train)])[1] 0.1166611attr(,"howgood")[1] "Highly Predictive">

3.7 FamilySize為2到4的乘客倖存率較高

SibSp和Parch兩個變數都說明，當乘客沒有親人時倖存率較低，當乘客有少數親人時，倖存率高於50%，而當乘客親人數過高時，倖存率反而降低。考慮到這兩個變數有著相似的特性，我們可以將SibSp和Parch相加，生成一個新的變數，FamilySize，也算是一種降維的手段。

> data$FamilySize <- data$SibSp + data$Parch + 1> ggplot(data = data[1:nrow(train),], mapping = aes(x = FamilySize, y = ..count.., fill=Survived)) + geom_bar(stat = "count", position="dodge") + xlab("FamilySize") + ylab("Count") + ggtitle("How FamilySize impact survivor") + geom_text(stat = "count", aes(label = ..count..), position=position_dodge(width=1), , vjust=-0.5) + theme(plot.title = element_text(hjust = 0.5), legend.position="bottom")>

從上圖可知，FamilySize為1的乘客，倖存率約為30%；FamilySize為2到4的乘客，倖存率高於50%；FamilySize大於等於4的乘客，辛存率大幅下降。

根據下列代碼，通過WOE和IV分析評估FamilySize變數的預測價值。

> WOETable(X=as.factor(data$FamilySize[1:nrow(train)]), Y=data$Survived[1:nrow(train)]) CAT GOODS BADS TOTAL PCT_G PCT_B WOE IV1 1 163 374 537 0.459154930 0.68123862 -0.3945249 0.08761755392 2 89 72 161 0.250704225 0.13114754 0.6479509 0.07746686163 3 59 43 102 0.166197183 0.07832423 0.7523180 0.06610840574 4 21 8 29 0.059154930 0.01457195 1.4010615 0.06246349985 5 3 12 15 0.008450704 0.02185792 -0.9503137 0.01274106436 6 3 19 22 0.008450704 0.03460838 -1.4098460 0.03687829407 7 4 8 12 0.011267606 0.01457195 -0.2571665 0.00084976658 8 6 6 6 0.016901408 0.01092896 0.4359807 0.00260387129 11 7 7 7 0.019718310 0.01275046 0.4359807 0.0030378497> IV(X=as.factor(data$FamilySize[1:nrow(train)]), Y=data$Survived[1:nrow(train)])[1] 0.3497672attr(,"howgood")[1] "Highly Predictive">

計算可得，FamilySize變數的IV為0.3497672，高於SibSp與Parch的IV，且「Highly Predictive」。因此，有必要可以用這個派生變數FamilySize代替原SibSp和Parch變數作為後續預測模型的特徵變數。

3.8 共票號乘客倖存率較高

Ticket變數重複度非常低，直接利用的價值不大。我們猜想票號相同的乘客，有可能是一家人，要麼同時倖存要麼同時遇難。首先統計出每張票對應的乘客數。

> ticket.count <- aggregate(data$Ticket, by=list(data$Ticket),function(x) sum(!is.na(x)))> summary(ticket.count) Group.1 x Length:929 Min. : 1.000 Class :character 1st Qu.: 1.000 Mode :character Median : 1.000 Mean : 1.409 3rd Qu.: 1.000 Max. :11.000

總共有929個不同的票號，最多有11個人共享票號，並且單獨票號的乘客占多數。

接下來將所有乘客按照Ticket分為兩組，一組是使用單獨票號，另一組是與他人共享票號，並統計繪製出各組的倖存與遇難人數。

#將統計好的同票號乘客數賦值給各個乘客> data$TicketCount <- apply(data,1,function(x) ticket.count[which(ticket.count[, 1] == x["Ticket"]), 2])#將TicketCount變數因子化為Share和Unique> data$TicketCount <- factor(sapply(data$TicketCount, function(x) ifelse(x > 1, "Share", "Unique")))#繪製圖形> ggplot(data = data[1:nrow(train),], mapping = aes(x = TicketCount, y = ..count.., fill=Survived)) + geom_bar(stat = "count", position="dodge") + xlab("TicketCount") + ylab("Count") + ggtitle("How TicketCount impact survivor") + geom_text(stat = "count", aes(label = ..count..), position=position_dodge(width=1), , vjust=-0.5) + theme(plot.title = element_text(hjust = 0.5), legend.position="bottom")

由上圖可見，單獨的乘客，只有130/(130+351)=27%倖存，而與他人同票號的乘客有212/(212+198)=51.7%倖存。

計算TicketCount變數的WOE與IV如下。

> WOETable(X=data$TicketCount[1:nrow(train)], Y=data$Survived[1:nrow(train)]) CAT GOODS BADS TOTAL PCT_G PCT_B WOE IV1 Share 212 198 410 0.619883 0.3606557 0.5416069 0.14039932 Unique 130 351 481 0.380117 0.6393443 -0.5199641 0.1347889> IV(X=data$TicketCount[1:nrow(train)], Y=data$Survived[1:nrow(train)])[1] 0.2751882attr(,"howgood")[1] "Highly Predictive">

其IV為0.2751882，且」Highly Predictive」。因此，也可以考慮將TicketCount作為後續模型的一個特徵變數。

3.8 支出船票費用越高倖存率越高

Fare變數是數值型變數，通常會猜想支出費用是否與倖存率之前有某種關係。

> ggplot(data = data[(!is.na(data$Fare)) & row(as.matrix(data[, "Fare"])) <= 891, ], aes(x = Fare, color=Survived)) + geom_line(aes(label=..count..), stat = "bin", binwidth=10) + labs(title = "How Fare impact survivor", x = "Fare", y = "Count", fill = "Survived")Warning: Ignoring unknown aesthetics: label

由上圖不難發現，Fare越大，倖存率越高。

3.9 不同倉位的乘客倖存率不同

對於Cabin變數，其值以字母開始，後面跟著數字。按照常理，我們可以大膽的猜想，字母代表某一個區域，數據則表示該區域內的序號，與火車票即有車廂號又有座位號類似。在此將Cabin的首字母提取出來，並分別統計繪製不同首字母倉位對應的乘客的倖存人數。

> ggplot(data[1:nrow(train), ], mapping = aes(x = as.factor(sapply(data$Cabin[1:nrow(train)], function(x) str_sub(x, start = 1, end = 1))), y = ..count.., fill = Survived)) + geom_bar(stat = "count", position="dodge") + xlab("Cabin") + ylab("Count") + ggtitle("How Cabin impact survivor") + geom_text(stat = "count", aes(label = ..count..), position=position_dodge(width=1), , vjust=-0.5) + theme(plot.title = element_text(hjust = 0.5), legend.position="bottom")

tips：將空值轉換為NA的方法：x[x==""]=NA

由上圖可知，B、C、D、E、F倉位的乘客倖存率高於50%，A、F倉位的乘客倖存率也接近50%，另外，還有大量的倉位數據缺失，在後續的建模之前可能有必要對缺失數據進行處理。Cabin變數的WOE和IV計算如下：

> data$Cabin <- sapply(data$Cabin, function(x) str_sub(x, start = 1, end = 1))> WOETable(X=as.factor(data$Cabin[1:nrow(train)]), Y=data$Survived[1:nrow(train)]) CAT GOODS BADS TOTAL PCT_G PCT_B WOE IV1 206 481 687 0.600583090 0.876138434 -0.3776231 0.1040560652 A 7 8 15 0.020408163 0.014571949 0.3368366 0.0019658513 B 35 12 47 0.102040816 0.021857923 1.5408094 0.1235465554 C 35 24 59 0.102040816 0.043715847 0.8476622 0.0494398735 D 25 8 33 0.072886297 0.014571949 1.6098023 0.0938745716 E 24 8 32 0.069970845 0.014571949 1.5689803 0.0869197767 F 8 5 13 0.023323615 0.009107468 0.9403716 0.0133684618 G 2 2 4 0.005830904 0.003642987 0.4703680 0.0010291269 T 1 1 1 0.002915452 0.001821494 0.4703680 0.000514563> > IV(X=as.factor(data$Cabin[1:nrow(train)]), Y=data$Survived[1:nrow(train)])[1] 0.4747148attr(,"howgood")[1] "Highly Predictive">

Cabin的IV為0.1866526，且「Highly Predictive」，因此也可以暫時作為模型的潛在特徵變數。

3.10 不同登船碼頭的乘客倖存率不同

Embarked變數表示登船碼頭，下列代碼統計並繪製不同登船地點的乘客倖存人數。

> ggplot(data[1:nrow(train), ], mapping = aes(x = Embarked, y = ..count.., fill = Survived)) + geom_bar(stat = "count", position="dodge") + xlab("Embarked") + ylab("Count") + ggtitle("How Embarked impact survivor") + geom_text(stat = "count", aes(label = ..count..), position=position_dodge(width=1), , vjust=-0.5) + theme(plot.title = element_text(hjust = 0.5), legend.position="bottom")>

從上圖可見，Embarked變數為S的乘客倖存率僅為217/(217+427)=33.7%，而Embarked為C或為NA的乘客倖存率均高於50%。初步判斷Embarked變數可用於預測乘客是否倖存。Embarked變數的WOE和IV計算如下：

> WOETable(X=as.factor(data$Embarked[1:nrow(train)]), Y=data$Survived[1:nrow(train)]) CAT GOODS BADS TOTAL PCT_G PCT_B WOE IV1 2 2 2 0.005847953 0.003629764 0.47692407 1.057908e-032 C 93 75 168 0.271929825 0.136116152 0.69203545 9.398788e-023 Q 30 47 77 0.087719298 0.085299456 0.02797385 6.769232e-054 S 217 427 644 0.634502924 0.774954628 -0.19996259 2.808509e-02> IV(X=as.factor(data$Embarked[1:nrow(train)]), Y=data$Survived[1:nrow(train)])[1] 0.1231986attr(,"howgood")[1] "Highly Predictive">

計算結果為IV為0.1227284，且「Highly Predictive」，驗證初步判斷。

4 缺失值處理

4.1 列出所有缺失數據

> attach(data)> missing <- list(Pclass=nrow(data[is.na(Pclass), ]))> missing$Name <- nrow(data[is.na(Name), ])> missing$Sex <- nrow(data[is.na(Sex), ])> missing$Age <- nrow(data[is.na(Age), ])> missing$SibSp <- nrow(data[is.na(SibSp), ])> missing$Parch <- nrow(data[is.na(Parch), ])> missing$Ticket <- nrow(data[is.na(Ticket), ])> missing$Fare <- nrow(data[is.na(Fare), ])> missing$Cabin <- nrow(data[is.na(Cabin), ])> missing$Embarked <- nrow(data[is.na(Embarked), ])> for (name in names(missing)) {+ if (missing[[name]][1] > 0) {+ print(paste("", name, " miss ", missing[[name]][1], " values", sep = ""))+ }+ }[1] "Age miss 263 values"[1] "Fare miss 1 values"[1] "Cabin miss 1014 values"[1] "Embarked miss 2 values"> detach(data)

對於缺失數據，通常可以分為三類：

缺失數據最簡單粗暴的辦法就是將包含缺失數據的觀測整個刪除，另外就是用中位數或者平均值進行替換，還有就是根據已有變數進行建模預測。

4.2 預測乘客年齡

年齡信息缺失的乘客數為263，缺失量比較大，不太適合使用中位數或者平均值進行填補，這裡通過使用其他變數進行預測缺失的年齡信息。模型方便使用rpart包中的決策樹模型，決策樹是一種有監督的機器學習方法，既可以用於分類，也可用於回歸。使用格式如下：

rpart(formula, data, weights, subset, na.action = na.rpart,method, model = FALSE, x = FALSE, y = TRUE, parms,control, cost, ...)

其中：

formula 回歸方程lm( )形式：y~x1+x2+x3+x4
data 包含前面公式的數據框na.action 缺失數據的處理辦法。默認辦法是刪除因變數缺失的觀測值，保留自變數缺失的觀測值
method 根據樹末端因變數的數據類型選擇分割方法：anova(連續型)、poisson（計數型）、class（離散型）、exp（生存型）
parms 設置3個參數：先驗概率、損失矩陣、分類純度
control 控制每個節點上的最小樣本量、交叉驗證的次數、複雜性參數

Age變數預測代碼如下：

> age.model <- rpart(Age ~ Pclass + Sex + SibSp + Parch + Fare + Embarked + Title + FamilySize, data=data[!is.na(data$Age), ], method="anova")> data$Age[is.na(data$Age)] <- predict(age.model, data[is.na(data$Age), ])>

4.3 登船碼頭處理

由於Embarked變數僅有2個缺失值，建議採用中位數填補缺失的Embarked值。

從如下數據可見，缺失Embarked信息的乘客的Pclass均為1，且Fare均為80。

> data[is.na(data$Embarked), c("PassengerId", "Pclass", "Fare", "Embarked")] PassengerId Pclass Fare Embarked62 62 1 80 <NA>830 830 1 80 <NA>>

另外，由下圖所見，Embarked為C且Pclass為1的乘客的Fare中位數為80。

> ggplot(data[!is.na(data$Embarked),], aes(x=Embarked, y=Fare, fill=factor(Pclass))) + geom_boxplot() + geom_hline(aes(yintercept=80), color="red", linetype="dashed", lwd=2) + scale_y_continuous(labels=dollar_format()) + theme_few()Warning message:Removed 1 rows containing non-finite values (stat_boxplot).

綜上，可以將缺失的Embarked值設置為『C』，並將其因子化。

> data$Embarked[is.na(data$Embarked)] <- "C"> data$Embarked <- as.factor(data$Embarked)

4.4 船票費用處理

由於缺失Fare值的記錄非常少，一般可直接使用平均值或者中位數填補該缺失值。這裡使用乘客的Fare中位數填補缺失值。代碼如下：

> data$Fare[is.na(data$Fare)] <- median(data$Fare, na.rm=TRUE)>

4.5 倉位數據處理

Cabin變數缺失的記錄數很多，無論是中位數或者平均數都不是特別合適，從其他變數來預測也不是特別合理，考慮到在上一節中，將NA單獨對待時，其IV已經比較高。因此這裡直接將缺失的Cabin設置為一個默認值。

#將Cabin缺失值設置為默認值X> data$Cabin <- as.factor(sapply(data$Cabin, function(x) ifelse(is.na(x), "X", str_sub(x, start = 1, end = 1))))>

5 模型訓練

這裡採用party包的cforest()函數，cforest()函數採用隨機森林模型。

#設置隨機數為1234> set.seed(1234)#建立模型，注意所採用的變數> model <- cforest(Survived ~ Pclass + Title + Sex + Age + SibSp + Parch + FamilySize + TicketCount + Fare + Cabin + Embarked, data = data[train.row, ], controls=cforest_unbiased(ntree=2000, mtry=3))>

6 交叉驗證

通常，應該將訓練數據分為兩部分，一部分用於訓練，另一部分用於驗證。或者使用k-fold交叉驗證。在此將所有訓練數據都用於訓練，然後隨機選取35%數據集用於驗證。

#創建交叉驗證信息顯示函數> cv.summarize <- function(data.true, data.predict) {+ print(paste("Recall:", Recall(data.true, data.predict)))+ print(paste("Precision:", Precision(data.true, data.predict)))+ print(paste("Accuracy:", Accuracy(data.predict, data.true)))+ print(paste("AUC:", AUC(data.predict, data.true)))+ }> set.seed(1234)> cv.test.sample <- sample(1:nrow(train), as.integer(0.35 * nrow(train)), replace = TRUE)> cv.test <- data[cv.test.sample,]> cv.prediction <- predict(model, cv.test, OOB=TRUE, type = "response")> cv.summarize(cv.test$Survived, cv.prediction)[1] "Recall: 0.953846153846154"[1] "Precision: 0.849315068493151"[1] "Accuracy: 0.864951768488746"[1] "AUC: 0.834681697612732">

7 預測

#預測結果> predict.result <- predict(model,data[(1+nrow(train)):(nrow(data)),],OOB=TRUE,type="response")#將需要預測的數據從預測結果中分離處理> output <- data.frame(PassengerId = test$PassengerId, Survived = predict.result)#將結果寫入文件> write.csv(output, file = "cit1.csv", row.names = FALSE)>

然後將結果上傳到Kaggle，該模型在Kaggle上的得分為0.80383，排名986名，前986/7070=13.9%。

8 調優

8.1 去掉關聯特徵

由於FamilySize結合了SibSp和Parch的信息，因此可以嘗試將SiSp和Parch從特徵變數中刪除。

> set.seed(1234)> model2 <- cforest(Survived ~ Pclass + Title + Sex + Age + FamilySize + TicketCount + Fare + Cabin + Embarked, data = data[train.row, ], controls=cforest_unbiased(ntree=2000, mtry=3))> predict.result <- predict(model2, data[test.row, ], OOB=TRUE, type = "response")> submit <- data.frame(PassengerId = test$PassengerId, Survived = predict.result)> write.csv(submit, file = "cit2.csv", row.names = FALSE)>

該模型預測結果在Kaggle的得分仍為0.80383。

8.2 去掉IV較低的變數

從前面的IV分析當中可知，Cabin的IV值相對較低，因此可以考慮將其從模型當中刪除。

> set.seed(1234)> model3 <- cforest(Survived ~ Pclass + Title + Sex + Age + FamilySize + TicketCount + Fare + Embarked, data = data[train.row, ], controls=cforest_unbiased(ntree=2000, mtry=3))> predict.result <- predict(model3, data[test.row, ], OOB=TRUE, type = "response")> submit <- data.frame(PassengerId = test$PassengerId, Survived = predict.result)> write.csv(submit, file = "cit3.csv", row.names = FALSE)>

該模型預測結果在Kaggle的得分仍為0.80383。

8.3 特徵挖掘

對於Name變數，前面已經派生出了Title變數，基於下列原因，可以推測乘客的姓氏可能具有一定的預測價值。

部分西方國家中人名的重複度較高，而姓氏重複度較低，姓氏具有一定辨識度

部分國家的姓氏具有一定的身份識別作用
姓氏相同的乘客，可能是一家人（這一點也基於西方國家姓氏重複度較低這一特點），而一家人同時倖存或遇難的可能性較高

另外，考慮到只出現一次的姓氏不可能同時出現在訓練集和測試集中，不具辨識度和預測作用，因此將只出現一次的姓氏均命名為』Small』。

> data$FamilyID <- paste(as.character(data$FamiliSize),data$Surname,sep="")> data$FamilyID[data$FamilySize <= 2] <- "Small"> famIDs <- data.frame(table(data$FamilyID))> famIDs <- famIDs[famIDs$Freq <= 2,]> data$FamilyID[data$FamilyID %in% famIDs$Var1] <- "Small"> data$FamilyID <- factor(data$FamilyID)>

將新特徵變數FamilyID加入模型繼續預測。

> set.seed(1234)> model4 <- cforest(as.factor(Survived) ~ Pclass + Sex + Age + Fare + Embarked + Title + FamilySize + FamilyID + TicketCount, data = data[train.row, ], controls=cforest_unbiased(ntree=2000, mtry=3))> predict.result <- predict(model4, data[test.row, ], OOB=TRUE, type = "response")> submit <- data.frame(PassengerId = test$PassengerId, Survived = predict.result)> write.csv(submit, file = "cit4.csv", row.names = FALSE)>

新模型的預測結果在Kaggle上的得分為0.82297，排名190名，前190/7068=2.7%。

9 總結

本文主要參考網路上已有的一些資料和方法對著走了一遍流程，是比較深度的模仿，不過對我而言，也算是對數據分析、特別是分類類型的數據分析問題的一個總結。

從如何通過數據預覽，探索式數據分析，缺失數據填補，刪除關聯特徵以及新特徵挖掘等步驟完整地完成預定的預測，特別是隨機森林模型在二元分類問題上將其準確度高的優勢發揮的淋漓盡致。關於隨機森林模型的介紹可以參看下文：

機器學習演算法之隨機森林（Random Forest）

後續想繼續提高排名，提高預測的準確率，則需要構建一些新的變數或是構建新的模型。