1. 选题介绍

　　泰坦尼克号是当时世界上最大的豪华客轮，被称为是“永不沉没的客轮”或是“梦幻客轮”。1912年4月10日，泰坦尼克号从英国南安普敦（Southampton）出发，途经法国瑟堡-奥克特维尔（Cherbourg-Octeville）以及爱尔兰（此时为英属）的昆士敦（Queenstown），计划中的目的地为美国的纽约（NewYork），开始了这艘“梦幻客轮”的处女航。4月14日晚11点40分，泰坦尼克号在北大西洋撞上冰山，两小时四十分钟后，4月15日凌晨2点20分沉没，由于缺少足够的救生艇，1500人葬生海底，造成了当时在和平时期最严重的一次航海事故，也是迄今为止最为人所知的一次海难。本次分析借助Rstudio软件编程，从事故生存率入手，探究影响乘客生还的主要因素，并以乘客是否存活为因变量进行预测分析。

2. 模型介绍

　　本次分析主要从乘客生还率相关描述性统计和生还预测建模两个角度展开，其中生还预测建模主要采用以下两个模型：二分类问题常用的logistic回归模型和包含多个决策树分类器的随机森林模型。

2.1 logistic回归

　　logistic回归是一种广义的线性回归分析模型，常用于数据挖掘，疾病自动诊断，经济预测等领域。在机器学习中，logistic回归是一种有监督的统计学习方法，主要用于对样本进行分类。如线性回归模型\(Y=\beta*X+\sigma\)中，\(Y\)的输出一般是连续的，对于每一个输入的\(X\)，都有一个对应的\(Y\)输出。模型的定义域和值域都可以是\([-\infty,+\infty]\)。但是对于逻辑回归，输入可以是连续的\([-\infty,+\infty]\)，但输出一般是离散的，即只有有限多个输出值。例如，当输出值域只有两个值\(\lbrace 0,1\rbrace\)，这两个值可以表示对样本的某种分类，高/低、患病/健康、阴性/阳性等，这就是最常见的二分类逻辑回归。因此，从整体上来说，通过逻辑回归模型，我们将在整个实数范围上的\(X\)映射到了有限个点上，这样就实现了对\(X\)的分类。因为每次拿过来一个\(X\)，经过逻辑回归分析，就可以将它归入某一类\(Y\)中。

　　logistic回归运用到的假设函数为sigmoid函数，其函数表达式为：\(Sigmoid(x)=\frac{1}{1+exp(x)}\)，其对应的函数图像如下图所示：

　　通过sigmoid函数的作用，logistic模型回归表达式，我们可以将输出的值限制在区间\([0,1]\)上，\(P(X)\)则可以用来表示概率\(P(Y=1|X)\)，即当一个\(X\)发生时，\(Y\)被分到1那一组的概率。在真实情况下，我们最终得到的\(Y\)的值是在\([0,1]\)这个区间上的一个数，然后我们可以选择一个阈值，通常是0.5，当\(y>0.5\)时，就将这个\(x\)归到1这一类，如果\(y<0.5\)就将\(x\)归到0这一类。但是阈值是可以调整的，比如说一个比较保守的人，可能将阈值设为0.9，也就是说有超过90%的把握，才相信这个\(x\)属于1这一类。本次分析中，以实际的生存率作为阈值，对logistic回归输出的概率进行分类。

2.2 随机森林

　　随机森林分类器是基于决策树分类而形成的一种分类器。由于决策树是一种单一的分类器模式，它得到的解往往是局部最优而非全局最优，而且训练时容易出现过拟合的情况。于是，研究分类技术的科研工作者把多个分类器放在一起训练,，使得集成学习一度成为研究的热点。在此期间涌现出许多著名的集成学习方法，比如Bagging、Random Subspace和Boosting等。随机森林是综合考虑多个决策树而形成的一种集成分类器模型，它不仅可用于分类还可以解决回归问题。随机森林的投票决策过程如下：

　　\(H(x)=arg max_{Y}\sum_{i=1}^{k} {I(h_{i}(x))=Y}\)

　　其中，\(H(x)\)表示组合分类模型，\(h_{i}\)表示单棵决策树， \(Y\)为输出变量,\(I(.)\)为指示性函数。算法会根据最大投票法则判断得票数最多的一类作为最后的分类结果。随机森林算法的构建过程主要包括训练集的产生、决策树的构建以及算法的形成和执行三个部分，可以用下图表示其算法过程：

　　按这种算法得到的随机森林中每一棵树都是很弱的，因为我们只从\(M\)个特征中随机选取\(m\)个让每一棵决策树进行学习。但是可以把此时的每棵树比喻成精通于某一小领域的专家，而随机森林则是综合了那些精通不同领域的专家。因此，对于一个分类或回归数据，各个专家可以从不同的角度加以判断，并通过投票得到结果。

