R语言实战-Kaggle-泰坦尼克号存活率分析(Titanic:Machine Learning from Disaster)
张英龙(Alvin Zhang)
2018年3月23日
项目背景
(翻译自kaggle)
泰坦尼克号的沉没是历史上最著名的沉船事件之一。1912年4月15日，泰坦尼克号在她的初次航行中，撞上冰山后沉没，导致2224名乘客和船员中有1502人丧生。这起轰动的惨剧震惊了国际社会，并因此建立了更好的船舶安全准则。
这起沉船事故造成如此规模的人员不幸遇难的原因之一是船上没有足够的救生艇供乘客和船员使用。尽管在沉船事故中幸存下来的人有一些运气成分，但有些人比其他人更容易生还，比如妇女、儿童和上层阶级。
项目目的
本文诣在分析泰坦尼克号上哪类人存活几率更大，并应用随机森林来预测哪些乘客能够存活。
一、数据导入
(项目数据集可以在kaggle上进行下载)
数据集简介 项目所提供的数据被分成两组：
训练集 (train.csv) 用于数据分析和构建机器学习模型，提供每位乘客是否存活的结果（提供survive变量）。
测试集 (test.csv) 用来查看模型在看不见的数据上的表现，不提供每个乘客的是否存活的结果（不提供survive变量）。
数据集中12个变量含义简介如下：
- PassengerId：乘客的身份编号，递增变量，对预测无帮助
- Survived：乘客是否幸存，0 = No（遇难）, 1 = Yes（幸存）
- Pclass：乘客的船舱等级，也是乘客社会经济地位的代表，1 = 1st(头等舱)=Upper（上流社会人士）, 2 = 2nd（二等舱）=Middle（中产阶层人士）, 3 = 3rd（三等舱）=Lower（底层社会人士）
- Name：乘客姓名
- Sex：乘客性别
- Age：乘客年龄
- SibSp：乘客的兄弟姐妹的个数及是否有配偶，Sib=Sibling（兄弟姐妹），Sp=Spouse（配偶）
- Parch：乘客的父母或子女的个数。Par=Parent（父母），Ch=Child（子女）
- Ticket：乘客的船票号
- Fare：乘客的船票价
- Cabin：乘客的舱位号
- Embarked：乘客开始出发港，C = Cherbourg, Q = Queenstown, S = Southampton
#加载包
library(readr) #读取数据
library(dplyr) #数据处理
library(ggplot2) #数据可视化
library(VIM) #缺失值处理
library(rpart) #回归树方法，建立模型预测缺失值
library(InformationValue) #计算数据的预测价值
library(stringr) #字符串处理
library(randomForest) #建立预测模型
#导入数据,这里使用绝对路径，可根据数据集所在位置进行修改
train_data <- read_csv("G:/Study_Resource/R/kaggle/titanic/train.csv") #导入训练集
test_data <- read_csv("G:/Study_Resource/R/kaggle/titanic/test.csv") #导入测试集二、数据预处理
#将两个数据集合并，方便进行数据预处理
all_data <- bind_rows(train_data, test_data)
#查看数据集数据类型等信息
str(all_data)## Classes 'tbl_df', 'tbl' and '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 NA "C85" NA "C123" ...
## $ Embarked : chr "S" "C" "S" "S" ...可以看出，合并之后的数据集共包含12个变量，1309条数据，其中，891条为训练数据，418条为测试数据。
#先将所有空值转变为缺失值，便于查看
all_data[all_data[1:nrow(all_data), ] == ""] <- NA
#查看哪些变量存在缺失值
aggr(all_data,plot=FALSE)##
## Missings in variables:
## Variable Count
## Survived 418
## Age 263
## Fare 1
## Cabin 1014
## Embarked 2由结果可知：
- Survived 变量缺失418条数据
- Fare 变量缺失1条数据
- Age 变量缺失263条数据
- Cabin 变量缺失1014条数据
- Embarked 变量缺失2条数据
缺失值处理
1.Survived 变量
此变量缺失值数量与测试集数据数量相同，都是需要预测的变量，所以缺失值无需补充
- Fare变量（缺失1条数据）
#查看Fare变量缺失值所在的行数
which(is.na(all_data$Fare))## [1] 1044#提取并查看该行所有变量信息
data_1044 <- all_data[1044, ]
show(data_1044)## # A tibble: 1 x 12
## PassengerId Survived Pclass Name Sex Age SibSp Parch
## <int> <int> <int> <chr> <chr> <dbl> <int> <int>
## 1 1044 NA 3 Storey, Mr. Thomas male 60.5 0 0
## # ... with 4 more variables: Ticket <chr>, Fare <dbl>, Cabin <chr>,
## # Embarked <chr>可以看出该名乘客的Pclass（船舱等级）为3，Embarked（开始出发港）为S，根据数据集中符合这两个条件的其他数据情况来推断并补充这个缺失值。
#提取所有满足条件的乘客们的信息并绘图
process_data <- subset(all_data, Pclass==3 & Embarked=="S")
ggplot(process_data, aes(x=Fare)) +
geom_density(fill="skyblue", alpha=.3) +
geom_vline(aes(xintercept=median(Fare, na.rm=TRUE)),
color="red", linetype=2, size=.5)
由图可知，满足Pclass=3并且Embarked=S这两个条件的乘客们的Fare变量的值集中在0~20之间，本文使用满足上述两个条件的Fare变量的中位数来替代缺失值。
#使用Fare变量中位数替代Fare变量缺失值
all_data$Fare[is.na(all_data$Fare)] <- median(process_data$Fare,na.rm=TRUE)- Age变量（缺失263条数据）
#使用除了Age变量之外的其他变量建立模型
age_model <- rpart(Age ~ Pclass + Sex + SibSp + Parch + Fare + Embarked,
data=all_data[!is.na(all_data$Age), ], method="anova")
#通过模型预测Age变量的缺失值并将缺失值补充完整
all_data$Age[is.na(all_data$Age)] <- predict(age_model, all_data[is.na(all_data$Age), ])- Cabin变量（缺失1014条数据） 数据集共有1309条数据，此变量缺失数据过多，暂定将缺失值填充为“unknown”
#将Cabin变量缺失值填充为“unknown”
all_data$Cabin[is.na(all_data$Cabin)] <- "unknown"- Embarked变量（缺失2条数据）
#查看Embarked变量缺失值所在的行数
which(is.na(all_data$Embarked))## [1] 62 830#提取并查看这两个缺失值所在行的所有变量信息
data_62 <- all_data[62, ]
show(data_62)## # A tibble: 1 x 12
## PassengerId Survived Pclass Name Sex Age SibSp Parch
## <int> <int> <int> <chr> <chr> <dbl> <int> <int>
## 1 62 1 1 Icard, Miss. Amelie female 38 0 0
## # ... with 4 more variables: Ticket <chr>, Fare <dbl>, Cabin <chr>,
## # Embarked <chr>data_830 <- all_data[830, ]
show(data_830)## # A tibble: 1 x 12
## PassengerId Survived Pclass Name
## <int> <int> <int> <chr>
## 1 830 1 1 Stone, Mrs. George Nelson (Martha Evelyn)
## # ... with 8 more variables: Sex <chr>, Age <dbl>, SibSp <int>,
## # Parch <int>, Ticket <chr>, Fare <dbl>, Cabin <chr>, Embarked <chr>可以看出这两名乘客的Pclass（船舱等级）均为1，Fare（船票价）均为80，绘制根据Pclass变量进行分类的Fare变量与Embarked关系的箱线图。
#提取数据并绘图
process_data <- all_data[!is.na(all_data$Embarked), ]
process_data$Pclass <- as.factor(process_data$Pclass)
ggplot(process_data, aes(x=Embarked, y=Fare, fill=Pclass)) +
geom_boxplot() +
geom_hline(aes(yintercept=80), color="red", linetype=2, size=.5)
由图可知，满足Pclass=1并且Embarked=C这两个条件的乘客们的Fare≈80，据此本文使用C替换Embarked变量的两个缺失值。
#使用“C”替代Embarked变量的缺失值
all_data$Embarked[is.na(all_data$Embarked)] <- "C"特征重构
1.Name 变量
首先，乘客的姓名对预测的意义不大，但是Name变量中还包括了一些乘客的称谓（例如Mr.和Miss等等），将这些称谓提取出来存放进新的变量Title中。
#提取乘客称谓
all_data$Title <- gsub("(.*, )|(\\..*)","",all_data$Name)
#将称谓按性别分类
table(all_data$Sex, all_data$Title)##
## Capt Col Don Dona Dr Jonkheer Lady Major Master Miss Mlle Mme
## female 0 0 0 1 1 0 1 0 0 260 2 1
## male 1 4 1 0 7 1 0 2 61 0 0 0
##
## Mr Mrs Ms Rev Sir the Countess
## female 0 197 2 0 0 1
## male 757 0 0 8 1 0#将出现频率较少的稀有称谓合并
rare_title <- c('Capt', 'Col', 'Don', 'Dona', 'Dr', 'Jonkheer',
'Major', 'Mlle','Mme', 'Rev', 'the Countess')
all_data$Title[all_data$Title == 'Lady'] <- 'Miss'
all_data$Title[all_data$Title == 'Ms'] <- 'Miss'
all_data$Title[all_data$Title == 'Sir'] <- 'Mr'
all_data$Title[all_data$Title %in% rare_title] <- 'Rare Title'2.SibSp 变量和Parch 变量
通过观察SibSp 和Parch变量的含义，不难发现两者都是家庭成员数的一部分，将它们概括为Familysize作为新的变量使用，Familysize变量定义如下：
Familysize=SibSp+Parch+1
（如果SibSp 和Parch变量均为0，则Familysize=1表示该乘客独自乘船）
#将SibSp和Parch变量合并为新的变量Familysize
all_data$Familysize <- all_data$SibSp + all_data$Parch + 13.Ticket 变量
通过观察数据，可以发现Ticket 变量有重复取值的情况，将该变量取值唯一的数据归类为“Unique”,重复取值的数据归类为“Share”
#统计每张票对应的乘客数
attach(all_data)
ticket_count <- aggregate(Ticket, by=list(Ticket), function(x) sum(!is.na(x)))
detach(all_data)
#将数据存入all_data数据集新建变量Tickettype中
all_data$Tickettype <- apply(all_data, 1,
function(x) ticket_count[which(ticket_count[ ,1]==x['Ticket']),2])
#将TicketCount=1和>1的数据分别标注为Unique和Share
all_data$Tickettype <- factor(sapply(all_data$Tickettype,
function(x) ifelse(x>1, "Share", "Unique")))三、数据分析
分别分析各变量对Survived变量（存活率）的影响
1.Pclass变量对Survived变量（存活率）的影响
#将Survived变量类型转换为因子
all_data$Survived <- factor(all_data$Survived)
#绘制Pclass-Survived分布图
ggplot(data=all_data[1:nrow(train_data), ], aes(x=Pclass, y=..count.., fill=Survived)) +
geom_bar(position="dodge") +
geom_text(stat = "count", aes(label=..count..), position=position_dodge(.9), vjust=-.3)
由图可知，Pclass=1的乘客幸存率最高（63.0%），幸存人数也是最多的（136人），Pclass=2的乘客幸存率也较高（47.3%），Pclass=3的乘客幸存率最低（24.2%）。据此本文认为船舱等级越高，存活几率就越大。为了更定量的计算Pclass变量的预测价值，可以算出Pclass的IV如下:
#计算IV判断Pclass变量的预测价值
IV(X=factor(all_data$Pclass[1:nrow(train_data)]),
Y=all_data$Survived[1:nrow(train_data)])## [1] 0.5009497
## attr(,"howgood")
## [1] "Highly Predictive"由结果可知，Pclass变量的IV=0.5009497，且“Highly Predictive”。因此，可将Pclass变量作为预测模型的特征变量之一。
2.Title变量对Survived变量（存活率）的影响
#绘制Title-Survived分布图
ggplot(data = all_data[1:nrow(train_data), ], aes(x=Title, y =..count.., fill=Survived)) +
geom_bar(position="dodge") +
geom_text(stat = "count", aes(label=..count..), position=position_dodge(.9), vjust=-.3)
由图可知，Title（乘客称谓）=“Mr”的乘客存活率较低（15.8%），而Title=“Master”、“Miss”、“Mrs”的乘客存活率均大于50%。据此本文认为称谓为“Master”、“Miss”、“Mrs”的乘客存活几率大。同样，使用IV来定量计算Title变量对于最终的预测是否有用。
#计算IV判断Title变量的预测价值
IV(X=factor(all_data$Title[1:nrow(train_data)]),
Y=all_data$Survived[1:nrow(train_data)])## [1] 1.49771
## attr(,"howgood")
## [1] "Highly Predictive"由结果可知，Title变量IV=1.49771，且”Highly Predictive”。因此，可将Title变量作为预测模型的特征变量之一。
3.Sex变量对Survived变量（存活率）的影响
#将Sex变量类型转换为因子
all_data$Sex <- as.factor(all_data$Sex)
#绘制Sex-Survived分布图
ggplot(data=all_data[1:nrow(train_data), ], aes(x=Sex, y=..count.., fill=Survived)) +
geom_bar(position="dodge") +
geom_text(stat="count", aes(label=..count..), position=position_dodge(.9), vjust=-.3)
由图可知，Sex=“female”的乘客幸存率（74.2%）远大于Sex=“male”的乘客幸存率（18.9%），女性幸存人数约为男性的2倍。据此本文认为女性存活几率比男性大。
#计算IV判断Sex变量的预测价值
IV(X=factor(all_data$Sex[1:nrow(train_data)]),
Y=all_data$Survived[1:nrow(train_data)])## [1] 1.341681
## attr(,"howgood")
## [1] "Highly Predictive"由结果可知，Sex变量的IV为1.341681且“Highly Predictive”。因此，可将Sex变量也作为预测模型的特征变量之一。
4.Age变量对Survived变量（存活率）的影响
#绘制Age-Survived分布图
ggplot(data=all_data[1:nrow(train_data), ], aes(x=Age, color=Survived)) +
geom_line(stat = "bin", binwidth=5)
由图可知，Age<18的乘客幸存率（62.0%）略大于Age<18的乘客幸存率（57.5%）。据此本文认为未成年人存活几率比成年人大。
#计算IV判断Age变量的预测价值
IV(X=factor(all_data$Age[1:nrow(train_data)]),
Y=all_data$Survived[1:nrow(train_data)])## [1] 0.3146341
## attr(,"howgood")
## [1] "Highly Predictive"由结果可知，Age变量的IV=0.3146341且“Highly Predictive”。因此，可将Age变量也作为预测模型的特征变量之一。
5.Familysize变量对Survived变量（存活率）的影响
#绘制Familysize-Survived分布图
ggplot(data=all_data[1:nrow(train_data), ],
aes(x=as.factor(Familysize), y=..count.., fill=Survived)) +
geom_bar(position="dodge") +
geom_text(stat = "count", aes(label=..count..),
position=position_dodge(.9),vjust=-.3)
由图可知，Familysize=2或3或4的乘客幸存率均大于50%，远大于Familysize为其他值的乘客。据此本文认为乘客的家庭成员数为2或3或4的乘客存活几率比较大。
#计算IV判断Familysize变量的预测价值
IV(X=factor(all_data$Familysize[1:nrow(train_data)]),
Y=all_data$Survived[1:nrow(train_data)])## [1] 0.3497672
## attr(,"howgood")
## [1] "Highly Predictive"由结果可知，Familysize变量的IV=0.3497672且“Highly Predictive”。因此，可将Familysize变量也作为预测模型的特征变量之一。
6.Tickettype变量对Survived变量（存活率）的影响
#绘制Tickettype-Survived分布图
ggplot(data=all_data[1:nrow(train_data), ],
aes(x=Tickettype, y=..count.., fill=Survived)) +
geom_bar(position="dodge") +
geom_text(stat = "count", aes(label=..count..),
position=position_dodge(.9),vjust=-.3)
由图可知，共享票号的乘客幸存率大于50%，远大于独立票号的乘客。据此本文认为共享票号的乘客存活几率比较大。
#计算IV判断Tickettype变量的预测价值
IV(X=factor(all_data$Tickettype[1:nrow(train_data)]),
Y=all_data$Survived[1:nrow(train_data)])## [1] 0.2751882
## attr(,"howgood")
## [1] "Highly Predictive"由结果可知，Tickettype变量的IV=0.2751882且“Highly Predictive”。因此，可将Tickettype变量也作为预测模型的特征变量之一。
7.Fare变量对Survived变量（存活率）的影响
#绘制Fare-Survived分布图
ggplot(data=all_data[1:nrow(train_data), ], aes(x=Fare, color=Survived)) +
geom_line(stat="bin", binwidth=10)
由图可知，船票价越高，存活几率越大。
#计算IV判断Fare变量的预测价值
IV(X=factor(all_data$Fare[1:nrow(train_data)]),
Y=all_data$Survived[1:nrow(train_data)])## [1] 0.6123083
## attr(,"howgood")
## [1] "Highly Predictive"由结果可知，Fare变量的IV=0.6123083且“Highly Predictive”。因此，可将Fare变量也作为预测模型的特征变量之一。
8.Cabin变量对Survived变量（存活率）的影响
#绘制Cabin-Survived分布图
ggplot(data=all_data[1:nrow(train_data), ],
aes(x=as.factor(sapply(all_data$Cabin[1:nrow(train_data)],
function(x) str_sub(x, end=1))),
y =..count.., fill = Survived)) +
geom_bar(position="dodge") +
geom_text(stat="count", aes(label=..count..), position=position_dodge(.9), vjust=-.3)
由图可知，舱位号=“B”、“C”、“D”、“E”、“F”的乘客存活率均大于50%，据此本文认为舱位号为“B”、“C”、“D”、“E”、“F”的乘客更容易存活。由于Cabin变量缺失值较多，暂不将其作为预测模型的特征变量。
9.Embarked变量对Survived变量（存活率）的影响
#绘制Embarked-Survived分布图
ggplot(data=all_data[1:nrow(train_data), ], aes(x=Embarked, y =..count.., fill=Survived)) +
geom_bar(position="dodge") +
geom_text(stat="count", aes(label=..count..), position=position_dodge(.9), vjust=-.3)
由图可知，开始出发港=“C”的乘客存活率最大（55.9%），开始出发港=“S”的乘客存活率最小（33.7%）据此本文认为开始出发港为“Cherbourg”的乘客存活几率大。
#计算IV判断Embarked变量的预测价值
IV(X=factor(all_data$Embarked[1:nrow(train_data)]),
Y=all_data$Survived[1:nrow(train_data)])## [1] 0.129404
## attr(,"howgood")
## [1] "Highly Predictive"由结果可知，Embarked变量的IV=0.129404且“Highly Predictive”。因此，可将Embarked变量也作为预测模型的特征变量之一。
四、建立预测模型
#查看数据类型
str(all_data)## Classes 'tbl_df', 'tbl' and 'data.frame': 1309 obs. of 15 variables:
## $ PassengerId: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Survived : Factor w/ 2 levels "0","1": 1 2 2 2 1 1 1 1 2 2 ...
## $ 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 : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...
## $ Age : num 22 38 26 35 35 ...
## $ 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 "unknown" "C85" "unknown" "C123" ...
## $ Embarked : chr "S" "C" "S" "S" ...
## $ Title : chr "Mr" "Mrs" "Miss" "Mrs" ...
## $ Familysize : num 2 2 1 2 1 1 1 5 3 2 ...
## $ Tickettype : Factor w/ 2 levels "Share","Unique": 2 1 2 1 2 2 1 1 1 1 ...#将特征变量中字符型变量转化为因子
all_data$Title <- as.factor(all_data$Title)
all_data$Embarked <- as.factor(all_data$Embarked)
#建立预测模型
set.seed(1234)
model <- randomForest(Survived ~ Pclass + Title + Sex + Age + Familysize +
Tickettype + Fare + Embarked, data=all_data[1:nrow(train_data), ])
#对test_data数据集Survived变量进行预测
prediction <- predict(model, all_data[(1+nrow(train_data)):(nrow(all_data)), ])
#输出结果
result <- data.frame(PassengerId = test_data$PassengerId,
Survived = prediction)
#将结果写入文件
write.csv(result, file="G:/Study_Resource/R/kaggle/titanic/prediction.csv", row.names = FALSE)将结果上传至Kaggle，排名3171，前30%
五、总结
本项目为作者第一次尝试使用R语言进行一次完整的数据分析流程，参考了知乎上一些其他同学的资料。总的来说，学习到了很多东西，也对书本中和网络上学到的知识有了更深的理解。通过这次实践也发现了自己很多的不足之处，有待提高。