導(dǎo)語(yǔ)：
1912年4月10日，號(hào)稱 “世界工業(yè)史上的奇跡”的豪華客輪泰坦尼克號(hào)開始了自己的處女航，從英國(guó)的南安普頓出發(fā)駛往美國(guó)紐約，4月14日晚，泰坦尼克號(hào)在北大西洋撞上冰山而傾覆，1502人葬生海底，705人得救。造成了當(dāng)時(shí)在和平時(shí)期最嚴(yán)重的一次航海事故，也是迄今為止最著名的一次海難。38歲的查爾斯·萊特勒是泰坦尼克二副，他是最后一個(gè)從冰冷的海水中被拖上救生船、職位最高的生還者。在他寫的回憶錄中，列舉幾個(gè)讓人震撼的情景：
- 在第一艘救生艇下水后，我對(duì)甲板上一名姓斯特勞的女人說：你能隨我一起到那只救生艇上去嗎？沒想到她搖了搖頭：不，我想還是呆在船上好。她的丈夫問：你為什么不愿意上救生艇呢？這名女人竟笑著回答：不，我還是陪著你。此后，我再也沒有見到過這對(duì)夫婦...
- 亞斯特四世（當(dāng)時(shí)世界第一首富）把懷著五個(gè)月身孕的妻子瑪?shù)铝账蜕?號(hào)救生艇后，站在甲板上，帶著他的狗，點(diǎn)燃一根雪茄煙，對(duì)劃向遠(yuǎn)處的小艇最后呼喊：我愛你們！一副默多克曾命令亞斯特上船，被亞斯特憤怒的拒絕：我喜歡最初的說法（保護(hù)弱者）！然后，把唯一的位置讓給三等艙的一個(gè)愛爾蘭婦女......
- 斯特勞斯是世界第二巨富，美國(guó)梅西百貨公司創(chuàng)始人。他無(wú)論用什么辦法，他的太太羅莎莉始終拒絕上八號(hào)救生艇，她說：多少年來，你去哪我去哪，我會(huì)陪你去你要去的任何地方。八號(hào)艇救生員對(duì)67歲的斯特勞斯先生提議：我保證不會(huì)有人反對(duì)像您這樣的老先生上小艇。斯特勞斯堅(jiān)定地回答：我絕不會(huì)在別的男人之前上救生艇。然后挽著63歲羅莎莉的手臂，一對(duì)老夫婦蹣姍地走到甲板的藤椅上坐下，等待著最后的時(shí)刻...
- 新婚燕爾的麗德帕絲同丈夫去美國(guó)渡蜜月，她死死抱住丈夫不愿獨(dú)自逃生，丈夫在萬(wàn) 般無(wú)奈中一拳將她打昏，麗德帕絲醒來時(shí)，她已在一條海上救生艇上了。此后，她終生未再嫁，以此懷念亡夫...
在這種生死存亡的緊要關(guān)頭，我們常常認(rèn)為社會(huì)等級(jí)越高、影響力越大，公眾認(rèn)可度越高的人物，生存的概率應(yīng)該越大，其次，乘客家庭成員多，成員間的協(xié)作和對(duì)求生的渴望度越高，生存的概率越高。然而，很多時(shí)候，事情產(chǎn)生這樣的結(jié)果的原因并非我們主觀臆測(cè)的那樣，我們需要通過對(duì)真實(shí)數(shù)據(jù)進(jìn)行科學(xué)的分析，才能發(fā)現(xiàn)很多事情并非我們想象的那樣簡(jiǎn)單，事情產(chǎn)生的本質(zhì)，往往隱藏在數(shù)據(jù)之中
下面我們就使用R語(yǔ)言根據(jù)已知存活情況的數(shù)據(jù)建立分析模型來預(yù)測(cè)其他一部分乘客的存活情況，其中，訓(xùn)練數(shù)據(jù)和測(cè)試數(shù)據(jù)均來源于：https://www.kaggle.com/c/titanic
本文的代碼和分析過程除部分修改外主體參考自Megan L. Risdal的文章：https://www.kaggle.com/mrisdal/exploring-survival-on-the-titanic
===========================第二次更新==============================
更新內(nèi)容:添加變量是否可以用于判斷存活預(yù)測(cè)的依據(jù)
1.乘客等級(jí)對(duì)生存率的影響
更新時(shí)間：2017年6月1號(hào)
=================================================================
一.數(shù)據(jù)的導(dǎo)入和查看
#有些包需要安裝,我們專門建立一個(gè)packagemanger.R文件來管理它門,在工程主入口文件中先進(jìn)行編譯后導(dǎo)入進(jìn)行使用
source('D:/R/RStudioWorkspace/titanic_test/utils/packageManager.R',encoding = 'UTF-8')
library(readr) # File read / write
library(ggplot2) # Data visualization
library(ggthemes) # Data visualization
library(scales) # Data visualization
library(plyr)
library(stringr) # String manipulation
library(InformationValue) # IV / WOE calculation
library(MLmetrics) # Mache learning metrics.e.g. Recall, Precision, Accuracy, AUC
library(rpart) # Decision tree utils
library(randomForest) # Random Forest
library(dplyr) # Data manipulation
library(e1071) # SVM
library(Amelia) # Missing value utils
library(party) # Conditional inference trees
library(gbm) # AdaBoost
library(class) # KNN
library(scales)
train <- read.csv('D:/R/RStudioWorkspace/Titanic dataset from Kaggle/train.csv',stringsAsFactors= FALSE)
test <- read.csv('D:/R/RStudioWorkspace/Titanic dataset from Kaggle/test.csv',stringsAsFactors= FALSE)
# 合并兩個(gè)數(shù)據(jù)框,查看相關(guān)變量名稱
total_data <- bind_rows(train,test)
str(total_data)
查看的數(shù)據(jù)結(jié)果如下:
'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" ...
$ 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" ...
我們觀察到一共有1309條數(shù)據(jù)，每一條數(shù)據(jù)有12個(gè)相關(guān)變量
$ PassengerId: 乘客編號(hào)
$ Survived :存活情況（存活：1 ; 死亡：0）
$ Pclass : 客場(chǎng)等級(jí)
$ Name : 乘客姓名
$ Sex : 性別
$ Age : 年齡
$ SibSp : 同乘的兄弟姐妹/配偶數(shù)
$ Parch : 同乘的父母/小孩數(shù)
$ Ticket : 船票編號(hào)
$ Fare : 船票價(jià)格
$ Cabin :客艙號(hào)
$ Embarked : 登船港口
二.特征工程
特征工程: 為了達(dá)到預(yù)測(cè)模型性能更佳，不僅要選取最好的算法，還要盡可能的從原始數(shù)據(jù)中獲取更多的信息。挖掘出更好的訓(xùn)練數(shù)據(jù)，就是特征工程建立的過程
2.1乘客社會(huì)等級(jí)越高,存活率越高
ggplot(total_data[!is.na(total_data$Survived),],aes(Pclass,fill=as.factor(Survived)))+
geom_bar(stat = 'count',position = 'dodge')+
xlab('Pclass')+ylab('Count')+
ggtitle(' how Pclass impact Survivor')+
scale_fill_manual(values=c("#FF6600", "#00C5CD")) +
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")
結(jié)果如下：
可以看到隨著乘客等級(jí)越低，在同一等級(jí)中的存活率越低，通過定量的計(jì)算Pclass的WOE（全稱是“Weight of Evidence”，即證據(jù)權(quán)重）和IV(Information Value，信息量,一個(gè)變量的IV 是一個(gè)可以定量衡量變量預(yù)測(cè)能力的指標(biāo),類似的指標(biāo)還有信息增益、基尼系數(shù)等等，可以參考這篇博客,詳細(xì)介紹了WOE和IV:http://blog.csdn.net/kevin7658/article/details/50780391)
可以算出Pclass的WOE和IV如下
WOETable(X=factor(total_data$Pclass[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
IV(X=factor(total_data$Pclass[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
#[1] 0.5009497
#attr(,"howgood")
#[1] "Highly Predictive"
從結(jié)果可以看出，Pclass的IV為0.5，且“Highly Predictive”，可以將Pclass 作為預(yù)測(cè)模型的特征變量
2.2 乘客頭銜Title 對(duì)生存率的影響
注意到在乘客名字（Name）中，有一個(gè)非常顯著的特點(diǎn)：乘客頭銜每個(gè)名字當(dāng)中都包含了具體的稱謂或者說是頭銜，將這部分信息提取出來后可以作為非常有用一個(gè)新變量，可以幫助我們進(jìn)行預(yù)測(cè)。此外也可以用乘客的姓代替家庭，生成家庭變量。
# 從名稱中挖掘
# 從乘客名字中提取頭銜
#R中的grep、grepl、sub、gsub、regexpr、gregexpr等函數(shù)都使用正則表達(dá)式的規(guī)則進(jìn)行匹配。默認(rèn)是egrep的規(guī)則，sub函數(shù)只實(shí)現(xiàn)第一個(gè)位置的替換，gsub函數(shù)實(shí)現(xiàn)全局的替換。
total_data$Title <- gsub('(.*, )|(\\..*)', '', total_data$Name)
# 查看按照性別劃分的頭銜數(shù)量
table(total_data$Sex, total_data$Title)
結(jié)果如下:
Capt Col Don Dona Dr Jonkheer Lady Major Master Miss Mlle Mme Mr Mrs Ms Rev Sir the Countess
female 0 0 0 1 1 0 1 0 0 260 2 1 0 197 2 0 0 1
male 1 4 1 0 7 1 0 2 61 0 0 0 757 0 0 8 1 0
我們發(fā)現(xiàn)頭銜的類別太多，并且好多出現(xiàn)的頻次是很低的，我們可以將這些類別進(jìn)行合并
# 合并低頻頭銜為一類
rare_title <- c('Dona', 'Lady', 'the Countess','Capt', 'Col', 'Don',
'Dr', 'Major', 'Rev', 'Sir', 'Jonkheer')
# 重命名稱呼
total_data$Title[total_data$Title == 'Mlle'] <- 'Miss'
total_data$Title[total_data$Title == 'Ms'] <- 'Miss'
total_data$Title[total_data$Title == 'Mme'] <- 'Mrs'
total_data$Title[total_data$Title %in% rare_title] <- 'Rare Title'
# 再次查看按照性別劃分的頭銜數(shù)量
table(total_data$Sex, total_data$Title)
得到如下結(jié)果：
Master Miss Mr Mrs Rare Title
female 0 264 0 198 4
male 61 0 757 0 25
下面來看看title 對(duì)生存率的影響，同樣的，使用圖形ggplot繪制
ggplot(total_data[!is.na(total_data$Survived),],aes(Title,fill=as.factor(Survived)))+
geom_bar(stat = 'count',position = 'dodge')+
ggtitle('how title impact survived')+
scale_fill_manual(values=c("#FF6600", "#00C5CD")) +
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")
得到如下結(jié)果
#查看Title 對(duì)Survived 的預(yù)測(cè)能力評(píng)估
WOETable(X=factor(total_data$Title[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
IV(X=factor(total_data$Title[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
#[1] 1.522418
#attr(,"howgood")
#[1] "Highly Predictive"
我們可以得出結(jié)論，Title 對(duì)survived有很好的預(yù)測(cè)效果，也需要把Title 添加到預(yù)測(cè)模型的特征變量中
最后，從名稱中獲取到姓氏
#sapply()函數(shù)：根據(jù)傳入?yún)?shù)規(guī)則重新構(gòu)建一個(gè)合理的數(shù)據(jù)類型返回
total_data$Surname <- sapply(total_data$Name, function(x) strsplit(x, split = '[,.]')[[1]][1])
2.3女性和小孩幸存概率應(yīng)該更大
作為弱者，女性和小孩在這種時(shí)刻應(yīng)該得到更好的照顧，生存率應(yīng)該會(huì)更高,
#性別對(duì)生存率的影響
ggplot(total_data[!is.na(total_data$Survived),],aes(Sex,fill=as.factor(Survived)))+
geom_bar(stat = 'count',position = 'dodge')+
ggtitle('how Sex impact survived')+
scale_fill_manual(values=c("#FF6600", "#00C5CD")) +
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")
得到如下結(jié)果：
同理我們可以得到性別的IV值為：1.341681 同樣也是“Highly Predictive”
#年齡對(duì)生存率的影響
#將年齡劃分成2個(gè)階段
total_data$AgeGroup[total_data$Age < 18] <- 'child'
total_data$AgeGroup[total_data$Age >= 18] <- 'adult'
table(total_data$AgeGroup,total_data$Survived)
ggplot(total_data[!is.na(total_data$Survived),],aes(x= AgeGroup,fill = as.factor(Survived)))+
geom_bar(stat = 'count',position = 'dodge')+
ggtitle('how AgeGroup impact survived')+
scale_fill_manual(values=c("#FF6600", "#00C5CD")) +
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")
得到如下結(jié)果
暫且不管NA （缺失數(shù)據(jù)）的存活情況，我們可以發(fā)現(xiàn)，小孩的成活概率是大于50%的
同樣的我們計(jì)算年齡組的IV值
WOETable(X=factor(total_data$AgeGroup[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
IV(X=factor(total_data$AgeGroup[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
#[1] 0.05655127
#attr(,"howgood")
#[1] "Somewhat Predictive"
發(fā)現(xiàn)預(yù)測(cè)能力為somewhat Predictive ：有些預(yù)測(cè)效果，暫且保留這個(gè)特征變量，到最后預(yù)測(cè)模型中對(duì)比加入和不加入這個(gè)變量對(duì)預(yù)測(cè)結(jié)果的影響大小再做結(jié)論
2.4 配偶及兄弟姐妹數(shù)適中的乘客更易幸存
我們來看看SibSp 這個(gè)變量對(duì)生存率的影響情況
#配偶及兄弟姐妹數(shù)適中的乘客更易幸存
ggplot(total_data[!is.na(total_data$Survived),],aes(x= SibSp,fill = as.factor(Survived)))+
geom_bar(stat = 'count',position = 'dodge')+
ggtitle('how sibsp impact survived')+
scale_fill_manual(values=c("#FF6600", "#00C5CD")) +
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")
得到如下結(jié)果
我們發(fā)現(xiàn)，配偶及兄弟姐妹數(shù)為1 或者2的生存率還是很高的，下面看看SibSp的IV值
WOETable(X=factor(total_data$SibSp[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
IV(X=factor(total_data$SibSp[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
#[1] 0.1448994
#attr(,"howgood")
#[1] "Highly Predictive"
高預(yù)測(cè)性，可以作為特征模型的一個(gè)預(yù)測(cè)變量
2.5 家庭成員數(shù)量的影響
既然我們已經(jīng)根據(jù)乘客的名字劃分成一些新的變量,我們可以把它進(jìn)一步做一些新的家庭變量。首先我們要做一個(gè)基于兄弟姐妹/配偶數(shù)量(s)和兒童/父母數(shù)量的家庭規(guī)模變量。
# 創(chuàng)建一個(gè)包含乘客自己的家庭規(guī)模變量
total_data$Fsize <- total_data$SibSp + total_data$Parch + 1
# Create a family variable
total_data$Family <- paste(total_data$Surname, total_data$Fsize, sep='_')
# 為了直觀顯示，我們可以用ggplot2 畫出家庭成員數(shù)量和生存家庭數(shù)情況的圖形
ggplot(total_data[!is.na(total_data$Survived),],aes(x= Fsize,fill = as.factor(Survived)))+
geom_bar(stat = 'count',position = 'dodge')+
ggtitle('how family size impact survived')+
scale_fill_manual(values=c("#FF6600", "#00C5CD")) +
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")
結(jié)果如下：
再來看看家庭成員數(shù)量對(duì)生存率的預(yù)測(cè)值IV
WOETable(X=factor(total_data$Fsize[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
IV(X=factor(total_data$Fsize[1:nrow(train)]),Y=total_data$Survived[1:nrow(train)])
[1] 0.3497672
attr(,"howgood")
[1] "Highly Predictive"
為高預(yù)測(cè)性
2.6支出船票價(jià)格對(duì)生存率的影響
船票價(jià)格是連續(xù)的，采用ggplot中 geom_line()進(jìn)行模擬顯示
#支出船票價(jià)格對(duì)生存率的影響
ggplot(total_data[1:nrow(train), ], aes(x = Fare, fill= as.factor(Survived),color = Survived)) +
geom_line(aes(label=..count..), stat = 'bin', binwidth=10) +
scale_fill_manual(values=c("#FF6600", "#00C5CD")) +
labs(title = "How Fare impact survivor", x = "Fare", y = "Count", fill = "Survived")
結(jié)果如下
觀察藍(lán)色存活數(shù)量的線條，我們可以發(fā)現(xiàn)，船票價(jià)格越高，生存率越高，我們?cè)賮砜纯创眱r(jià)格的IV值，結(jié)果為高預(yù)測(cè)性
#[1] 0.6123083
#attr(,"howgood")
#[1] "Highly Predictive"
2.7 客艙位置的影響
可以發(fā)現(xiàn)在乘客客艙變量 passenger cabin 也存在一些有價(jià)值的信息如客艙層數(shù) deck,但是這個(gè)變量的缺失值太多，無(wú)法做出新的有效的變量，暫時(shí)放棄這個(gè)變量的挖掘
三.缺失數(shù)據(jù)的處理
觀察文件中的數(shù)據(jù)，我們會(huì)發(fā)現(xiàn)有些乘客的信息參數(shù)并不完整，由于所給的數(shù)據(jù)集并不大，我們不能通過刪除一行或者一列來處理缺失值，因而對(duì)于我們關(guān)注的一些字段參數(shù)，我們需要根據(jù)統(tǒng)計(jì)學(xué)的描述數(shù)據(jù)（平均值、中位數(shù)等等）來合理給出缺失值
3.1 列出所有缺失值
我們可以通過函數(shù)查看缺失數(shù)據(jù)的變量在第幾條數(shù)據(jù)出現(xiàn)缺失和總共缺失的個(gè)數(shù)
3.1 年齡的缺失和填補(bǔ)
#統(tǒng)計(jì)年齡的缺失個(gè)數(shù)
age_null_count <- sum(is.na(total_data$Age))
#age_null_count = 263
通常我們會(huì)使用 rpart (recursive partitioning for regression) 包來做缺失值預(yù)測(cè) 在這里我將使用 mice 包進(jìn)行處理。我們先要對(duì)因子變量（factor variables）因子化，然后再進(jìn)行多重插補(bǔ)法。
#統(tǒng)計(jì)年齡的缺失處理
age_null_count <- sum(is.na(total_data$Age))
# 使自變量因子化
factor_vars <- c('PassengerId','Pclass','Sex','Embarked',
'Title','Surname','Family','Fsize')
#lapply()返回一個(gè)長(zhǎng)度與X一致的列表，每個(gè)元素為FUN計(jì)算出的結(jié)果，且分別對(duì)應(yīng)到X中的每個(gè)元素。
total_data[factor_vars] <- lapply(total_data[factor_vars],function(x) as.factor(x))
# 設(shè)置隨機(jī)值
set.seed(129)
# 執(zhí)行多重插補(bǔ)法，剔除一些沒什么用的變量:
mice_mod <- mice(total_data[, !names(total_data) %in% c('PassengerId','Name','Ticket','Cabin','Family','Surname','Survived')], method='rf')
# 保存完成的輸出
mice_output <- complete(mice_mod)
讓我們來比較一下我們得到的結(jié)果與原來的乘客的年齡分布以確保沒有明顯的偏差
# 繪制直方圖
par(mfrow=c(1,2))
hist(total_data$Age, freq=F, main='Age: Original Data',
col='darkgreen', ylim=c(0,0.04))
hist(mice_output$Age, freq=F, main='Age: MICE Output',
col='lightgreen', ylim=c(0,0.04))
結(jié)果如下,右邊圖和左邊圖有很高的相似度
所以，我們可以用mice模型的結(jié)果對(duì)原年齡數(shù)據(jù)進(jìn)行替換。
# 用mice模型數(shù)據(jù)替換原始數(shù)據(jù)
full$Age <- mice_output$Age
# 再次查看年齡的缺失值數(shù)據(jù)
sum(is.na(full$Age))
# 0
3.2票價(jià)的缺失處理
#查看票價(jià)的缺失值
getFareNullID <- function(total_data){
count <- 0
for(i in 1:nrow(total_data))
if(is.na(total_data$Fare[i])){
#打印缺失票價(jià)的具體行數(shù)
print(i);
count <- count+1
}
return(count)
}
fare_null_count <- getFareNullID(total_data)
#fare_null_count = 1
得到票價(jià)缺失個(gè)數(shù)為1 ，缺失行數(shù)為第1044行
查看這一行我們會(huì)發(fā)現(xiàn)
total_data[1044,]
PassengerId Survived Pclass Name Sex Age SibSp Parch Ticket Fare Cabin
1044 1044 NA 3 Storey, Mr. Thomas male 60.5 0 0 3701 NA
Embarked Title Surname Fsize Family
1044 S Mr Storey 1 Storey_1
我們發(fā)現(xiàn)港口和艙位是完整的，我們可以根據(jù)相同的港口和相同的艙位來大致估計(jì)該乘客的票價(jià)，我們?nèi)∵@些類似乘客的中位數(shù)來替換缺失的值
#從港口Southampton ('S')出發(fā)的三等艙乘客。 從相同港口出發(fā)且處于相同艙位的乘客數(shù)目
same_farenull <- sum(total_data$Pclass == '3' & total_data$Embarked == 'S')
# 基于出發(fā)港口和客艙等級(jí)，替換票價(jià)缺失值
total_data$Fare[1044] <- median(total_data[total_data$Pclass == '3' & total_data$Embarked == 'S', ]$Fare, na.rm = TRUE)
3.3登船港口號(hào)的缺失
#登船港口號(hào)的缺失值函數(shù)
getEmbarkedNullCount <- function(total_data) {
count0 <- 0
count <- 0
for(i in 1:nrow(total_data))
if(total_data$Embarked[i] == ""){
#可以打印出缺失的所在行數(shù)
print(i);
count <- count +1
}
return(count)
}
#登船港口號(hào)的缺失個(gè)數(shù)
embarked_null_count <- getEmbarkedNullCount(total_data)
#embarked_null_count =2
得到登船港口號(hào)缺失的個(gè)數(shù)為2 ，分別為 62 、830，我估計(jì)對(duì)于有相同艙位等級(jí)（passenger class）和票價(jià)（Fare）的乘客也許有著相同的 登船港口位置embarkment .我們可以看到他們支付的票價(jià)分別為: $ 80 和 $ 80 同時(shí)他們的艙位等級(jí)分別是： 1 和 1 . 我們可以用箱線圖繪制出這三者之間關(guān)系圖
從港口 ('C')出發(fā)的頭等艙支付的票價(jià)的中位數(shù)正好為80。因此我們可以放心的把處于頭等艙且票價(jià)在$80的乘客62和830 的出發(fā)港口缺失值替換為'C'
total_data$Embarked[c(62, 830)] <- 'C'
我們基本上完成了重要參數(shù)缺失值的處理，我們的數(shù)據(jù)集變得更加完整了呢，接下來，需要根據(jù)新的數(shù)據(jù)集創(chuàng)建出新的特征工程
四.新特征工程的建立
通過上面缺失值填補(bǔ)的完成，我們?cè)囍谛碌臄?shù)據(jù)集中挖掘出對(duì)乘客的存活有影響的一些因素，根據(jù)文章剛開始的幾段真實(shí)場(chǎng)景預(yù)測(cè)，我們考慮在這種災(zāi)難性的時(shí)刻，小孩和老人相對(duì)于青年或者中年人應(yīng)該會(huì)得到更好的照顧，生存的概率應(yīng)該更高，其次，如果你是一位母親，你相比于其他成年女性是否會(huì)有更高的存活可能？其實(shí)我還有一個(gè)想法，那就是乘客的社會(huì)地位或者說階層 和當(dāng)時(shí)的收入水平層次可能對(duì)生存有一定的影響，當(dāng)然這兩個(gè)因素對(duì)于現(xiàn)在的我們來說非常難以獲取，畢竟事情發(fā)生在100多年前，或許當(dāng)時(shí)的政府，也需要很長(zhǎng)時(shí)間才能準(zhǔn)確的獲取到覺大部分人的這些側(cè)面信息。
4.1年齡的劃分
我們考慮將年齡劃分成三個(gè)階段，小于18歲算小孩，18歲及以上至50歲為青壯年，50歲以上為老年人
#將年齡劃分成3個(gè)階段
total_data$AgeGroup[total_data$Age < 18] <- 'child'
total_data$AgeGroup[total_data$Age >= 18 & total_data$Age <= 50] <- 'young'
total_data$AgeGroup[total_data$Age > 50] <- 'old'
table(total_data$AgeGroup,total_data$Survived)
#得到如下結(jié)果
0 1
child 70 63
old 51 24
young 428 255
相比于成人，小孩的生存概率接近50%，小孩得到的照顧比成年高的多
4.2是否為母親
我們從性別和頭銜中提煉出一位成年女性是否為一位母親，看看她的生存概率如何
# Adding Mother variable
total_data$IsMother <- 'Not'
total_data$IsMother[total_data$Sex == 'female' & total_data$Parch > 0 & total_data$Age > 18 & total_data$Title != 'Miss'] <- 'Yes'
# Show counts
table(total_data$IsMother, total_data$Survived)
#結(jié)果如下：
0 1
Not 534 303
Yes 15 39
我們發(fā)現(xiàn)，如果是一位母親，那么你生存下來的概率高達(dá)70%,之后，我們整合上面兩個(gè)新變量到原數(shù)據(jù)集
# 完成因子化
total_data$AgeGroup <- factor(total_data$AgeGroup)
total_data$IsMother <- factor(total_data$IsMother)
#mice 包中顯示缺失數(shù)據(jù)的一種模式。
md.pattern(total_data)
五.預(yù)測(cè)
到了最激動(dòng)人心的時(shí)刻了有沒有，前面四個(gè)步驟都是為了預(yù)測(cè)在做前期準(zhǔn)備，如何進(jìn)行預(yù)測(cè)呢？
5.1.拆分測(cè)試和訓(xùn)練數(shù)據(jù)集
#拆分?jǐn)?shù)據(jù)集
train <- total_data[1:891,]
test <- total_data[892:1309,]
5.2 構(gòu)建訓(xùn)練模型
我們使用隨機(jī)森林法則作用于訓(xùn)練數(shù)據(jù)集來構(gòu)建我們需要的預(yù)測(cè)模型
#拆分?jǐn)?shù)據(jù)集
train <- total_data[1:891,]
test <- total_data[892:1309,]
set.seed(754)
# 構(gòu)建預(yù)測(cè)模型
rf_model <- randomForest(factor(Survived) ~ Pclass + Sex + Age + Fare+ Embarked + Title + Fsize,data = train)
交叉驗(yàn)證
一般情況下，應(yīng)該將訓(xùn)練數(shù)據(jù)分為兩部分，一部分用于訓(xùn)練，另一部分用于驗(yàn)證?；蛘呤褂胟-fold交叉驗(yàn)證。本文將所有訓(xùn)練數(shù)據(jù)都用于訓(xùn)練，然后隨機(jī)選取30%數(shù)據(jù)集用于驗(yàn)證。
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(415)
cv.test.sample <- sample(1:nrow(train), as.integer(0.3 * nrow(train)), replace = TRUE)
cv.test <- total_data[cv.test.sample,]
cv.prediction <- predict(rf_model, cv.test, OOB=TRUE, type = "response")
cv.summarize(cv.test$Survived, cv.prediction)
#"Recall: 0.982658959537572"
#[1] "Precision: 0.904255319148936"
#[1] "Accuracy: 0.921348314606742"
#[1] "AUC: 0.895584798917722"
5.3 相關(guān)性檢測(cè)
通過隨機(jī)森林中所有決策樹的Gini 計(jì)算出其他變量相對(duì)于生存變量的相關(guān)性排行，我們可以看出那些因素對(duì)生存率影響較大
# 重要性系數(shù)
importance <- importance(rf_model)
varImportance <- data.frame(Variables = row.names(importance),
Importance = round(importance[ ,'MeanDecreaseGini'],2))
# 創(chuàng)建基于重要性系數(shù)排列的變量
rankImportance <- varImportance %>%
mutate(Rank = paste0('#',dense_rank(desc(Importance))))
# 使用 ggplot2 繪出重要系數(shù)的排名
ggplot(rankImportance, aes(x = reorder(Variables, Importance),
y = Importance, total_data = Importance)) +
geom_bar(stat='identity') +
geom_text(aes(x = Variables, y = 0.5, label = Rank),
hjust=0, vjust=0.55, size = 4, colour = 'red') +
labs(x = 'Variables') +
coord_flip() +
theme_few()
結(jié)果如下：
5.4 預(yù)測(cè)
最后，我們使用訓(xùn)練好的特征模型作用于測(cè)試數(shù)據(jù)上，得到我們的預(yù)測(cè)結(jié)果
prediction <- predict(rf_model, test)
# 保存數(shù)據(jù)結(jié)果passagerId 和survived參數(shù)
solution <- data.frame(PassengerID = test$PassengerId, Survived = prediction)
# 保存到文件
write.csv(solution, file = 'D:/R/RStudioWorkspace/titanic_test/output/predict_Solution.csv', row.names = F)
得到預(yù)測(cè)結(jié)果文件，我們可以上傳到Kaggle,查看自己的排名情況，
第二次預(yù)測(cè)進(jìn)行了特征變量的刪減，刪了 AgeGroup 和 IsMother,SibSp,，我們第一次選擇特征變量的時(shí)候認(rèn)為的小孩、老人和是否為母親 這幾個(gè)特征應(yīng)該有很大的生存幾率，但是結(jié)果并不是這樣，現(xiàn)實(shí)還是比較殘酷！
就先分析到這吧，感謝你的時(shí)間，后面靈感涌現(xiàn)挖掘到新的特征變量再添加到特征工程中，這樣預(yù)測(cè)結(jié)果應(yīng)該會(huì)更加準(zhǔn)確。排名也會(huì)更加靠前，加油！