　　对于一个规模大小为\(N\)的森林，算法需要训练\(N\)棵决策树，为此需要产生对应数量的训练集。为了使决策树不产生局部最优解，随机森林采用有放回的Bagging随机抽样技术来产生\(N\)个训练样本集合，而这样的操作将导致抽样的训练样本中不可避免的出现重复的情况。

　　单棵决策树的生成过程包括两个：节点分裂和随机特征变量的随机选取。相比于随机特征变量的随机选取，节点分裂是构建决策树的过程中更重要的一步，因为只有通过节点分裂才能产生一棵完整的决策树。由于随机森林采用CART作为基本分类器，所以这里选择的节点分裂是依据Gini不纯度最小准则。计算公式为：\(Gini(t)=1-sum_{j=1}^{e} [p(j|t)]^2 (j=1,\ldots,c)\) 这里\(p(j|t)\)表示节点\(t\)上类别为\(j\)的概率，当节点\(t\)的所有样本属于同一类时，基尼指数取到最小值0，此时样本类别最纯；当基尼指数取最大值1时，处在当前节点的样本类别是最不纯的，即类别不一。

　　随机特征变量是指节点分裂中用于比较最佳属性的特征数目，决策树在节点分裂时会在所有特征中随机选择一部分用于计算最佳分裂属性，即是说并不是所有的特征属性都参与计算，一般取\(m=\log_{2} (M+1)\)或者\(m=\sqrt{M}\)，\(M\)为输入变量的个数，这样可以减少决策树之间的相关性，同时提升每棵决策树的分类精度。待所有决策树按上述方法训练好之后，就生成了随机森林。对一个新的测试样本进行分类，随机森林将综合每棵决策树的分类结果，并采取Gini公式所描述的多数投票法则实现算法的输出结果。

3 数据预处理

3.1 数据说明

　　原始数据共有两个数据集：train训练数据集和test测试数据集。训练集共包含891条观测样本，10个变量。测试集共包含891条待测数据。各变量的具体含义如下表所示：

变量字典
index	meaning
PassengerId	乘客ID
Survived	是否存活
Pclass	乘客等级
Name	乘客姓名
Sex	乘客性别
Age	乘客年龄
SibSp	乘客随行兄弟姐妹数量
Parch	乘客随行父母数量
Ticket	票号
Fare	票价
Cabin	船舱号
Embarked	登船港口

乘客等级字典
index	meaning
1	头等
2	二等
3	三等

乘客等级字典
index	meaning
C	Cherbourg
S	Southampton
Q	Queenstown

　　其中Survived变量为本次的被解释变量，其余为解释变量。通过对变量含义的分析，可以知道Sex性别变量、Pclass乘客等级和Survived是否存活的取值存在一定关系，故将该三个变量进行因子化处理，转化为因子数据，便于后续的数据处理以及建模。由于乘客编号、票号、姓名这三个变量对于乘客是否生还并无影响，只是简单的序列索引号或者标记，故在此对其进行删除操作。至此，训练集和测试集的最终解释变量为8个。

3.2 缺失情况

　　正式分析前，首先对数据的缺失情况进行查看。利用map函数，求解各变量缺失比例，并借助相关包，对数据各变量缺失情况进行统计。

变量缺失比例
Pclass	Name	Sex	Age	SibSp	Parch	Ticket	Fare	Cabin	Embarked
0	0	0	0.2	0	0	0	0	0.77	0.0022

##     Pclass Name Sex SibSp Parch Ticket Embarked Fare Age Cabin    
## 87       1    1   1     1     1      1        1    1   1     1   0
## 244      1    1   1     1     1      1        1    1   1     0   1
## 4        1    1   1     1     1      1        1    1   0     1   1
## 82       1    1   1     1     1      1        1    1   0     0   2
## 1        1    1   1     1     1      1        1    0   1     0   2
##          0    0   0     0     0      0        0    1  86   327 414

变量缺失比例
Pclass	Name	Sex	Age	SibSp	Parch	Ticket	Fare	Cabin	Embarked
0	0	0	0.21	0	0	0	0.0024	0.78	0

　　通过查看训练集的缺失情况可以发现，Cabin船舱号变量的缺失最为严重，缺失率超过77%，其次为Age乘客年龄变量，缺失率接近20%，Embarked登船港口变量仅缺失2条，其余变量均无缺失。就单个样本观测而言，仅有183条样本数据没有任何缺失，仅缺失Cabin船舱号变量的样本高达529个。测试集中，船舱号和年龄缺失依旧较大，登船港口没有缺失，票价仅缺失1条。从实际角度出发，船舱号变量应该对生存与否没有太大影响，且暂时没有较好的方法对其进行填充，故直接删除该变量。

