题外话：在文章正式开始之前，我还是想先写一点题外话，一是为了引出写作这篇博客的目的，二则是希望能够记录下现在的所思所想为以后留个纪念。首先介绍一下自己，毕业3年多的小硕一枚，大学期间学的专业是高分子材料，毕业后也一直在从事化工行业方面的工作。最近由于公司变动，再加上本身自己对于未来规划以及个人兴趣的原因，使我产生了转行的想法。而数据分析作为我在学生时代就已经有所接触的方向，确实对我有着很大的吸引，所以在近半年的时间里我开始了对数据分析的学习。期间从统计学原理开始，到最常用的软件Excel，再到SPSS、R语言、MySQL等数据分析和数据库软件，最后还简单的学习了如何使用python做网络爬虫以及基于python的Scrapy框架。下面我将尝试利用R语言对Kaggle上的泰坦尼克号项目进行分析，作为初入数据分析的入门者，希望能够通过这个项目对自己的学习进行一次检验。
利用R语言对泰坦尼克号沉没事件幸存者的数据分析与预测
1. 背景简介
泰坦尼克号(RMS Titanic)作为英国白星航运公司下的一艘奥林匹克级邮轮，是20世纪最著名的邮轮之一，这不仅因为其在建成后是当时世界上最大的豪华客运轮船，还因为这艘被誉为‘永不沉没’的轮船竟然在其处女航中就遭遇了厄运。在1912年4月15号，载有2224人的泰坦尼克号因与冰山相撞而沉没，这次事故造成了1502人死亡，同时也震惊了世界。在这次事故中，造成如此多的伤亡主要原因之一是由于船上没有配备足够的救生船，但是对于幸存者而言还有一些其他的因素导致了他们更有可能比别人活下来。在这里本文将对其他可能的影响因素进行分析，寻找出那些人在泰坦尼克号上更有可能成为幸存者。
该项目的来源为Kaggle上的一篇名为“Titanic: Machine Learning from Disaster”的Competitions，文中所有资料与数据均来源于此。
2. 数据的导入与整理
在这里Kaggle提供了三组数据分别为train、test和gender_submission，均为csv文件。在这里可以使用R语言直接读取，不过我还是将数据先导入到了MySQL中，然后再用R来读取，主要是为了在读取和整理数据时方便一些。导入MySQL的方法就不在这里介绍了，下面是用R读取数据的语句。
library(DBI)
library(RMySQL)
conn <- dbConnect(RMySQL::MySQL(), dbname = "titanic", username = "root",password = "root", port = 3306)
train <- dbSendQuery(conn, "SELECT Passenger, Survival, Pclass, Sex, Age, Sibsp, Parch, fare, Embarked FROM train")
test <- dbSendQuery(conn, "SELECT test.Passenger, Survival, Pclass, Sex, Age, Sibsp, Parch, fare, Embarked FROM test, gender_submission WHERE test.Passenger = gender_submission.Passenger")
traindata <- dbFetch(train, n = -1)
testdata <- dbFetch(test, n = -1)> str(traindata)
'data.frame': 891 obs. of 9 variables:
$ Passenger: int 1 2 3 4 5 6 7 8 9 10 ...
$ Survival : int 0 1 1 1 0 0 0 0 1 1 ...
$ Pclass : int 3 1 3 1 3 3 1 3 3 2 ...
$ Sex : chr "male" "female" "female" "female" ...
$ Age : chr "22" "38" "26" "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 ...
$ fare : num 7.25 71.28 7.92 53.1 8.05 ...
$ Embarked : chr "S" "C" "S" "S" ... > str(testdata)
'data.frame': 417 obs. of 9 variables:
$ Passenger: int 892 893 894 895 896 897 898 899 900 901 ...
$ Survival : int 0 1 0 0 1 0 1 0 1 0 ...
$ Pclass : int 3 3 2 3 3 3 3 2 3 3 ...
$ Sex : chr "male" "female" "male" "male" ...
$ Age : chr "34.5" "47" "62" "27" ...
$ Sibsp : int 0 1 0 0 1 0 0 1 0 2 ...
$ Parch : int 0 0 0 0 1 0 0 1 0 0 ...
$ fare : num 7.83 7 9.69 8.66 12.29 ...
$ Embarked : chr "Q" "S" "Q" "S" ...
其中需要注意的是在Age和Embarked中或多或少的都有一部分空值，而如果从MySQL直接读取到R中这部分空值会形成list格式，使用na.omit()或complete.cases()都不能去除空值，所以在导入R之前，需要现在MySQL中将这些没有数据的位置都替换为NULL。
在数据表train中，Kaggle给出了12组变量，其中Passenger作为编号可作为数据主键没有实际意义。Name虽有一定规律，但其所表达的家庭联系可以由Parch和Sibsp来替代，并且姓名的分析复杂并且可靠性较低，所以在这里不对姓名进行分析；而Ticket的编号随机性较大，且没有规律，所以也不做分析；Cabin的空值过多，无法代表整体数据，所以也不做分析。下面对剩下的7组变量与Survival之间的关系进行分析。
3. 什么因素决定着船上人的生与死？
通过上一部分的数据整理后，剩下的七种因素可能影响着一个人是否成为幸存者，而这七种因素到底是如何影响着一个人的命运，又是那种因素的作用更大呢，下面对各个因素进行逐一分析。
3.1 头等舱的重要性
首先我们来看乘客所坐的船舱是如何影响生还率的。在train中的891个样本中全部记录了乘客所乘坐的船舱等级，其中船舱被分为三个等级，具体分析如下：
> library(ggplot2)
> Pclass_S <- table(traindata$Survival, traindata$Pclass)
> Pclass_S_prop <- prop.table(Pclass_S, 2)
> ggplot(data = traindata, aes(x = Pclass, fill = factor(Survival)))+geom_bar(stat='count', position='dodge') + scale_x_continuous(breaks=c(1:3)) + labs(x = 'Pclass')> Pclass_S
1 2 3
0 80 97 372
1 136 87 119
> Pclass_S_prop
1 2 3
0 0.3703704 0.5271739 0.7576375
1 0.6296296 0.4728261 0.2423625
由上图可以看出，在三个不同等级的船舱中，3等舱中人数最多为491人，1等和2等舱人数相近，分别为216人和184人。然而在生还率上1等舱为62.96%，2等舱为47.28%，3等舱为24.24%。可以明显的看出随着船舱等级的下降乘客的生还率呈现出显著的下降趋势，也就是说在泰坦尼克号事件中，乘坐1等舱可以有效的帮助其成为幸存者。
3.2 女士优先
进入现代社会后，女士优先已经不再是男士们表现绅士风度的行为，女士优先已经成为了一种类似于道德标准的社会行为准则，而在近百年前的泰坦尼克号上，人们是否也会遵循女士优先的准则呢？
> Sex_S <- table(traindata$Survival, traindata$Sex)
> Sex_S_prop <- prop.table(Sex_S, 2)
> ggplot(data = traindata, aes(x = Sex, fill = factor(Survival)))+geom_bar(stat='count', position='dodge')> Sex_S
female male
0 81 468
1 233 109
> Sex_S_prop
female male
0 0.2579618 0.8110919
1 0.7420382 0.1889081
从以上数据可以看出，在train中所记录的女性有314人，而成功生还的有233人，生还率为74.20%，而相对的男性的生还率则只有18.89%。在面对重大灾难时，男士们表现出了高尚的品质，把更多的生的希望给了女性。所以身为女性可以有效的提高成为幸存者的概率。
3.3 年龄的影响
下面我们来看一下年龄对生还率的影响，是否身强体壮有利于成为幸存者。
> Agedata <- as.numeric(unlist(traindata$Age))
> Age_S <- table(traindata$Survival, cut(Agedata, breaks = c(0, 15, 30, 45, 60, 75, 90), labels = c('kids', 'teenagers', 'prime', 'middle', 'agedness', 'senium' )))
> Age_S_prop <- prop.table(Age_S, 2)
> ggplot(data = data.frame(traindata$Survival, Agedata), aes(x = cut(Agedata, breaks = c(0, 15, 30, 45, 60, 75, 90)), fill = factor(traindata.Survival)))+geom_bar(stat='count', position='dodge') + labs(x = 'Age') + scale_x_discrete(labels = c('kids', 'teenagers', 'prime', 'middle', 'agedness', 'senium'))> Age_S
kids teenagers prime middle agedness senium
0 34 209 116 48 17 0
1 49 117 86 33 4 1
> Age_S_prop
kids teenagers prime middle agedness senium
0 0.4096386 0.6411043 0.5742574 0.5925926 0.8095238 0.0000000
1 0.5903614 0.3588957 0.4257426 0.4074074 0.1904762 1.0000000
在这里将年龄以15岁一个节点分为六个区间，并根据年龄大小对每个区间进行了分类描述，在排除75岁到90岁的一个个例后，在各区间中生还率最高的是年龄小于15岁的未成年人，生还率达到了59%；而生还率最低的则为年龄在60岁到75岁之间的老年人，其生还率为19.05%。对于年龄在15岁到60岁之间的成年人来说，虽然年龄跨度较大，但在生还率上并没有显著差异，最高的仅比最低的高出大约6.5%左右。所以从对年龄的分析可以看出，在面对重大灾难时，未成年人是人们首先保护的对象，其最有可能成为幸存者。
3.4 亲情的好与坏
在这一节我们对Sibsp和Parch两个变量进行讨论，Sibsp指的是一个人拥有的旁系亲属的数量，而Parch则代表直系亲属的数量。由于都是比较亲密的关系，并且这种关系在面对灾难时会对人们的生还造成影响，所以在这里对这两个变量一同进行讨论。
Sibsp_S <- table(traindata$Survival, traindata$Sibsp)
Parch_S <- table(traindata$Survival, traindata$Parch)
Sibsp_S_prop <- prop.table(Sibsp_S, 2)
Parch_S_prop <- prop.table(Parch_S, 2)
ggplot(data = traindata, aes(x = Sibsp, fill = factor(Survival)))+geom_bar(stat='count', position='dodge') + scale_x_continuous(breaks=c(0:8)) + labs(x = 'Sibsp')
ggplot(data = traindata, aes(x = Parch, fill = factor(Survival)))+geom_bar(stat='count', position='dodge') + scale_x_continuous(breaks=c(0:6)) + labs(x = 'Parch')
Families <- traindata$Sibsp +traindata$Parch
ggplot(data = traindata, aes(x = Families, fill = factor(Survival)))+geom_bar(stat='count', position='dodge') + scale_x_continuous(breaks=c(0:10)) + labs(x = 'Families')> Sibsp_S
0 1 2 3 4 5 8
0 398 97 15 12 15 5 7
1 210 112 13 4 3 0 0
> Sibsp_S_prop
0 1 2 3 4 5 8
0 0.6546053 0.4641148 0.5357143 0.7500000 0.8333333 1.0000000 1.0000000
1 0.3453947 0.5358852 0.4642857 0.2500000 0.1666667 0.0000000 0.0000000
> Parch_S
0 1 2 3 4 5 6
0 445 53 40 2 4 4 1
1 233 65 40 3 0 1 0
> Parch_S_prop
0 1 2 3 4 5 6
0 0.6563422 0.4491525 0.5000000 0.4000000 1.0000000 0.8000000 1.0000000
1 0.3436578 0.5508475 0.5000000 0.6000000 0.0000000 0.2000000 0.0000000
3.5 现金和港口
一个乘客上船时所带的现金，以及他所登船的港口会对他成为幸存者有影响么？这两个看似和成为幸存者毫无关系的因素，可能正从侧面表现出了幸存者所拥有的属性。那么还是首先从简单的单因素统计绘图开始。其中将Fare这一变量分为三个区间，第一个区间为(0, 50]标签为poor，第二个区间为(50, 100]标签为middle，第三个区间为(100, 600]标签为rich。
Faredata <- as.numeric(unlist(traindata$fare))
Fare_S <- table(traindata$Survival, cut(Faredata, breaks = c(0, 50, 100, 600), labels = c('poor', 'middle', 'rich')))
Fare_S_prop <- prop.table(Fare_S, 2)
ggplot(data = data.frame(traindata$Survival, Faredata), aes(x = cut(Faredata, breaks = c(0, 50, 100, 600)), fill = factor(traindata.Survival)))+geom_bar(stat='count', position='dodge') + labs(x = 'Fare') + scale_x_discrete(labels = c('poor', 'middle', 'rich'))
Embarked_S <- table(traindata$Survival, traindata$Embarked)
Embarked_S_prop <- prop.table(Embarked_S, 2)
ggplot(data = traindata, aes(x = Embarked, fill = factor(Survival)))+geom_bar(stat='count', position='dodge')> Fare_S
poor middle rich
0 484 37 14
1 232 70 39
> Fare_S_prop
poor middle rich
0 0.6759777 0.3457944 0.2641509
1 0.3240223 0.6542056 0.7358491
> Embarked_S
C Q S
0 75 47 427
1 93 30 217
> Embarked_S_prop
C Q S
0 0.4464286 0.6103896 0.6630435
1 0.5535714 0.3896104 0.3369565
由以上 数据可以看出，对于Fare这个变量而言，登船时所带的现金越多成为幸存者的概率也就越高。而对于登船港口而言， 从Cherbourg-Octeville登船的乘客要比从Southampton或Queenstown登船的乘客活下来的概率高。那么下面先对Fare这个变量进行分析，看看带有大量现金的乘客他们的其他属性如何？
plot_FP <- ggplot(data = data.frame(traindata$Pclass, Faredata), aes(x = cut(Faredata, breaks = c(0, 50, 100, 600)), fill = factor(traindata.Pclass)))+geom_bar(stat='count', position='dodge') + labs(x = 'Fare') + scale_x_discrete(labels = c('poor', 'middle', 'rich')) + scale_fill_brewer(palette = "Blues") + theme_bw()
plot_FS <- ggplot(data = data.frame(traindata$Sex, Faredata), aes(x = cut(Faredata, breaks = c(0, 50, 100, 600)), fill = factor(traindata.Sex)))+geom_bar(stat='count', position='dodge') + labs(x = 'Fare') + scale_x_discrete(labels = c('poor', 'middle', 'rich')) + scale_fill_brewer(palette = "Blues") + theme_bw()
plot_FA <- ggplot(data = data.frame(Agedata, Faredata), aes(x = cut(Faredata, breaks = c(0, 50, 100, 600)), fill = factor(Age_Fare)))+geom_bar(stat='count', position='dodge') + labs(x = 'Fare') + scale_x_discrete(labels = c('poor', 'middle', 'rich')) + scale_fill_brewer(palette = "Blues") + theme_bw()
plot_FF <- ggplot(data = data.frame(Families, Faredata), aes(x = cut(Faredata, breaks = c(0, 50, 100, 600)), fill = factor(Families)))+geom_bar(stat='count', position='dodge') + labs(x = 'Fare') + scale_x_discrete(labels = c('poor', 'middle', 'rich')) + scale_fill_brewer(palette = "Blues") + theme_bw()library(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))
vplayout = function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(plot_FP, vp = vplayout(1, 1))
print(plot_FS, vp = vplayout(1, 2))
print(plot_FA, vp = vplayout(2, 1))
print(plot_FF, vp = vplayout(2, 2))
从上图可以看出，标签为rich和middle的乘客，他们普遍乘坐的是1等舱，并且这部分乘客中女性的比例要高于男性，也就是说带有大量现金的都是一些乘坐1等舱的女性乘客。而根据之前的分析，乘坐1等舱和身为女性都会极大的提升成为幸存者的概率，所以在这里可以理解为携带大量现金并不能帮助乘客成为幸存者，而是活下来的幸存者中有很大一部分恰好喜欢携带大量现金而已。同理可对登船港口Embarked进行如下分析：
plot_EP <- ggplot(data = traindata, aes(x = (Embarked), fill = factor(Pclass)))+geom_bar(stat='count', position='dodge') + labs(x = 'Embarked') + scale_x_discrete(labels = c('Cherbourg-Octeville', 'Queenstown', 'Southampton')) + scale_fill_brewer(palette = "Blues") + theme_bw()
plot_ES <- ggplot(data = traindata, aes(x = (Embarked), fill = factor(Sex)))+geom_bar(stat='count', position='dodge') + labs(x = 'Embarked') + scale_x_discrete(labels = c('Cherbourg-Octeville', 'Queenstown', 'Southampton')) + scale_fill_brewer(palette = "Blues") + theme_bw()
plot_EA <- ggplot(data = data.frame(Agedata, traindata$Embarked), aes(x = (traindata.Embarked), fill = factor(Age_cut)))+geom_bar(stat='count', position='dodge') + labs(x = 'Embarked') + scale_x_discrete(labels = c('Cherbourg-Octeville', 'Queenstown', 'Southampton')) + scale_fill_brewer(palette = "Blues") + theme_bw()
plot_EF <- ggplot(data = data.frame(Families, traindata$Embarked), aes(x = (traindata.Embarked), fill = factor(Families)))+geom_bar(stat='count', position='dodge') + labs(x = 'Embarked') + scale_x_discrete(labels = c('Cherbourg-Octeville', 'Queenstown', 'Southampton')) + scale_fill_brewer(palette = "Blues") + theme_bw()
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))
vplayout = function(x, y) viewport(layout.pos.row = x, layout.pos.col = y)
print(plot_EP, vp = vplayout(1, 1))
print(plot_ES, vp = vplayout(1, 2))
print(plot_EA, vp = vplayout(2, 1))
print(plot_EF, vp = vplayout(2, 2))
可以看出从 Cherbourg-Octeville港登船的乘客多为乘坐1等舱的女性，所以当进行统计分析时，从Cherbourg-Octeville港登船的乘客成为幸存者的概率要远远高于其他两港的乘客。
4. 基于分类算法的幸存者预测
由于该项目希望能够通过数据进行机器学习，并得到可以用于预测的模型，所以本文中尝试采用三种分类算法模型对所给出的数据进行学习，采用的三种算法分别为K临近法(K-mean)、随机森林法(Random Forest)和支持向量机法(SVM)。下面对着三种算法进行逐一讨论。
4.1 K临近法则的模型训练及预测
在对模型的训练过程中，所选用的训练数据为Kaggle给出的train数据，用于预测的数据为test数据。由于在数据矩阵中存在大量空值，所以在训练前应先将具有空值的对象剔除，并且由于我在运算时，输入的数据类型为数据框格式，但函数在计算时自动转换为矩阵的过程中依然会出现"Warning in data.matrix(x): NAs introduced by coercion"这样的警告，所以在进行训练之前对训练数据和测试数据均进行了手动转换为矩阵格式，下面为进行KNN训练的代码：
library(class)
result <- M_testdata[,1]
accuracy <- vector()
for(i in 1:30){
KNN_M <- knn(train = M_traindata[, 2:8], test = M_testdata[, 2:8], cl = M_traindata[, 1], k = i)
CT <- table(result, KNN_M)
accuracy <- c(accuracy, sum(diag(CT))/sum(CT)*100)
}
> accuracy
[1] 57.70393 60.12085 59.81873 59.21450 62.83988 61.32931 62.53776 63.14199 62.53776 61.02719 60.72508
[12] 59.81873 57.70393 57.70393 58.61027 60.42296 60.42296 60.72508 61.93353 61.02719 61.93353 62.23565
[23] 62.53776 62.53776 61.93353 62.53776 61.93353 60.72508 60.42296 61.32931set.seed(123)
random <- sample(1:1043, 104)
result_KNN_R <- Random_data[random, 1]
accuracy_KNN_R <- vector()
for(i in 1:30){
KNN_M <- knn(train = Random_data[-random, 2:8], test = Random_data[random, 2:8], cl = Random_data[-random, 1], k = i)
CT_KNN_R <- table(result_KNN_R, KNN_M)
accuracy_KNN_R <- c(accuracy_KNN_R, sum(diag(CT_KNN_R))/sum(CT_KNN_R)*100)
}
accuracy_KNN_R_Max <- max(accuracy_KNN_R)> accuracy_KNN_R
[1] 69.23077 70.19231 71.15385 67.30769 71.15385 75.00000 71.15385 70.19231 70.19231 70.19231 68.26923
[12] 66.34615 68.26923 64.42308 62.50000 62.50000 64.42308 64.42308 66.34615 65.38462 65.38462 66.34615
[23] 67.30769 66.34615 67.30769 66.34615 66.34615 66.34615 66.34615 66.346154.2 随即森林的模型训练及预测
在上一节中train和test数据被合并，并重新通过随机抽取的方式重新分组了训练数据和测试数据，所以在这一节以及之后的训练中都采用重新分组的样本数据。对样本数据进行随即森林模型训练，并使用测试数据进行预测。
> set.seed(123)
> RF <- randomForest(factor(Survival_train) ~ ., data = Random_data[-random,], importance = TRUE)> RF
Call:
randomForest(formula = factor(Survival_train) ~ ., data = Random_data[-random, ], importance = TRUE)
Type of random forest: classification
Number of trees: 500
No. of variables tried at each split: 2
OOB estimate of error rate: 14.59%
Confusion matrix:
0 1 class.error
0 515 52 0.09171076
1 85 287 0.22849462
RF_tree <- plot(RF)
tree <- c(1:500)
OOB <- data.frame(tree, RF_tree)
ggplot(data = OOB, aes(x = tree))+geom_line(aes(y = OOB), colour = "black", size = 0.8)+geom_line(aes(y = X0), colour = "red", size = 0.8)+geom_line(aes(y = X1), colour = "green", size = 0.8) + labs(y = "Error.rate") + theme_bw()
下面再来看一下，数据中各个变量对于结果的重要性。
importance <- importance(RF, type = 1)
gini <- importance(RF, type = 2)> importance
MeanDecreaseAccuracy
Pclass_train 28.182566
Sex_train 112.706638
Age_train 18.021392
Sibsp_train 13.734762
Parch_train 11.630621
Fare_train 25.760278
Embarked_train 5.642108
> gini
MeanDecreaseGini
Pclass_train 23.251678
Sex_train 169.793011
Age_train 48.646257
Sibsp_train 13.081019
Parch_train 13.158235
Fare_train 60.050685
Embarked_train 7.047434
P_RF <- predict(RF, Random_data[random,])
CT_RF <- table(Random_data[random,1], P_RF)
accuracy_RF <- sum(diag(CT_RF))/sum(CT_RF)*100> accuracy_RF
[1] 91.34615
4.3 支撑向量机的模型训练及预测
支撑向量机分为支撑向量分类器和支撑向量回归机，在这里采用支撑向量分类器进行机器学习。当对全部数据进行训练时，支撑向量机的预测正确率较低，为了提高预测正确率，在这里把Fare和Embarked这两个由其他输入变量所影响的冗余变量去掉，经过训练得出的模型可以大幅提高预测的正确率。
library(e1071)
Svm <- tune.svm(data.frame(Random_data[-random,2:6]), factor(Random_data[-random,1]), data = data.frame(Random_data[-random,]), type = "C-classification", kernel = "radial", cost = c(0.001,0.01,0.1,1,5,10,100,1000), gamma = c(0.5,1,2,3,4), scale = FALSE)> Svm
Parameter tuning of ‘svm’:
- sampling method: 10-fold cross validation
- best parameters:
gamma cost
0.5 1
- best performance: 0.1767445
Svmfit <- svm(data.frame(Random_data[-random,2:6]), factor(Random_data[-random,1]), data = data.frame(Random_data[-random,]), type = "C-classification", kernel = "radial", cost = 1, gamma = 0.5, scale = FALSE)
Svm_P <- predict(Svmfit, data.frame(Random_data[random,2:6]))
CT_svm <- table(Random_data[random,1], Svm_P)
accuracy_svm <- sum(diag(CT_svm))/sum(CT_svm)*100> accuracy_svm
[1] 92.30769
6987
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
