【Kaggle】泰坦尼克号生存训练
2017-11-08 本文已影响69人
奔跑的蜈蚣
就是这篇文章，知乎称“您的帐号由于存在异常行为暂时被知乎反作弊系统限制使用”，然后任凭我申诉多久，都恢复不了了！！！最可恶的是，在你发布文章的时候一点儿提示都没有，显示发布成功，但就是擅自删除，辛苦写的文章就找不到了！！！
1.引言
本文数据来自Kaggle中知名的数据集Titanic Machine Learning from Disaster，是利用训练集训练模型，来预测测试集中的乘客能否在沉船事件存活。
说明：本文借鉴了 Megan L. Risdal 的文章，在细节处略有改动。
2.数据导入与观察
加载需要用到的包：
> pkgs <- c("dplyr","ggplot2","ggthemes","scales","mice","randomForest")
> install.packages(pkgs,dependencies = TRUE)
> library(dplyr) # 操作数据的包
> library(ggplot2) # 绘图包
> library(ggthemes) # ggplot2的主题修改包
> library(scales) # 可视化的包
> library(VIM) # 查看缺失数据的包
> library(mice) # 插补数据的包
> library(randomForest) # 随机森林
导入并初步观察数据：
> train <- read.csv(file.choose(),stringsAsFactors = FALSE)
> test <- read.csv(file.choose(),stringsAsFactors = FALSE)
> str(train)
略
> str(test)
略
两个数据集除Survived字段不同外，其他字段均相同。为了方便数据清洗，我们合并训练集与测试集。
> full <- bind_rows(train, test)
> str(full)
'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(full)
PassengerId Survived Pclass Name Sex Age
Min. : 1 Min. :0.0000 Min. :1.000 Length:1309 Length:1309 Min. : 0.17
1st Qu.: 328 1st Qu.:0.0000 1st Qu.:2.000 Class :character Class :character 1st Qu.:21.00
Median : 655 Median :0.0000 Median :3.000 Mode :character Mode :character Median :28.00
Mean : 655 Mean :0.3838 Mean :2.295 Mean :29.88
3rd Qu.: 982 3rd Qu.:1.0000 3rd Qu.:3.000 3rd Qu.:39.00
Max. :1309 Max. :1.0000 Max. :3.000 Max. :80.00
NA's :418 NA's :263
SibSp Parch Ticket Fare Cabin Embarked
Min. :0.0000 Min. :0.000 Length:1309 Min. : 0.000 Length:1309 Length:1309
1st Qu.:0.0000 1st Qu.:0.000 Class :character 1st Qu.: 7.896 Class :character Class :character
Median :0.0000 Median :0.000 Mode :character Median : 14.454 Mode :character Mode :character
Mean :0.4989 Mean :0.385 Mean : 33.295
3rd Qu.:1.0000 3rd Qu.:0.000 3rd Qu.: 31.275
Max. :8.0000 Max. :9.000 Max. :512.329
NA's :1
合并后的数据中，共有1309个观测，其中训练集891个，测试集418个，生存情况（Survived）中缺失值418个（正常，需要预测的部分），年龄（Age）中缺失值263个，船票费用（Fare）中缺失值1个。
解释一下各个变量对应的含义：
3.数据清洗
3.1 观察姓名变量
我们注意到乘客名字（Name）有一个非常显著的特点：每个名字当中都包含了具体的称谓或者头衔。我们可以将这部分信息提取出来，或许可以作为非常有用的新变量。
> # 乘客姓名提取头衔
> full$Title <- gsub("(.*, )|(\\..*)", "", full$Name)
> # 按照性别划分头衔数量
> table(full$Sex, full$Title)
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
> # 对于那些出现频率较低的头衔合并为一类
> rare_title <- c("Capt","Col","Don","Dona","Dr", "Jonkheer", "Lady","Major", "Rev", "Sir","the Countess")
> # 数据量不多，我想试试手工调整头衔
> filter(full,full$Title %in% rare_title)%>%
+ select(Name,Sex,Age,SibSp,Parch,Title)%>%
+ arrange(Title)
Name Sex Age SibSp Parch Title
1 Crosby, Capt. Edward Gifford male 70 1 1 Capt
2 Simonius-Blumer, Col. Oberst Alfons male 56 0 0 Col
3 Weir, Col. John male 60 0 0 Col
4 Gracie, Col. Archibald IV male 53 0 0 Col
5 Astor, Col. John Jacob male 47 1 0 Col
6 Uruchurtu, Don. Manuel E male 40 0 0 Don
7 Oliva y Ocana, Dona. Fermina female 39 0 0 Dona
8 Minahan, Dr. William Edward male 44 2 0 Dr
9 Moraweck, Dr. Ernest male 54 0 0 Dr
10 Pain, Dr. Alfred male 23 0 0 Dr
11 Stahelin-Maeglin, Dr. Max male 32 0 0 Dr
12 Frauenthal, Dr. Henry William male 50 2 0 Dr
13 Brewe, Dr. Arthur Jackson male NA 0 0 Dr
14 Leader, Dr. Alice (Farnham) female 49 0 0 Dr
15 Dodge, Dr. Washington male 53 1 1 Dr
16 Reuchlin, Jonkheer. John George male 38 0 0 Jonkheer
17 Duff Gordon, Lady. (Lucille Christiana Sutherland) ("Mrs Morgan") female 48 1 0 Lady
18 Peuchen, Major. Arthur Godfrey male 52 0 0 Major
19 Butt, Major. Archibald Willingham male 45 0 0 Major
20 Byles, Rev. Thomas Roussel Davids male 42 0 0 Rev
21 Bateman, Rev. Robert James male 51 0 0 Rev
22 Carter, Rev. Ernest Courtenay male 54 1 0 Rev
23 Kirkland, Rev. Charles Leonard male 57 0 0 Rev
24 Harper, Rev. John male 28 0 1 Rev
25 Montvila, Rev. Juozas male 27 0 0 Rev
26 Lahtinen, Rev. William male 30 1 1 Rev
27 Peruschitz, Rev. Joseph Maria male 41 0 0 Rev
28 Duff Gordon, Sir. Cosmo Edmund ("Mr Morgan") male 49 1 0 Sir
29 Rothes, the Countess. of (Lucy Noel Martha Dyer-Edwards) female 33 0 0 the Countess
> # 我将Lady调整为Mrs，Sir调整为Mr
> rare_title2 <- c("Capt","Col","Don","Dona","Dr", "Jonkheer","Major", "Rev","the Countess")
> # 对于一些称呼进行重新指定（按含义） 如mlle, ms指小姐, mme 指女士
> full$Title[full$Title %in% c("Mlle","Ms")]<- "Miss"
> full$Title[full$Title %in% c("Mme","Lady")]<- "Mrs"
> full$Title[full$Title =="Sir"]<- "Mr"
> full$Title[full$Title %in% rare_title2] <- "Rare Title"
> # 查看重新调整后的情况
> table(full$Sex, full$Title)
Master Miss Mr Mrs Rare Title
female 0 264 0 199 3
male 61 0 758 0 24
> # 从乘客姓名中，提取姓氏,作为家庭变量
> full$Surname <- sapply(full$Name, function(x) strsplit(x, split = '[,.]')[[1]][1])
> length(unique(full$Surname))
[1] 875
Megan L. Risdal 的文章中，将乘客姓氏也提取出来，提示可以发掘乘客姓氏之间的联系，但没有进行进一步操，我们这里也就不探讨了。
3.2 观察家庭情况
处理完乘客姓名这一变量，我们再考虑衍生一些家庭相关的变量。基于已有变量SubSp和Parch生成家庭人数family size 这一变量。
> # 生成家庭人数变量，包括自己在内
> full$Fsize <- full$SibSp + full$Parch + 1
> # 生成一个家庭变量：格式为姓_家庭人数
> full$Family <- paste(full$Surname, full$Fsize, sep="_")
> # 使用 ggplot2 绘制家庭人数与生存情况之间的关系
> ggplot(full[1:891,], aes(x = Fsize, fill = factor(Survived))) +
+ geom_bar(stat="count", position="dodge") +
+ scale_x_continuous(breaks=c(1:11)) +
+ labs(x = "Family Size") +
+ theme_few()
通过图形我们可以明显发现以下特点：
- 个人和家庭人数>4的家庭，存活人数小于死亡人数
- 家庭人数在[2:4]的家庭，存活人数大于死亡人数
因此，我们可以将家庭人数变量进行分段合并，明显的可以分为3段：个人，小家庭，大家庭，由此生成新变量。
> # 离散化
> full$FsizeD[full$Fsize == 1] <- "single"
> full$FsizeD[full$Fsize < 5 & full$Fsize > 1]<- "small"
> full$FsizeD[full$Fsize > 4] <- "large"
> # 通过马赛克图查看家庭规模与生存情况之间的关系
> mosaicplot(table(full$FsizeD,full$Survived), main="Family Size by Survival", shade=TRUE)
显而易见，个人与大家庭不利于在沉船事故中生存，而小家庭当中生存率相对较高。
3.3 尝试生成更多变量
在乘客客舱变量Cabin中，也存在一些有价值的信息，如客舱层数deck。
> # 可以看出这一变量有很多缺失值,有单个客户多个客舱，格式基本为字母+数字
> head(full$Cabin,30)
[1] "" "C85" "" "C123" "" "" "E46"
[8] "" "" "" "G6" "C103" "" ""
[15] "" "" "" "" "" "" ""
[22] "D56" "" "A6" "" "" "" "C23 C25 C27"
[29] "" ""
> # 假设第一个字母即为客舱层数，建立一个层数变量（Deck），取值范围A - z:
> full$Deck<-factor(sapply(full$Cabin, function(x) strsplit(x, NULL)[[1]][1]))
> summary(full$Deck)
A B C D E F G T NA's
22 65 94 46 41 21 5 1 1014
上面只是对变量进行了基本处理，还有很多可以进一步操作的地方，如有些乘客名下包含很多间房 (如28行， "C23 C25 C27"), 但是这一变量有1014 个缺失值，占比太高了。 后面就不再进一步考虑。
4.处理缺失值
我们先找到所有的缺失数据看看情况：
- 数值型数据：age缺失263个，fare缺失1个
- 字符型数据：cabin缺失1014个，embarked缺失2个
- 因子型数据：deck缺失1014个
> sapply(full,function(x) {sum(is.na(x))})
PassengerId Survived Pclass Name Sex Age SibSp Parch Ticket
0 418 0 0 0 263 0 0 0
Fare Cabin Embarked Title Surname Fsize Family FsizeD Deck
1 0 0 0 0 0 0 0 1014
> sapply(full,function(x) {sum(x == "",na.rm=TRUE)})
PassengerId Survived Pclass Name Sex Age SibSp Parch Ticket
0 0 0 0 0 0 0 0 0
Fare Cabin Embarked Title Surname Fsize Family FsizeD Deck
0 1014 2 0 0 0 0 0 0
处理缺失值的方法有很多种，考虑到数据集本身较小，样本数也不多，因而不能直接整行或者整列删除缺失值样本。我们考虑对于缺失值较少的，用均值或中位数填补，缺失值较多的通过现有数据和变量进行预估填补。
4.1 登船港口缺失
找到缺失数据在哪一行，观察情况。
> filter(full,full$Embarked=="")
PassengerId Survived Pclass Name Sex Age SibSp Parch Ticket Fare
1 62 1 1 Icard, Miss. Amelie female 38 0 0 113572 80
2 830 1 1 Stone, Mrs. George Nelson (Martha Evelyn) female 62 0 0 113572 80
Cabin Embarked Title Surname Fsize Family FsizeD Deck
1 B28 Miss Icard 1 Icard_1 single B
2 B28 Mrs Stone 1 Stone_1 single B
我们可以看到他们支付的票价都是$ 80，舱位等级都是1，我们可以大胆推论有相同舱位等级（passenger class）和票价（Fare）的乘客，也许有着相同的登船港口？
> # 去除缺失值乘客的ID
> embark_fare <- full %>% filter(PassengerId != 62 & PassengerId != 830)
> # 用 ggplot2 绘制embarkment, passenger class, & median fare 三者关系图
> ggplot(embark_fare, aes(x = Embarked, y = Fare, fill = factor(Pclass))) +
+ geom_boxplot() +
+ geom_hline(aes(yintercept=80),
+ colour="red", linetype="dashed", lwd=2) +
+ scale_y_continuous(labels=dollar_format()) +
+ theme_few()
很明显，港口C出发的头等舱支付的票价中位数为80，因此我们可以把两个缺失值替换为C。
> full$Embarked[c(62, 830)] <- "C"
4.1 票价缺失
找到缺失数据在哪一行，观察情况。
> filter(full,is.na(full$Fare))
PassengerId Survived Pclass Name Sex Age SibSp Parch Ticket Fare Cabin Embarked Title Surname
1 1044 NA 3 Storey, Mr. Thomas male 60.5 0 0 3701 NA S Mr Storey
Fsize Family FsizeD Deck
1 1 Storey_1 single <NA>
这是从港口S出发的三等舱乘客，根据上面登船港口缺失值填补的推论，我们可以作图看看：
> ggplot(full[full$Pclass == "3" & full$Embarked == "S", ],
+ aes(x = Fare)) +
+ geom_density(fill = "#99d6ff", alpha=0.4) +
+ geom_vline(aes(xintercept=median(Fare, na.rm=T)),colour="red", linetype="dashed", lwd=1) +
+ scale_x_continuous(labels=dollar_format()) +
+ theme_few()
从得到的图形上看,将缺失值用中位数进行替换是合理的。替换数值为$8.05。
> full$Fare[1044] <- median(full[full$Pclass==3&full$Embarked=="S",]$Fare,na.rm=TRUE)
> full$Fare[1044]
[1] 8.05
4.3 年龄缺失
用户年龄（Age） 中有大量缺失存在，简单用中位数或均值肯定不妥，这里我们用mice包的随机插补，将基于其他变量构建一个预测模型对年龄缺失值进行填补。
> # 使变量因子化
> factor_vars <- c("PassengerId","Pclass","Sex","Embarked",
+ "Title","Surname","Family","FsizeD")
> full[factor_vars] <- lapply(full[factor_vars],function(x) as.factor(x))
> # 设置随机种子
> set.seed(123)
> # 执行多重插补法，剔除一些没什么用的变量
> mice_mod <- mice(full[, !names(full) %in%
+ c("PassengerId","Name","Ticket","Cabin",
+ "Family","Surname","Survived")],
+ method="rf")
iter imp variable
1 1 Age Deck
1 2 Age Deck
1 3 Age Deck
1 4 Age Deck
1 5 Age Deck
2 1 Age Deck
2 2 Age Deck
2 3 Age Deck
2 4 Age Deck
2 5 Age Deck
3 1 Age Deck
3 2 Age Deck
3 3 Age Deck
3 4 Age Deck
3 5 Age Deck
4 1 Age Deck
4 2 Age Deck
4 3 Age Deck
4 4 Age Deck
4 5 Age Deck
5 1 Age Deck
5 2 Age Deck
5 3 Age Deck
5 4 Age Deck
5 5 Age Deck
> # 保存完整输出
> mice_output <- complete(mice_mod)
> # 绘制年龄分布图
> par(mfrow=c(1,2))
> hist(full$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))
从上面的图来看，数据填补前后，并没有发生明显的偏移，随机插补应该是有效的，那么下面可以用mice模型的结果对原年龄数据进行替换。
> full$Age <- mice_output$Age
> # 检查缺失值是否被完全替换了
> sum(is.na(full$Age))
[1] 0
5.特征工程
现在我们知道每一位乘客的年龄，那么基于“妇女与儿童优先”的原则，我们可以生成一些变量，如儿童（Child）和 母亲（Mother）。
划分标准：
- 儿童 ： 年龄Age < 18
- 母亲 ： 女性+年龄 > 18+拥有超过1个子女+头衔不是'Miss'
> # 首先我们来看年龄与生存情况之间的关系
> ggplot(full[1:891,], aes(Age, fill = factor(Survived))) +
+ geom_histogram() +
+ # 分性别来看，因为我们知道性别对于生存情况有重要影响
+ facet_grid(.~Sex) +
+ theme_few()
生成child变量, 并且基于此划分儿童child与成人adult。
> full$child[full$Age < 18] <- "Child"
> full$child[full$Age >= 18] <- "Adult"
> # 展示对应人数
> table(full$child, full$Survived)
0 1
Adult 479 270
Child 70 72
下面来生成母亲这个变量。
> #增加mother变量
> full$Mother <- "NotMother"
> full$Mother[full$Sex == "female" & full$Parch > 0 & full$Age > 18 & full$Title != "Miss"] <- "Mother"
> table(full$Mother, full$Survived)
0 1
Mother 15 39
NotMother 534 303
将child和mother变量转化成因子类型
> full$child <- factor(full$child)
> full$Mother <- factor(full$Mother)
至此，所有我们需要的变量都已经生成，并且其中没有缺失值。 为了保险起见，我们进行二次确认。
> sapply(full,function(x){sum(is.na(x))})
PassengerId Survived Pclass Name Sex Age SibSp Parch Ticket
0 418 0 0 0 0 0 0 0
Fare Cabin Embarked Title Surname Fsize Family FsizeD Deck
0 0 0 0 0 0 0 0 1014
child Mother
0 0
> sapply(full,function(x){sum(x=="",na.rm=T)})
PassengerId Survived Pclass Name Sex Age SibSp Parch Ticket
0 0 0 0 0 0 0 0 0
Fare Cabin Embarked Title Surname Fsize Family FsizeD Deck
0 1014 0 0 0 0 0 0 0
child Mother
0 0
6.模型设定与预测
在完成上面的工作之后，我们进入到最后一步：预测泰坦尼克号上乘客的生存状况。 在这里我们使用随机森林分类算法(The RandomForest Classification Algorithm) 。
第一步，拆分训练集与测试集
> train <- full[1:891,]
> test <- full[892:1309,]
第二步， 建立模型
> # 设置随机种子
> set.seed(1234)
> # 建立模型 (注意: 不是所有可用变量全部加入)
> rf_model <- randomForest(factor(Survived) ~ Pclass + Sex + Age + SibSp +Parch +
+ Fare + Embarked + Title + FsizeD + child + Mother,
+ data = train)
> # 显示模型误差
> plot(rf_model, ylim=c(0,0.36))
> legend("topright", colnames(rf_model$err.rate), col=1:3, fill=1:3)
黑色那条线表示：整体误差率（the overall error rate）保持在20% 左右
红色和绿色线分别表示：遇难与生还的误差率，红线保持在10%，绿线保持在30%左右。
我们的模型，在预测死亡情况时更准确一些。
第三步， 查看变量重要性
> importance <- importance(rf_model)
> varImportance <- data.frame(Variables = row.names(importance), Importance = round(importance[ ,"MeanDecreaseGini"],2))
> rankImportance <- varImportance %>%
+ mutate(Rank = paste0("#",dense_rank(desc(Importance))))
> # 作图
> ggplot(rankImportance, aes(x = reorder(Variables, Importance), y = Importance, fill = 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()
我们从图上可以看出，头衔和票价对于生存情况影响最大，其次是性别和年龄，而乘客舱位排第五。 而母亲和孩子对于生存与否的影响最小，真是有点出乎意料。
第四步， 预测
> # 将模型带入测试集
> prediction <- predict(rf_model, test)
> # 保存结果
> solution <- data.frame(PassengerID = test$PassengerId, Survived = prediction)
> # 输出结果到CSV文件格式
> write.csv(solution, file = "rf_mod_Solution.csv", row.names = F)
7.总结
本次数据分析是基于Megan L. Risdal 的文章，算是一次深度模仿，而且对于其中的一些操作还有不太理解的部分，还应该多看多学，争取早日独立完成一项数据分析工作。