3.3 缺失处理

　　根据变量的类型以及各自的缺失情况，以下将结合数据本身的特点，综合考虑变量间的相关性，进行缺失值处理。

3.3.1 Age年龄缺失值处理

　　从实际情况出发，在西方的文化中，出于对人的尊重，不同年龄和不同身份的人有着不同的称谓，比如Miss指的是年轻女性，Mrs指的是已婚女性，Master指的是少爷，Mr是对男性的尊称。基于年龄与名字之间所存在的一定关系，考虑使用名字称谓来对Age年龄变量缺失进行填补。单纯借助称谓来填充年龄有一定的缺陷，比如Miss和Mr的年龄范围较大，如果缺失值为幼小的女性或者男性，利用年龄均值进行缺失值填充可能会有较大偏差。因此，还可以将乘客随行父母的数量和乘客随行兄弟姐妹的数量一起纳入分析，一般而言，小孩外出需要父母的陪同，不过另一方面来说，未婚成年女子也有可能陪伴父母出游，这一点也是很难分辨的。具体操作为先找出数据集中年龄缺失对应的称谓以及是否有随行父母两个变量，接着对各类称谓的年龄平均值进行计算，最后把计算结果填充到相应的年龄缺失中。对于称谓寻找的过程，一般而言称谓后紧跟着标点符号.，故对Name姓名变量按照.进行称谓提取，这一部分需要用到正则表达式。

　　由于年龄变量在训练集和测试集的缺失都较为严重，为更好地利用数据，简化填充过程，将两个数据集合并之后统一进行年龄缺失值填充，并在填充完成后重新拆分为训练集与测试集。

　　整体称谓统计，年龄缺失对应的称谓统计和各类称谓的年龄均值结果依次如下：

AGE缺失对应的称谓统计
Var1	Freq
Dr.	1
Master.	8
Miss.	50
Mr.	176
Mrs.	27
Ms.	1

AGE缺失称谓平均年龄统计
mr0	mr1	mrs	dr	miss0	miss1	master	ms
32	35	37	44	27	12	5.5	28

　　利用对应称谓的平均年龄作为缺失填充，对所有数据进行年龄缺失填充，年龄变量缺失值填充完成后，拆分合并数据集为测试集和训练集。

3.3.2 Embarked缺失值处理

　　通过前述的缺失值分析，Embarked登船港口变量缺失较少，仅在训练集中有2条缺失，查看具体缺失的样本情况，发现缺失对应数据的其他变量除年龄不一样，其他关键变量均相同。

登船港口统计
Var1	Freq
C	168
Q	77
S	644
NA	2

缺失样本
Pclass	Sex	Age	SibSp	Parch	Fare	Embarked	Survived
1	female	38	0	0	80	NA	1
1	female	62	0	0	80	NA	1

登船港口等级统计
Embarked	mean(Fare)	median(Fare)	n()
C	105	78	85
Q	90	90	2
S	70	52	127
NA	80	80	2

##  /\     /\
## {  `---'  }
## {  O   O  }
## ==>  V <==  No need for mice. This data set is completely observed.
##  \  \|/  /
##   `-----'

##     Pclass Sex Age SibSp Parch Fare Embarked Survived  
## 891      1   1   1     1     1    1        1        1 0
##          0   0   0     0     0    0        0        0 0

　　因为两个样本的等级和票价相同，基于同一登陆地的票价和等级应该有一定规律来反推缺失值。此外票号也可能包含登船港口信息的猜测，可以考虑使用该缺失数据对应票号附近并结合等级票价的情况，对Embarked变量取值对其进行缺失值填补。通过绘制登船港口与票价的箱线图以及与乘客等级的相关图，发现登船港口与票价和乘客等级有一定的关系。经过对乘客等级和票价的分组分析后，选择将Embarked变量缺失值填补为S。

票价缺失统计
Pclass	Sex	Age	SibSp	Parch	Fare	Embarked
3	male	60	0	0	NA	S

登船港口票价统计
Embarked	mean(Fare, na.rm = T)	median(Fare, na.rm = T)	n()
C	11	7.2	35
Q	9	7.8	41
S	14	8.1	142

##  /\     /\
## {  `---'  }
## {  O   O  }
## ==>  V <==  No need for mice. This data set is completely observed.
##  \  \|/  /
##   `-----'

##     Pclass Sex Age SibSp Parch Fare Embarked  
## 418      1   1   1     1     1    1        1 0
##          0   0   0     0     0    0        0 0

