Titanic生存预测
——数据模型汇总报告
摘要
R语言多元统计分析课程是一门综合理论和实践的大课程，既需要我们掌握基本的多元统计分析技术理论，又需要针对具体问题在R的环境中实现。
本文在基本的多元统计分析技术理论基础上，结合机器学习基本模型，挑选了Kaggle（数据建模竞赛网站）的入门赛——Titanic生存预测作为实战演练，较为完整地呈现了数据建模的基本流程和思路。
本文实现的模型包括简单基础模型——递归分割树，条件推理树，K近邻分类器，逻辑回归，朴素贝叶斯分类器以及进阶的集成学习方法——Bagging、Boosting和Random Forest。
建立了以上多种模型后，基于混淆矩阵的模型评估方法给出了Titanic生存预测的基本结论。
关键词： 机器学习 R语言 Kaggle入门赛
目录
Titanic生存预测
摘要
一 问题重述
二 数据建模流程
三 数据说明
四 数据预处理
4.1 缺失值处理
4.1.1 Age缺失值处理
4.1.2 Embarked缺失值处理
4.2 异常值处理
4.2.1 Fare异常值处理
五 描述性分析
5.1 整体数据情况
5.2 生存因素分析
六 模型与评价
6.1 基本模型
6.1.1 递归分割树
6.1.2 条件推理树
6.1.3 K近邻分类器
6.1.4 逻辑回归
6.1.5 朴素贝叶斯分类器
6.1.6 模型评价
6.2 集成学习
6.2.1 Bagging
6.2.2 Boosting
6.2.3 Random Forest
6.2.4 模型评价
七 参考文献
八 程序代码
一 问题重述
泰坦尼克号沉没是世界历史上最臭名昭著的沉船事件之一。 1912年4月15日，泰坦尼克号的处女航期间，在大西洋上撞上一座冰山并于2小时后沉没，2224名乘客和船员中有1502人遇难。 这一耸人听闻的悲剧震撼了国际社会，并直接催生了更好的船舶安全条例。
沉船原因之一是乘客和船员没有足够的救生艇。 虽然幸存下来的人中运气占了一些比例，但毫无疑问有些人更有可能生存下来，比如妇女，儿童和上层阶级。
在这个挑战中，我们要求你分析什么样的人可能生存。 特别是，我们要求你运用机器学习的方法来预测哪些乘客更容易在事故中幸存下来。
二 数据建模流程
三 数据说明
|
变量名 |
含义 |
|
Survived |
是否幸存 |
|
Name |
乘客姓名 |
|
Sex |
乘客性别 |
|
Age |
乘客年龄 |
|
SibSp |
乘客随行的兄弟姐妹数量 |
|
Parch |
乘客随行的父母/兄弟数量 |
|
Ticket |
票号 |
|
Fare |
票价 |
|
Cabin |
船舱 |
|
Pclass |
乘客等级 1=头等 2=二等 3=三等 |
|
Embarked |
登船港口 C = Cherbourg S = Southampton Q = Queenstown |
四 数据预处理
4.1 缺失值处理
在RStudio中读入数据后，直接加载Amelia包可视化缺失值情况，如下：
可以发现：Cabin值几乎都缺失，Age有一部分缺失，Embarked缺失两条。
4.1.1 Age缺失值处理
我发现乘客姓名一列的信息给出了对该乘客的称呼，如下图：
我们知道一个人的称呼不仅可以反映一个人的性别，更可以在一定程度上反映一个人的年龄情况，因此我们先找出缺失年龄的观测所对应的称呼，如下表：
|
缺失值对应称呼 |
缺失数目 |
|
Dr |
1 |
|
Miss |
36 |
|
Master |
4 |
|
Mr |
119 |
|
Mrs |
17 |
因此我采用的方法就是用相应缺失值对应称呼的年龄均值去插补
4.1.2 Embarked缺失值处理
Embarked变量有三个类型，而缺失值只有2个，所以直接采用简单的众数插补法
4.2 异常值处理
4.2.1 Fare异常值处理
在查看各个变量的过程中，无意中发现了Fare即票价变量中的数据有0，总共有6个0数据，在这个场景中我认为0数据是异常的，但不是做简单的删除该观测，而是选择了与年龄插补类似的方法，用乘客等级的均值去插补0数据。
五 描述性分析
5.1 整体数据情况
1)
整体生还情况
2)
出行状况
3)
乘客年龄分布
5.2 生存因素分析
1)
生存——性别关系
可以发现女性生还的比例远远大于男性，这和当时绅士们的谦让有着密不可分的联系
2) 生存——乘客等级关系
可以发现随着乘客等级的下降，遇难率在升高，所以可以反映出当时还是阶层高的人有着更多的优先权
3) 生存——登船港口关系
图中反映出在
S港登船的乘客遇难率更高，而C港登船乘客的幸存率更高。但我们知道泰坦尼克号是在大西洋沉没的，理论上应该与乘客在哪里登船没有直接关系。当然了，这只是一种描述性分析，并不能反映百分之百的真实情况。
4) 生存——出行情况关系
图是根据是否有子女/父母/兄弟姐妹陪同的数据画的，可以看出独自出行的危险性更高。
5) 生存——年龄关系，生存——票价关系
可以看出，年龄越小的乘客，票价越高的乘客，生还概率都越大。
六 模型建立与评价
6.1 基本模型
6.1.1 递归分割树
也称经典决策树，以一个二元输出变量和一组预测变量为基础。具体算法如下：
(1)选定一个最佳预测变量将全部样本单元分为两类，实现两类中的纯度最大。如果预测变量连续，则选定一个分割点进行分类，使得两类纯度最大化；如果预测变量为分类变量，则对各类别进行合并再分类。
(2) 对每一个子类别继续执行步骤(1)。
(3) 重复步骤(1)~(2)，直到子类别中所含的样本单元数过少，或者没有分类法能将不纯度下降到一个给定阈值以下。最终集中的子类别即终端节点（terminal node）。根据每一个终端节点中样本单元的类别数众数来判别这一终端节点的所属类别。
(4) 对任一样本单元执行决策树，得到其终端节点，即可根据步骤3得到模型预测的所属类别。
不过，上述算法通常会得到一棵过大的树，从而出现过拟合现象。结果就是，对于训练集外单元的分类性能较差。为解决这一问题，可采用10折交叉验证法选择预测误差最小的树。这一剪枝后的树即可用于预测。
R中的rpart包支持rpart()函数构造决策树， prune()函数对决策树进行剪枝。
6.1.2 条件推理树
条件推断树与传统决策树类似，但变量和分割的选取是基于显著性检验的，而不是纯净度或同质性一类的度量。显著性检验是置换检验。
条件推断树的算法如下。
(1) 对输出变量与每个预测变量间的关系计算p值。
(2) 选取p值最小的变量。
(3) 在因变量与被选中的变量间尝试所有可能的二元分割（通过排列检验），并选取最显著的分割。
(4) 将数据集分成两群，并对每个子群重复上述步骤。
(5) 重复直至所有分割都不显著或已到达最小节点为止。
条件推断树可由party包中的ctree()函数获得。
6.1.3 K近邻分类器
K最近邻(k-Nearest Neighbor，KNN)分类算法，是一个理论上比较成熟的方法，也是最简单的机器学习算法之一。该方法的思路是：如果一个样本在特征空间中的k个最相似(即特征空间中最邻近)的样本中的大多数属于某一个类别，则该样本也属于这个类别。
K 近邻算法使用的模型实际上对应于对特征空间的划分。K 值的选择，距离度量和分类决策规则是该算法的三个基本要素：
(1)K 值的选择会对算法的结果产生重大影响。K值较小意味着只有与输入实例较近的训练实例才会对预测结果起作用，但容易发生过拟合；如果 K 值较大，优点是可以减少学习的估计误差，但缺点是学习的近似误差增大，这时与输入实例较远的训练实例也会对预测起作用，使预测发生错误。在实际应用中，K 值一般选择一个较小的数值，通常采用交叉验证的方法来选择最优的 K 值。随着训练实例数目趋向于无穷和 K=1 时，误差率不会超过贝叶斯误差率的2倍，如果K也趋向于无穷，则误差率趋向于贝叶斯误差率。
(2)该算法中的分类决策规则往往是多数表决，即由输入实例的 K 个最临近的训练实例中的多数类决定输入实例的类别
(3)距离度量一般采用 Lp 距离，当p=2时，即为欧氏距离，在度量之前，应该将每个属性的值规范化，这样有助于防止具有较大初始值域的属性比具有较小初始值域的属性的权重过大。
KNN算法不仅可以用于分类，还可以用于回归。通过找出一个样本的k个最近邻居，将这些邻居的属性的平均值赋给该样本，就可以得到该样本的属性。更有用的方法是将不同距离的邻居对该样本产生的影响给予不同的权值(weight)，如权值与距离成反比。 该算法在分类时有个主要的不足是，当样本不平衡时，如一个类的样本容量很大，而其他类样本容量很小时，有可能导致当输入一个新样本时，该样本的K个邻居中大容量类的样本占多数。 该算法只计算“最近的”邻居样本，某一类的样本数量很大，那么或者这类样本并不接近目标样本，或者这类样本很靠近目标样本。无论怎样，数量并不能影响运行结果。可以采用权值的方法（和该样本距离小的邻居权值大）来改进。
该方法的另一个不足之处是计算量较大，因为对每一个待分类的文本都要计算它到全体已知样本的距离，才能求得它的K个最近邻点。目前常用的解决方法是事先对已知样本点进行剪辑，事先去除对分类作用不大的样本。该算法比较适用于样本容量比较大的类域的自动分类，而那些样本容量较小的类域采用这种算法比较容易产生误分。
实现 K 近邻算法时，主要考虑的问题是如何对训练数据进行快速 K 近邻搜索，这在特征空间维数大及训练数据容量大时非常必要。
6.1.4 逻辑回归
逻辑回归（logistic regression）是广义线性模型的一种，可根据一组数值变量预测二元输出，R中的基本函数glm()可用于拟合逻辑回归模型。 glm()函数自动将预测变量中的分类变量编码为相应的虚拟变量。
6.1.5 朴素贝叶斯分类器
贝叶斯分类算法是统计学的一种分类方法，它是一类利用概率统计知识进行分类的算法。在许多场合，朴素贝叶斯(Na?ve Bayes，NB)分类算法可以与决策树和神经网络分类算法相媲美，该算法能运用到大型数据库中，而且方法简单、分类准确率高、速度快。
由于贝叶斯定理假设一个属性值对给定类的影响独立于其它属性的值，而此假设在实际情况中经常是不成立的，因此其分类准确率可能会下降。为此，就衍生出许多降低独立性假设的贝叶斯分类算法，如TAN(tree augmented Bayes network)算法。
6.1.6 模型评价
上述5种模型优缺点，如下表：
|
算法 |
优势 |
不足 |
|
递归 分割树 |
无参 能同时解决分类和回归问题 灵活易于理解 |
易产生偏置及过度适应 |
|
条件 推理树 |
无参 能同时解决分类和回归问题 灵活易于理解 偏置问题不如递归分割树严重 |
易产生过度适应 |
|
K近邻 分类器 |
无参 可处理任意类型数据 |
分类结果不容易解释 不能处理连续变量缺失 算法性能依赖于数据维度 |
|
逻辑回归 |
容易理解 可知置信区间和逻辑概率 更新速度快 |
不能处理多重共线性 不能处理连续变量缺失 对连续变量异常值比较敏感 |
|
朴素贝叶斯分类器 |
算法简单，容易理解 适合训练数据集较小的情况 可以处理一些噪声和缺失数据 方便获得预测概率 |
独立性假设不太可能成立 当训练样例增加时，分类会逐渐产生偏离 |
利用混淆矩阵得到5种模型的精确度，如下表：
|
算法 |
精确度 |
|
递归分割树 |
0.8125 |
|
条件推理树 |
0.8125 |
|
K近邻分类器 |
0.6691 |
|
逻辑回归 |
0.625 |
|
朴素贝叶斯分类器 |
0.7757 |
6.2 集成学习
6.2.1 Bagging
Bagging是一种用来提高学习算法准确度的方法，这种方法通过构造一个预测函数系列，然后以一定的方式将它们组合成一个预测函数。Bagging要求“不稳定”（不稳定是指数据集的小的变动能够使得分类结果的显著的变动）的分类方法。
基本思想：
(1)给定一个弱学习算法,和一个训练集
(2)单个弱学习算法准确率不高
(3)将该学习算法使用多次,得出预测函数序列,进行投票
(4)最后结果准确率将得到提高
算法：
(1)For t = 1, 2, …, T Do
从数据集S中取样（放回选样）
训练得到模型Ht
对未知样本X分类时,每个模型Ht都得出一个分类，得票最高的即为未知样本X的分类
(2)也可通过得票的平均值用于连续值的预测
6.2.2 Boosting
Boosting方法是一种用来提高弱分类算法准确度的方法,这种方法通过构造一个预测函数系列,然后以一定的方式将他们组合成一个预测函数。Boosting是一种提高任意给定学习算法准确度的方法。它的思想起源于 Valiant提出的 PAC ( Probably Approxi mately Correct)学习模型。
Boosting方法是一种用来提高弱分类算法准确度的方法,这种方法通过构造一个预测函数系列,然后以一定的方式将他们组合成一个预测函数。他是一种框架算法,主要是通过对样本集的操作获得样本子集,然后用弱分类算法在样本子集上训练生成一系列的基分类器。他可以用来提高其他弱分类算法的识别率,也就是将其他的弱分类算法作为基分类算法放于Boosting 框架中,通过Boosting框架对训练样本集的操作,得到不同的训练样本子集,用该样本子集去训练生成基分类器;每得到一个样本集就用该基分类算法在该样本集上产生一个基分类器,这样在给定训练轮数 n 后,就可产生 n 个基分类器,然后Boosting框架算法将这 n个基分类器进行加权融合,产生一个最后的结果分类器,在这 n个基分类器中,每个单个的分类器的识别率不一定很高,但他们联合后的结果有很高的识别率,这样便提高了该弱分类算法的识别率。在产生单个的基分类器时可用相同的分类算法,也可用不同的分类算法,这些算法一般是不稳定的弱分类算法。
6.2.3 Random Forest
随机森林（random forest）是一种组成式的有监督学习方法。在随机森林中，我们同时生成多个预测模型，并将模型的结果汇总以提升分类准确率。
随机森林的算法涉及对样本单元和变量进行抽样，从而生成大量决策树。对每个样本单元来说，所有决策树依次对其进行分类。所有决策树预测类别中的众数类别即为随机森林所预测的这一样本单元的类别。
假设训练集中共有N个样本单元， M个变量，则随机森林算法如下。
(1) 从训练集中随机有放回地抽取N个样本单元，生成大量决策树。
(2) 在每一个节点随机抽取m<M个变量，将其作为分割该节点的候选变量。每一个节点处的变量数应一致。
(3) 完整生成所有决策树，无需剪枝（最小节点为1）。
(4) 终端节点的所属类别由节点对应的众数类别决定。
(5) 对于新的观测点，用所有的树对其进行分类，其类别由多数决定原则生成。
生成树时没有用到的样本点所对应的类别可由生成的树估计，与其真实类别比较即可得到袋外预测（out-of-bag， OOB）误差。无法获得验证集时，这是随机森林的一大优势。随机森林算法可计算变量的相对重要程度。
randomForest包中的randomForest()函数可用于生成随机森林。函数默认生成500棵树，并且默认在每个节点处抽取sqrt(M)个变量，最小节点为1。
6.2.4 模型评价
分别建模，可以得到各个模型的误差情况，如下表：
|
算法 |
误差 |
|
Bagging |
0.1971 |
|
Boosting |
0.1761 |
|
Random Forest |
0.1777 |
[1]Robert I. Kabacoff. R语言实战[M]. 北京：人民有电出版社.
[2]蔡桂宏，洪锦魁. R语言迈向大数据之路[M]. 北京：清华大学出版社
[3]丘祐玮. 机器学习与R语言实战[M]. 北京：机械工业出版社
八 程序代码
1) 训练数据处理代码
rm(list=ls())
setwd("D:\\R\\R with Multivariate Statistical Analysis\\Titanic生存预测")
#读取数据
train.data <- read.csv("原始数据\\train.csv",na.strings = c("NA",""))
#数据类型转换
train.data$Survived <- as.factor(train.data$Survived)
train.data$Pclass <- as.factor(train.data$Pclass)
str(train.data)
#查看各变量缺失值情况
sapply(train.data,function(df)
{sum(is.na(df) == TRUE)/length(df);})
dim(train.data)[1]
sum(is.na(train.data$Embarked))
sum(is.na(train.data$Age))
sum(is.na(train.data$Cabin))
#缺失数据可视化
library("Amelia")
missmap(train.data,main = "Missing Map")
library(VIM)
aggr(train.data,numbers = TRUE)
matrixplot(train.data)
#Embarked缺失值处理
table(train.data$Embarked,useNA = "always")
train.data$Embarked[which(is.na(train.data$Embarked))] <- ‘S‘
table(train.data$Embarked,useNA = "always")
#Age缺失值处理---根据称呼插补年龄
library(stringr)
train.data$Name <- as.character(train.data$Name)
table_words <- table(unlist(strsplit(train.data$Name,"\\s+")))
sort(table_words [grep(‘\\.‘,names(table_words))],decreasing = TRUE)
tb <- cbind(train.data$Age,str_match(train.data$Name,"[a-zA-Z]+\\."))
table(tb[is.na(tb[,1]),2])
mean.mr <- mean(train.data$Age[grepl("Mr\\.",train.data$Name) & !is.na(train.data$Age)])
mean.mrs <- mean(train.data$Age[grepl("Mrs\\.",train.data$Name) & !is.na(train.data$Age)])
mean.dr <- mean(train.data$Age[grepl("Dr\\.",train.data$Name) & !is.na(train.data$Age)])
mean.miss <- mean(train.data$Age[grepl("Miss\\.",train.data$Name) & !is.na(train.data$Age)])
mean.master <- mean(train.data$Age[grepl("Master\\.",train.data$Name) & !is.na(train.data$Age)])
train.data$Age[grepl("Mr\\.",train.data$Name) & is.na(train.data$Age)] <- mean.mr
train.data$Age[grepl("Mrs\\.",train.data$Name) & is.na(train.data$Age)] <- mean.mrs
train.data$Age[grepl("Dr\\.",train.data$Name) & is.na(train.data$Age)] <- mean.dr
train.data$Age[grepl("Miss\\.",train.data$Name) & is.na(train.data$Age)] <- mean.miss
train.data$Age[grepl("Master\\.",train.data$Name) & is.na(train.data$Age)] <- mean.master
train.data$Name <- as.factor(train.data$Name)
#处理后缺失值可视化
missmap(train.data,main = "Missing Map")
aggr(train.data,numbers = TRUE)
matrixplot(train.data)
#异常值处理
boxplot(train.data$Fare,main = "Fare distribution",ylab = "Fare")
library(vioplot)
vioplot(train.data$Fare,col = "gold")
title("Fare distribution of Vioplot",ylab = "Fare")
summary(train.data$Fare)
#Fare为0
##查看票价和乘客等级关系
x1 <- train.data$Fare[train.data$Pclass == 1]
x2 <- train.data$Fare[train.data$Pclass == 2]
x3 <- train.data$Fare[train.data$Pclass == 3]
vioplot(x1,x2,x3, col = "gold",
names = c("Pclass 1","Pclass 2","Pclass 3"))
title("Correlation of Fare and Pclass",
xlab = "Class of Pclass",
ylab = "Fare")
boxplot(Fare~Pclass,data = train.data,
main = "Correlation of Fare and Pclass",
xlab = "Class of Pclass",ylab = "Fare",
col = c("gold","darkgreen"),
varwidth = TRUE)
##票价为0的用各自等级的均值重填
train.data$Fare[train.data$Fare < 1 &
train.data$Pclass == 1] <- mean(x1)
train.data$Fare[train.data$Fare < 1 &
train.data$Pclass == 2] <- mean(x2)
train.data$Fare[train.data$Fare < 1 &
train.data$Pclass == 3] <- mean(x3)
##再看关系图
y1 <- train.data$Fare[train.data$Pclass == 1]
y2 <- train.data$Fare[train.data$Pclass == 2]
y3 <- train.data$Fare[train.data$Pclass == 3]
vioplot(y1,y2,y3, col = "gold",
names = c("Pclass 1","Pclass 2","Pclass 3"))
title("Correlation of Fare and Pclass",
xlab = "Class of Pclass",ylab = "Fare")
boxplot(Fare~Pclass,data = train.data,
main = "Correlation of Fare and Pclass",
xlab = "Class of Pclass",ylab = "Fare",
col = c("gold","darkgreen"),
varwidth = TRUE)
#异常值缺失值处理结束后查看
summary(train.data)
#没有问题数据了
#保存文件
str(train.data)
setwd("D:\\R\\R with Multivariate Statistical Analysis\\Titanic生存预测\\处理后数据")
write.csv(train.data,file = "trained.csv")
2) 训练数据处理代码
rm(list=ls())
setwd("D:\\R\\R with Multivariate Statistical Analysis\\Titanic生存预测")
#读取数据
test.data <- read.csv("原始数据\\test.csv",na.strings = c("NA",""))
#数据类型转换
test.data$Pclass <- as.factor(test.data$Pclass)
#查看各变量缺失值情况
sapply(test.data,function(df)
{sum(is.na(df) == TRUE)/length(df);})
#Age缺失值处理--根据称呼插补年龄
library(stringr)
test.data$Name <- as.character(test.data$Name)
table_words <- table(unlist(strsplit(test.data$Name,"\\s+")))
sort(table_words [grep(‘\\.‘,names(table_words))],decreasing = TRUE)
tb <- cbind(test.data$Age,str_match(test.data$Name,"[a-zA-Z]+\\."))
table(tb[is.na(tb[,1]),2])
mean.mr1 <- mean(test.data$Age[grepl("Mr\\.",test.data$Name) & !is.na(test.data$Age)])
mean.mrs1 <- mean(test.data$Age[grepl("Mrs\\.",test.data$Name) & !is.na(test.data$Age)])
mean.dr1 <- mean(test.data$Age[grepl("Dr\\.",test.data$Name) & !is.na(test.data$Age)])
mean.miss1 <- mean(test.data$Age[grepl("Miss\\.",test.data$Name) & !is.na(test.data$Age)])
mean.master1 <- mean(test.data$Age[grepl("Master\\.",test.data$Name) & !is.na(test.data$Age)])
test.data$Age[grepl("Mr\\.",test.data$Name) & is.na(test.data$Age)] <- mean.mr1
test.data$Age[grepl("Mrs\\.",test.data$Name) & is.na(test.data$Age)] <- mean.mrs1
test.data$Age[grepl("Dr\\.",test.data$Name) & is.na(test.data$Age)] <- mean.dr1
test.data$Age[grepl("Miss\\.",test.data$Name) & is.na(test.data$Age)] <- mean.miss1
test.data$Age[grepl("Master\\.",test.data$Name) & is.na(test.data$Age)] <- mean.master1
####插补test唯一Age缺失值
sum(is.na(test.data$Age))
class(test.data$Age)
test.data[is.na(test.data$Age),]
test.data$Age[is.na(test.data$Age)] <- mean.miss1
#Fare缺失值处理
mean.1 <- mean(test.data$Pclass == 1 & !is.na(test.data$Fare))
mean.2 <- mean(test.data$Pclass == 2 & !is.na(test.data$Fare))
mean.3 <- mean(test.data$Pclass == 3 & !is.na(test.data$Fare))
test.data$Fare[test.data$Pclass == 1 & is.na(test.data$Fare)] <- mean.1
test.data$Fare[test.data$Pclass == 2 & is.na(test.data$Fare)] <- mean.2
test.data$Fare[test.data$Pclass == 3 & is.na(test.data$Fare)] <- mean.3
aggr(test.data,numbers = TRUE)
#保存文件
setwd("D:\\R\\R with Multivariate Statistical Analysis\\Titanic生存预测\\处理后数据")
write.csv(test.data,file = "tested.csv")
3) 描述性分析代码
rm(list=ls())
setwd("D:\\R\\R with Multivariate Statistical Analysis\\Titanic生存预测\\处理后数据")
#读取数据
train.data <- read.csv("trained.csv",na.strings = c("NA",""))
str(train.data)
train.data$Survived <- as.factor(train.data$Survived)
train.data$Pclass <- as.factor(train.data$Pclass)
#总体生还情况
library(plotrix)
s <- sum(train.data$Survived == 1)
p <- sum(train.data$Survived == 0)
a <- c(s,p)
pct <- round(a/(dim(train.data)[1])*100)
label <- c("Survived","Perished")
label2 <- paste(label," ",pct,"%",sep = "")
pie3D(a,labels = label2)
title("Passanger Survival")
#乘客年龄分布
summary(train.data$Age)
d <- density(train.data$Age)
plot(d,main = "Distribution of Age",font.main = 3)
polygon(d,col = sample(colors(),1),border = "white")
rug(train.data$Age,col = "black")
#亲属情况
train.data$Relative <- train.data$SibSp + train.data$Parch
summary(train.data$Relative)
barplot(table(train.data$Relative),
main = "Passenger Relatives")
train.data$Relative <- as.character(train.data$Relative)
train.data$Relative[train.data$Relative == "0"] <- "Alone"
train.data$Relative[train.data$Relative != "0"
&train.data$Relative != "Alone"] <- "More"
train.data$Relative <- as.factor(train.data$Relative)
summary(train.data$Relative)
barplot(table(train.data$Relative),
main = "Passenger Relatives")
s1 <- sum(train.data$Relative == "Alone")
p1 <- sum(train.data$Relative != "Alone")
a1 <- c(s1,p1)
pct1 <- round(a1/(dim(train.data)[1])*100)
lab1 <- c("Alone","More")
lab2 <- paste(lab1," ",pct1,"%",sep = "")
pie3D(a1,labels = lab2)
title("Passanger Accompany")
#生还因素描述分析
##性别
library(vcd)
spineplot(table(train.data$Sex,train.data$Survived),
main = "Survival by Sex")
##乘客等级
spineplot(table(train.data$Pclass,train.data$Survived),
main = "Survival by Pclass")
##港口
spineplot(table(train.data$Embarked,train.data$Survived),
main = "Survival by Embarked")
##是否单独
spineplot(table(train.data$Relative,train.data$Survived),
main = "Survival by Relative")
##年龄
par(mfrow = c(1:2))
boxplot(train.data$Age~train.data$Survived,
main = "Survival by Age",
xlab = "Survived",ylab = "Age")
z1 <- train.data$Age[train.data$Survived == 0]
z2 <- train.data$Age[train.data$Survived == 1]
vioplot(z1,z2,col = "gold",
names = c("Perished","Survived"))
title("Survival by Age",ylab = "Age")
##票价
boxplot(train.data$Fare~train.data$Survived,
main = "Survival by Fare",
xlab = "Survived",ylab = "Fare")
z3 <- train.data$Fare[train.data$Survived == 0]
z4 <- train.data$Fare[train.data$Survived == 1]
vioplot(z3,z4,col = "gold",
names = c("Perished","Survived"))
title("Survival by Fare",ylab = "Fare")
4) 模型代码
rm(list=ls())
setwd("D:\\R\\R with Multivariate Statistical Analysis\\Titanic生存预测\\处理后数据")
#读取数据
train.data <- read.csv("trained.csv",na.strings = c("NA",""))
str(train.data)
train.data$Survived <- as.factor(train.data$Survived)
train.data$Pclass <- as.factor(train.data$Pclass)
test.data <- read.csv("tested.csv",na.strings = c("NA",""))
test.data$Pclass <- as.factor(test.data$Pclass)
#划分数据
data <- train.data[,!names(train.data) %in% c("X","PassengerId","Name","Ticket","Cabin")]
set.seed(2)
index <- sample(2,nrow(data),replace = TRUE,prob = c(0.7,0.3))
trainset <- data[index == 1,]
testset <- data[index == 2,]
#递归分割树
library(rpart)
data.rp <- rpart(Survived ~ .,data = trainset)
data.rp
printcp(data.rp)
plotcp(data.rp)
summary(data.rp)
#可视化树
plot(data.rp,uniform = TRUE,branch = 0.6,margin = 0.2)
text(data.rp,all = TRUE,use.n = TRUE)
#评估树
predictions <- predict(data.rp,testset,type = "class")
table(testset$Survived,predictions)
library(caret)
confusionMatrix(table(predictions,testset$Survived))
#条件推理树
library(party)
train.ctree <- ctree(Survived~.,data = trainset)
train.ctree
plot(train.ctree,main = "Tree")
#评估条件推理树
ctree.predict <- predict(train.ctree,testset)
confusionMatrix(ctree.predict,testset$Survived)
library(ROCR)
train.ctree.pred <- predict(train.ctree,testset)
train.ctree.prob <- 1 - unlist(treeresponse(train.ctree,testset),
use.names = F)[seq(1,nrow(testset)*2,2)]
train.ctree.prob.rocr <- prediction(train.ctree.prob,testset$Survived)
train.ctree.perf <- performance(train.ctree.prob.rocr,"tpr","fpr")
train.ctree.auc.perf <- performance(train.ctree.prob.rocr,
measure = "auc",x.measure = "cutoff")
plot(train.ctree.perf,col = 2,colorize = T,
main = paste("AUC:",[email protected]))
#Logistic回归分类
fit <- glm(Survived~.,data = trainset,
family = binomial())
summary(fit)
fit <- glm(Survived~Pclass + Sex + Age + SibSp,
data = trainset,
family = binomial())
summary(fit)
pred <- predict(fit,testset,type = "response")
Class <- pred > .5
summary(Class)
tb <- table(testset$Survived,Class)
tb
fit.mod <- ifelse(testset$Survived == "1",1,0)
pred.class <- fit.mod
pred.class[pred>=.5] <- 1-pred.class[pred>=.5]
ctb <- table(fit.mod,pred.class)
ctb
confusionMatrix(ctb)
#贝叶斯分类器
library(e1071)
classfier <- naiveBayes(trainset[,!names(trainset)%in%c("Survived")],
trainset$Survived)
classfier
bayes.table <-table(predict(classfier,testset[,!names(trainset)%in%c("Survived")]),
testset$Survived)
bayes.table
confusionMatrix(bayes.table)
#K近邻
library(class)
levels(trainset$Sex) <- list("0" = "male","1" = "female")
levels(trainset$Embarked) <- list("1" = "C","2" = "Q","3" = "S")
levels(testset$Sex) <- list("0" = "male","1" = "female")
levels(testset$Embarked) <- list("1" = "C","2" = "Q","3" = "S")
s.knn <- knn(trainset[,!names(trainset)%in%c("Survived")],
testset[,!names(trainset)%in%c("Survived")],
trainset$Survived,k = 3)
summary(s.knn)
table(testset$Survived,s.knn)
confusionMatrix(table(testset$Survived,s.knn))
#bagging
library(adabag)
set.seed(2)
s.bagging <- bagging(Survived~.,data = trainset,mfinal = 10)
s.bagging$importance
s.pre <- predict.bagging(s.bagging,newdata = testset)
s.pre$confusion
s.pre$error
s.bagging.cv <- bagging.cv(Survived~.,v = 10,
data = trainset,mfinal = 10)
s.bagging.cv$confusion
s.bagging.cv$error
#boosting
set.seed(2)
s.boosting <- boosting(Survived~.,data = trainset,
mfinal = 10,coeflearn = "Freund",
boos = FALSE,
control = rpart.control(maxdepth = 3))
s.boosting.pre <- predict.boosting(s.boosting,newdata = testset)
s.boosting.pre$confusion
s.boosting.pre$error
s.boosting.cv <- boosting.cv(Survived~.,v=10,data = trainset,
mfinal = 5,control = rpart.control(cp=0.01))
s.boosting.cv$confusion
s.boosting.cv$error
#随机森林
library(randomForest)
s.rf <- randomForest(Survived~.,data = trainset, importance = T)
s.rf
s.rf.pre <- predict(s.rf,testset)
table(s.rf.pre,testset$Survived)
plot(s.rf)
importance(s.rf)
varImpPlot(s.rf)
margin.rf <- margin(s.rf,trainset)
plot(margin.rf)
library(tibble)
library(ggplot2)
library(caret)
importance <- varImp(s.boosting,scale = FALSE)
plot(importance)
importance
##############
library(ada)
library(ipred)
library(randomForest)
suv.bagging <- errorest(Survived~.,data = trainset,
model = bagging)
suv.bagging
suv.boosting <- errorest(Survived~.,data = trainset,
model = ada)
suv.boosting
suv.rf <- errorest(Survived~.,data = trainset,
model = randomForest)
suv.rf
suv.pre <- function(object,newdata){
predict(object,newdata = newdata,type = "class")
}
suv.ctree <- errorest(Survived~.,data = trainset,
model = rpart,predict = suv.pre)
suv.ctree
R语言多元统计分析课程是一门综合理论和实践的大课程，既需要我们掌握基本的多元统计分析技术理论，又需要针对具体问题在R的环境中实现。
本文在基本的多元统计分析技术理论基础上，结合机器学习基本模型，挑选了Kaggle（数据建模竞赛网站）的入门赛——Titanic生存预测作为实战演练，较为完整地呈现了数据建模的基本流程和思路。
本文实现的模型包括简单基础模型——递归分割树，条件推理树，K近邻分类器，逻辑回归，朴素贝叶斯分类器以及进阶的集成学习方法——Bagging、Boosting和Random Forest。
建立了以上多种模型后，基于混淆矩阵的模型评估方法给出了Titanic生存预测的基本结论。
关键词： 机器学习 R语言 Kaggle入门赛