　　针对测试集的数据预处理，通过缺失值分析可以发现，数据缺失情况与训练集较为相似，故采用上述相同的办法对其进行缺失值处理。此外，通过对测试集的缺失值分析可以发现，在测试集中Fare变量存在一个缺失值，采取的缺失值补充办法为，以同等级同登船港口的平均值作为缺失值的填充。对于Fare变量中出现0值的情况，不能排除有免费票的情况，故现阶段不对Fare变量进行过多的处理。

　　至此，数据预处理暂时告一段落，缺失值已经全部解决。不同的缺失填充方法可能对于之后的建模有一定的差别，就分析该问题而言，前述的缺失值填充是有一定可行性的。

4. 描述性统计

　　本节主要从描述性统计入手，初步探究影响乘客生还的因素。

4.1 生存率分析

　　首先，通过绘制总体生存率柱状图可以发现，在此次Titanic遇难事件中，乘客的整体存活率为42.75%，故在之后模型中，设定阈值与0.4275。下面以训练集为分析对象，区分离散型变量和连续型变量，通过绘制柱形图和直方图查看各自变量的分布情况，以便确定正确的变量形式，并初步探究各变量与存活率之间的关系。

4.2 离散变量

　　通过Pclass与存活情况的分析可以看到，等级越高的人存活率越高，等级低的人数虽然多，但存活率并不高，这点是符合正常逻辑的，社会等级地位越高的人，受到的服务是更好的，在遇难时也是有一定的空间优势（数字越大，等级越低）。针对性别变量而言，女性的存活率明显高于男性，说明在此次遇难事故中，欧洲的人们普遍秉持着女士优先的传统，绅士作风得到充分体现。

　　此外，还可以继续探究不同乘客等级中，性别在生存优势上的差异。从Pclass-Sex柱形图可以明显的看到，各等级中，男性死亡人数均超过女性；乘客等级3中女性存活人数和死亡人数相差不大，可以猜测等级越低的人群中，在生存面前女性优势相对较小。而对于登船口岸的分析可以发现，登船港口与存活率间有一定的关联，其中在英国始发站Southampton登船的人最多，同时存活人数也最多；在法国Cherbourg登船的乘客存活人数大于死亡人数。

4.3 连续变量

　　通过对年龄绘制直方图可以看到，年龄分布大致符合正态分布，结合生存情况来看，小于10岁的孩子存活率较高，与女性存活较高保持一致，可以发现妇女儿童在事故中有优先施救的可能；而作为中坚力量的30岁阶段的人死亡人数明显高于存活人数。从票价的直方图来看，票价分布明显呈现右偏的趋势，且高价票的生还率明显高于低票价的。

　　从整体来看，陪同人数较少，绝大多数情况下取值为零，故仿照分类变量绘制柱形图进行分析。可以看到，没有同行人员的死亡率高于有同行人员的，同行人员数量过多的生存率近乎为0。不难想象，2人结伴而行在遇到灾难时可以相互帮助，而随行人员过多则可能是一个家庭，救援时主力的男主人并没有办法解救所有人，且在多项决策中更容易同时遇难。

5. 模型建立

　　利用训练集数据进行预测建模分析，为了验证模型的预测效果以及是否存在过拟合，依照交叉验证的原则，首先将训练集按七比三的比例拆分为更细的训练集和测试集。此外，logistic模型和随机森林模型的预测结果均为概率，在转化为是否存活时，以前面设定的阈值0.4275作为划分标准，得到最终的生还预测结果。

5.1 logistic模型

## Error: variable 'Pclass' was fitted with type "factor" but type "numeric" was supplied

## 
## Call:
## glm(formula = Survived ~ ., family = binomial(link = "logit"), 
##     data = trainset)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.807  -0.614  -0.430   0.633   2.426  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.35611    0.59529    7.32  2.5e-13 ***
## Pclass2     -1.09263    0.36475   -3.00   0.0027 ** 
## Pclass3     -1.97023    0.35466   -5.56  2.8e-08 ***
## Sexmale     -2.59638    0.23764  -10.93  < 2e-16 ***
## Age         -0.05068    0.00982   -5.16  2.4e-07 ***
## SibSp       -0.39200    0.13178   -2.97   0.0029 ** 
## Parch       -0.00777    0.13621   -0.06   0.9545    
## Fare         0.00365    0.00310    1.17   0.2401    
## EmbarkedQ   -0.20311    0.46404   -0.44   0.6616    
## EmbarkedS   -0.47850    0.28359   -1.69   0.0915 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 826.64  on 618  degrees of freedom
## Residual deviance: 552.82  on 609  degrees of freedom
## AIC: 572.8
## 
## Number of Fisher Scoring iterations: 5

　　logistic建模过程和具体的模型拟合结果如上代码所示。从logistic回归的结果中可以发现，乘客等级，乘客性别，乘客年龄以及随行人数对于生还率的影响是统计显著的。得到的结论与描述性统计分析基本符合，在其他条件不变的情况下，年龄越小获救率越高，乘客等级越高获救率越高，女性的获救率高于男性，其中性别和等级对获救率的影响贡献度是最大的。

　　从模型拟合的角度来看，模型准确率大致为0.82,模型对于预测新数据而言，准确率也在0.79。总体来说，logistic模型对于预测生还情况效果良好。具体而言，在样本内无论是预测存活还是预测死亡，预测的准确率都略高于样本外的预测；无论是样本内还是样本外，对于死亡的预测准确率是高于预测存活的准确率的，说明模型预测效果较好。

logic模型准确率
logic_accuracy	logic_accuracy2
0.82	0.79

logic预测情况1
	0	1
0	349	84
1	30	156

logic预测情况2
	0	1
0	151	37
1	19	65

5.2 随机森林

## Error in predict.randomForest(rf_mod, test): Type of predictors in new data do not match that of the training data.

## 
## Call:
##  randomForest(formula = Survived ~ ., data = trainset, importance = T) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 18%
## Confusion matrix:
##     0   1 class.error
## 0 344  35       0.092
## 1  75 165       0.312

##             0    1 MeanDecreaseAccuracy MeanDecreaseGini
## Pclass   12.2 24.9                   27             20.5
## Sex      48.7 65.9                   73             67.7
## Age      18.5 21.8                   29             42.2
## SibSp    21.3  4.4                   21             13.5
## Parch    13.2  9.7                   17             11.5
## Fare     18.6 23.8                   30             50.2
## Embarked  6.2 12.4                   14              9.8

　　从随机森林的建模结果可以发现，在影响生存的因素中，性别是最为重要的，其次位乘客等级，再来是票价和年龄。同logistic模型结果进行比较发现，票价的影响是较大的，而陪同人数相较而言甚至比登船港口的重要性还低。通过分析可以知道，票价和乘客等级以及登船地点有较强的相关性，故在建模中会略微凸显。

　　从混淆矩阵中可以发现，随机森林模型对于死亡的预测准确率是高于存活的预测率。

rf模型准确率
rf_accuracy	rf_accuracy2
0.92	0.82

rf预测情况1
	0	1
0	370	41
1	9	199

rf预测情况2
	0	1
0	157	36
1	13	66

　　可以看到随机森林建模的拟合效果,模型预测准确率从logistic模型的82%和79%分别提升到了92%和82%，预测效果显著增加，凸显了随机森林的优势。此外，具体到生存和死亡的预测准确情况，无论是样本内还是样本外的预测准确率明显较Logistic模型有明显的提升，且样本内各情况的预测准确率略高于样本外的预测准确率，对于死亡的预测准确率高于存活的预测准确率。

6. 总结

　　通过对乘客生还率的分析，我们可以清晰地看到，在事故发生时，女性的获救率明显高于男性，说明欧洲人普遍秉持着女士优先的传统；其次，年幼的乘客获救率也明显较高，说明在事故面前，人们尊重女性关爱幼者是言行一致的；再者，从乘客等级来看，社会地位越高的人获救率也更高，这也符合我们的日常认知。

　　本次课程论文充分利用课程中所学的map循环函数，ggplot绘图函数，modelr包建模函数，Rmarkdown文本编辑，对数据进行了相应分析，一定程度上简化了编程内容，提升了编程效率。当然，在本次课程论文的撰写过程中，还是发现仍存在许多的不足，例如对于图片输出的格式控制上未能统一大小，格式化图文排版相对较差，需要通过输出word文档后进行微调，还有进一步提升的空间。

　　总的来说，通过此次的实战演练，熟悉了R语言编程在实际数据分析中的运用，熟练掌握相关函数的使用以及数据可视化的基本能力，循环迭代以及运行速率的提升以及自动化文本编辑能力是日后需要努力的方向。　

附录

library(tidyverse)
library(ggplot2)
library(gridExtra)
library(modelr)
#sigma函数
sigmoid <- function(x){
  1/(1+exp(-x))
}
x <- seq(-5,5,0.01)
s <- tibble(x=x,y=sigmoid(x))
ggplot(s,aes(x=x,y=y))+
  geom_vline(xintercept=0,color='black')+
  geom_smooth(se=F)+
  labs(title='sigmoid函数图像',y=NULL)+
  theme(plot.title = element_text(hjust = 0.5))

#readdata
train <- read_csv("D:/Titanic/train.csv")
test <- read_csv("D:/Titanic/test.csv")
trainx <- train[,-c(1,2)]
test <- test[,-1]
alldata <- rbind(trainx,test)
trainy <- train[,2]
#变量字典
dictionary <- tibble(index=names(train),meaning=c('乘客ID','是否存活','乘客等级','乘客姓名','乘客性别','乘客年龄','乘客随行兄弟姐妹数量','乘客随行父母数量','票号','票价','船舱号','登船港口'))
plclass <- tibble(index=c(1,2,3),meaning=c('头等','二等','三等'))
embarked <- tibble(index=c('C','S','Q'),meaning=c('Cherbourg','Southampton','Queenstown'))
knitr::kable(dictionary,caption='变量字典',align='c')
knitr::kable(plclass,caption='乘客等级字典',align='c')
knitr::kable(embarked,caption='乘客等级字典',align='c')
#缺失分析
library(mice)  
library(Amelia)
library(VIM)
options(digits=2)
ms1 <- map_dbl(trainx,function(df) {sum(is.na(df==TRUE))/length(df)})#缺失率计算
knitr::kable(t(ms1),caption='变量缺失比列',align='c',digits = 4)
missmap(trainx,main = "Missing Map")#整体缺失图

aggr(trainx,numbers = TRUE)#缺失比例图

#test缺失值分析
ms2 <- map_dbl(test,function(df) {sum(is.na(df==TRUE))/length(df)})#缺失率计算
knitr::kable(t(ms2),caption='变量缺失比列',align='c',digits = 4)
md.pattern(test)#缺失矩阵
aggr(test,numbers = TRUE)#缺失比例图

#缺失处理
#Age（alldata填充）
library(stringr)
table_words <- table(unlist(strsplit(alldata$Name,"\\s+")))
#称谓提取
chengwei <- sort(table_words [grep('\\.',names(table_words))],decreasing = TRUE)
#AGE缺失对应的称谓
cheng <- str_match(alldata$Name,"[a-zA-Z]+\\.")
ageque <- cbind(alldata$Age,cheng)
ageque2 <- cbind(age=alldata$Age,touxian=cheng,parch=alldata$Parch)
#缺失称谓的均值
mean_mr0 <- mean(alldata$Age[grepl("Mr\\.",alldata$Name) & !is.na(alldata$Age) & alldata$Parch==0])
mean_mr1 <- mean(alldata$Age[grepl("Mr\\.",alldata$Name) & !is.na(alldata$Age) & alldata$Parch!=0])
mean_mrs <- mean(alldata$Age[grepl("Mrs\\.",alldata$Name) & !is.na(alldata$Age)])
mean_dr <- mean(alldata$Age[grepl("Dr\\.",alldata$Name) & !is.na(alldata$Age)])
mean_miss0 <- mean(alldata$Age[grepl("Miss\\.",alldata$Name) & !is.na(alldata$Age) & alldata$Parch==0])
mean_miss1 <- mean(alldata$Age[grepl("Miss\\.",alldata$Name) & !is.na(alldata$Age) & alldata$Parch!=0])
mean_master <- mean(alldata$Age[grepl("Master\\.",alldata$Name) & !is.na(alldata$Age)])
mean_ms <- mean(alldata$Age[grepl("Ms\\.",alldata$Name) & !is.na(alldata$Age)])
mean_que <- tibble(mr0=mean_mr0,mr1=mean_mr1,mrs=mean_mrs,dr=mean_dr,miss0=mean_miss0,miss1=mean_miss1,master=mean_master,ms=mean_ms)
knitr::kable(table(ageque[is.na(ageque[,1]),2]),caption='AGE缺失对应的称谓统计',align='c')
knitr::kable(mean_que,caption='AGE缺失称谓平均年龄统计',align='c')
#均值填充
alldata$Age[grepl("Mr\\.",alldata$Name) & is.na(alldata$Age) & alldata$Parch==0] <- mean_mr0
alldata$Age[grepl("Mr\\.",alldata$Name) & is.na(alldata$Age) & alldata$Parch!=0] <- mean_mr1
alldata$Age[grepl("Mrs\\.",alldata$Name) & is.na(alldata$Age)] <- mean_mrs
alldata$Age[grepl("Dr\\.",alldata$Name) & is.na(alldata$Age)] <- mean_dr
alldata$Age[grepl("Miss\\.",alldata$Name) & is.na(alldata$Age) & alldata$Parch==0] <- mean_miss0
alldata$Age[grepl("Miss\\.",alldata$Name) & is.na(alldata$Age) & alldata$Parch!=0] <- mean_miss1
alldata$Age[grepl("Master\\.",alldata$Name) & is.na(alldata$Age)] <- mean_master
alldata$Age[grepl("Ms\\.",alldata$Name) & is.na(alldata$Age)] <- mean_ms
alldata <- alldata[,-c(2,7,9)]
trainx <- alldata[1:891,]
trainx <- cbind(trainx,trainy)
test <- alldata[892:1309,]
#Embarked(trainx)
knitr::kable(table(trainx$Embarked,useNA = "always"),caption='登船港口统计',align='c')
knitr::kable(filter(trainx,is.na(Embarked)),caption='整体称谓统计',align='c')
e1 <- ggplot(trainx)+
  geom_boxplot(aes(y=Fare,x=Embarked))+
  labs(title='Embarked-Fare箱线图',y='票价',x='登船港口')+
  theme(plot.title = element_text(hjust = 0.5))
e2 <- ggplot(trainx)+
  geom_count(aes(x=Pclass,y=Embarked))+
  labs(title='Embarked-Pclass瓦图',x='乘客等级',y='登船港口')+
  theme(plot.title = element_text(hjust = 0.5))
grid.arrange(e1,e2,ncol=2)

mc <- trainx%>%filter(Pclass==1)%>%group_by(Embarked)%>%summarise(mean(Fare),median(Fare),n())
knitr::kable(mc,caption='登船港口等级统计',align='c')
trainx$Embarked[which(is.na(trainx$Embarked))] = 'C'
md.pattern(trainx)
#fare(test)
knitr::kable(filter(test,is.na(Fare)),caption='票价缺失统计',align='c')
mc2 <- test%>%filter(Pclass==3)%>%group_by(Embarked)%>%summarise(mean(Fare,na.rm=T),median(Fare,na.rm=T),n())
knitr::kable(mc2,caption='登船港口票价统计',align='c')
test$Fare[which(is.na(test$Fare))] =as.numeric(13.9)
md.pattern(test)

#变量因子化
trainx$Pclass <- as.numeric(trainx$Pclass)
trainx$Sex <- as.factor(trainx$Sex)
trainx$Embarked <- as.factor(trainx$Embarked)
trainx$Survived <- as.factor(trainx$Survived)
#整体存活率
sum(trainx$Survived==1)/length(trainx)
ggplot(trainx)+
  geom_bar(aes(x=Survived,fill=Survived))+
  labs(title='总体生存率柱形图',x='生存情况',y='生还人数')+
  theme(plot.title = element_text(hjust = 0.5))

r1 <- ggplot(trainx)+
  geom_bar(aes(x=Survived,fill=Pclass),position = 'dodge')+
  labs(title='Pclass柱形图',x='生存情况',y='生还人数')+
  theme(plot.title = element_text(hjust = 0.5))
r2 <- ggplot(trainx)+
  geom_bar(aes(x=Survived,fill=Sex),position = 'dodge')+
  labs(title='Sex柱形图',x='生存情况',y='生还人数')+
  theme(plot.title = element_text(hjust = 0.5))
r3 <- ggplot(trainx)+
  geom_bar(aes(x=Survived,fill=Sex),position = 'dodge')+
  facet_wrap(~Pclass)+
  labs(title='Pclass-Sex柱形图',x='生存情况',y='生还人数')+
  theme(plot.title = element_text(hjust = 0.5))
r4 <- ggplot(trainx)+
  geom_bar(aes(x=Survived,fill=Embarked),position = 'dodge')+
  labs(title='Embarked柱形图',x='生存情况',y='生还人数')+
  theme(plot.title = element_text(hjust = 0.5))
grid.arrange(r1,r2,r3,r4,ncol=2,nrow=2)

trainx%>%
  count(Survived,Embarked)%>%
  ggplot()+
  geom_tile(aes(x=Survived,y=Embarked,fill=n))+
  labs(title='Embarked-Survived瓦图',x='生存情况',y='登船港口')+
  theme(plot.title = element_text(hjust = 0.5))

f1 <- ggplot(trainx)+
  geom_freqpoly(aes(x=Age,color=Survived))+
  labs(title='年龄分布图',x='年龄',y='生还人数')+
  theme(plot.title = element_text(hjust = 0.5))
f2 <- ggplot(trainx)+
  geom_boxplot(aes(y=Age,x=Survived))+
  labs(title='生存年龄箱线图',y='年龄',x='生存情况')+
  theme(plot.title = element_text(hjust = 0.5))
f3 <- ggplot(trainx)+
  geom_freqpoly(aes(x=Fare,color=Survived))+
  labs(title='票价分布图',x='票价',y='生还人数')+
  theme(plot.title = element_text(hjust = 0.5))
f4 <- ggplot(trainx)+
  geom_boxplot(aes(y=Fare,x=Survived))+
  labs(title='生存票价箱线图',y='票价',x='生存情况')+
  theme(plot.title = element_text(hjust = 0.5))
grid.arrange(f1,f2,f3,f4,ncol=2,nrow=2)

fb1 <- ggplot(trainx)+
  geom_bar(aes(x=Parch,fill=Survived),position = 'dodge')+
  labs(title='Parch分布图',x='Parch',y='生还人数')+
  theme(plot.title = element_text(hjust = 0.5))
fb2 <- ggplot(trainx)+
  geom_bar(aes(x=SibSp,fill=Survived),position = 'dodge')+
  labs(title='SibSp分布图',x='SibSp',y='生还人数')+
  theme(plot.title = element_text(hjust = 0.5))
grid.arrange(fb1,fb2,ncol=2)

set.seed(2)
index <- sample(2,nrow(trainx),replace = TRUE,prob = c(0.7,0.3))
trainset <- trainx[index == 1,]
testset <- trainx[index == 2,]
#train
logic_mod <- glm(Survived ~.,family=binomial(link='logit'),data=trainset)
grid_logic <- data_grid(data=trainset,trainset[,1:7])
pre_logic <- grid_logic%>%
  add_predictions(logic_mod)%>%
  mutate(
    tpred=ifelse(pred>0.4275,1,0),
    true=trainset$Survived
    )
logic_accuracy <- 1-mean(pre_logic$tpred != trainset$Survived)
#test
ttpre_logic <- predict(logic_mod,test)
summary(logic_mod)
library(ROCR)
pr = prediction(pre_logic$pred, trainset$Survived)
prf = performance(pr, measure = "tpr", x.measure = "fpr")
plot(prf,main='logic-train-ROC曲线图')
#train-test
tpre_logic <- data_grid(data=testset,testset[,1:7])%>%
  add_predictions(logic_mod)%>%
  mutate(
    tpred=ifelse(pred>0.4275,1,0),
    true=testset$Survived
    )
logic_accuracy2 <- 1-mean(tpre_logic$tpred != testset$Survived)
prt = prediction(tpre_logic$pred, testset$Survived)
prft = performance(prt, measure = "tpr", x.measure = "fpr")
plot(prft,main='logic-test-ROC曲线图')

knitr::kable(cbind(logic_accuracy,logic_accuracy2),caption='logic模型准确率',align='c')
knitr::kable(table(pre_logic$tpred,trainset$Survived),caption='logic预测情况1',align='c')
knitr::kable(table(tpre_logic$tpred,testset$Survived),caption='logic预测情况2',align='c')
library(randomForest)
rf_mod <- randomForest(Survived~.,data = trainset, importance = T)
pre_rf <- data_grid(data=trainset,trainset[,1:7])%>%
  add_predictions(rf_mod)%>%
  mutate(
    true=trainset$Survived
  )
rf_mod
plot(rf_mod)
importance(rf_mod)
#train
rf_accuracy <- 1-mean(pre_rf$pred != trainset$Survived)
#train-test
tpre_rf <- data_grid(data=testset,testset[,1:7])%>%
  add_predictions(rf_mod)%>%
  mutate(
    true=testset$Survived
  )
rf_accuracy2 <- 1-mean(tpre_rf$pred != testset$Survived)
knitr::kable(cbind(rf_accuracy,rf_accuracy2),caption='rf模型准确率',align='c')
knitr::kable(table(pre_rf$pred,trainset$Survived),caption='rf预测情况1',align='c')
knitr::kable(table(tpre_rf$pred,testset$Survived),caption='rf预测情况2',align='c')

泰坦尼克号生存分析

李娅萍学号:218025200065

1. 选题介绍

2. 模型介绍

2.1 logistic回归

2.2 随机森林

3 数据预处理

3.1 数据说明

3.2 缺失情况

3.3 缺失处理

3.3.1 Age年龄缺失值处理

3.3.2 Embarked缺失值处理

4. 描述性统计

4.1 生存率分析

4.2 离散变量

4.3 连续变量

5. 模型建立

5.1 logistic模型

5.2 随机森林

6. 总结

附录

泰坦尼克号生存分析

李娅萍 学号:218025200065

1. 选题介绍

2. 模型介绍

2.1 logistic回归

2.2 随机森林

3 数据预处理

3.1 数据说明

3.2 缺失情况

3.3 缺失处理

3.3.1 Age年龄缺失值处理

3.3.2 Embarked缺失值处理

4. 描述性统计

4.1 生存率分析

4.2 离散变量

4.3 连续变量

5. 模型建立

5.1 logistic模型

5.2 随机森林

6. 总结

附录

李娅萍学号:218025200065