The Titanic Dataset 数据挖掘

简介

Titanic，就是当年第一航行便失事的超级大船——泰坦尼克号，大家可能对 Rose 与 Jack 的爱情故事念念不忘，但同时 Titanic 留下的乘客数据也是数据分析一笔宝贵的财富，很多新人都会拿这套数据集练练手。本想做一下 Titanic 的数据分析练练手，然后上 kaggle 一搜，一堆 Titanic 数据分析的 Paper，于是决定偷个懒，找了一篇 Fork 数最高的 Paper 翻译一下，当做练习自己的英语：）

附原文目录：

1 Introduction
- 1.1 Load and check data
2 Feature Engineering
- 2.1 What’s in a name?
- 2.2 Do families sink or swim together?
- 2.3 Treat a few more variables …
3 Missingness
- 3.1 Sensible value imputation
- 3.2 Predictive imputation
- 3.3 Feature Engineering: Round 2
4 Prediction
- 4.1 Split into training & test sets
- 4.2 Building the model
- 4.3 Variable importance
- 4.4 Prediction!
5 Conclusion

此文一共有三部分组成：

特征工程（Feature engineering）
缺失值填补（Missing value imputation）
预测（Prediction）（译者注：这点最重要！）

装载数据并检验

# Load packages
library('ggplot2') # 数据可视化
library('ggthemes') # 数据可视化
library('scales') # 数据可视化
library('dplyr') # 数据操作
library('mice') # 缺失值处理
library('randomForest') # 随机森林分类算法

以上的包安装完之后，可以开始读入数据了，顺便再简单看一下数据的内容。

train <- read.csv('../input/train.csv', stringsAsFactors = F)
test  <- read.csv('../input/test.csv', stringsAsFactors = F)

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" ...

通过上面打印的数据，可以看见变量名、变量的类型、以及他们开头一小部分的值。还有，通过第一行打印的内容我们也可以知道，我们操作的数据集总共有1309个个案与12个变量。因为上面打印的变量是缩写，表意不是很明确，所以这里我做一个表格，让我们的变量语义更加清晰。

特征工程

通过姓名我们可以知道什么？

第一个最值得关注的变量就是乘客姓名（Name），我们可以将其分解成额外的、有意义的变量，而这些变量恰恰有助于我们接下来的预测。比如，一个人的姓氏代表着他的家庭。现在，开始动手分析吧！

# 使用正则从乘客姓名中匹配它名字中的title部分
full$Title <- gsub('(.*, )|(\\..*)', '', full$Name)

# Show title counts by sex
table(full$Sex, full$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

# Titles with very low cell counts to be combined to "rare" level
rare_title <- c('Dona', 'Lady', 'the Countess','Capt', 'Col', 'Don', 
                'Dr', 'Major', 'Rev', 'Sir', 'Jonkheer')

# Also reassign mlle, ms, and mme accordingly
full$Title[full$Title == 'Mlle']        <- 'Miss' 
full$Title[full$Title == 'Ms']          <- 'Miss'
full$Title[full$Title == 'Mme']         <- 'Mrs' 
full$Title[full$Title %in% rare_title]  <- 'Rare Title'

# Show title counts by sex again
table(full$Sex, full$Title)

##         
##          Master Miss  Mr Mrs Rare Title
##   female      0  264   0 198          4
##   male       61    0 757   0         25

# Finally, grab surname from passenger name
full$Surname <- sapply(full$Name,  
                      function(x) strsplit(x, split = '[,.]')[[1]][1])

# 译者注：上面这行代码其实可以简化成下面这行
# full$Surname <- gsub(',.*',"",full$Name)

# count the unique surnames
nlevels(factor(full$Surname))

## [1] 875

我们得到了875个不同的姓名，其实基于这些姓名，我们可以推断出他们的种族关系，但是基于篇幅关系，这里就不做这项研究了。

家庭成员是一起溺死或幸存的吗？

既然我们已经把乘客的名字分割成了一些新的变量，我们可以进一步研究并做出一些新的家庭变量。首先，我们将根据兄弟姐妹/配偶的数量或者儿童/家长的人数（有的人可能拥有不止一个配偶）来创建一个家庭规模的变量（family size）。

# Create a family size variable including the passenger themselves
full$Fsize <- full$SibSp + full$Parch + 1

# Create a family variable 
full$Family <- paste(full$Surname, full$Fsize, sep='_')

family size这个变量究竟长什么样子？为了更好地理解它与幸存率的关系，先把他丢进训练数据里看看。

# Use ggplot2 to visualize the relationship between family size & survival
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以上，幸存率几乎为0。现在把family size这个变量分解成三个层次，创建一个离散的家庭规模变量。

# Discretize family size
full$FsizeD[full$Fsize == 1] <- 'singleton'
full$FsizeD[full$Fsize < 5 & full$Fsize > 1] <- 'small'
full$FsizeD[full$Fsize > 4] <- 'large'

# Show family size by survival using a mosaic plot
mosaicplot(table(full$FsizeD, full$Survived), main='Family Size by Survival', shade=TRUE)

可以看到，家庭规模越小，幸存率越高。

姓名（Name）这个变量基本已经挖掘完了。除了姓名，其实年龄变量（age）也有很大的挖掘价值。但是它有263条缺失值，我们得先处理它。

再处理几个变量……

在客舱变量（Cabin）里可能还有一些有用的信息，包括它们所在的甲板，让我们先看一看。

# This variable appears to have a lot of missing values
full$Cabin[1:28]

##  [1] ""            "C85"         ""            "C123"        ""           
##  [6] ""            "E46"         ""            ""            ""           
## [11] "G6"          "C103"        ""            ""            ""           
## [16] ""            ""            ""            ""            ""           
## [21] ""            "D56"         ""            "A6"          ""           
## [26] ""            ""            "C23 C25 C27"

# The first character is the deck. For example:
strsplit(full$Cabin[2], NULL)[[1]]

## [1] "C" "8" "5"

# Create a Deck variable. Get passenger deck A - F:
full$Deck<-factor(sapply(full$Cabin, function(x) strsplit(x, NULL)[[1]][1]))

其实到了这里，还有更多可以做的，比如看一下客舱的位置排列信息（e.g., row 28: “C23 C25 C27”），但是鉴于篇幅问题，我们这里的分析到此为止。

缺失值处理

现在，我们准备开始研究那些缺失的数据，并想办法纠正它，这其中有许多不同的修补方式。鉴于数据集是小尺寸的，我们不可以选择直接删除包含缺失值的整个个案（行）或变量（列）。这样子，我们留下的选项，要么是用一个合理的数据去填补缺失值，例如使用均值、中位数或众数去修补缺失值；要么，去预测。我们将混合使用这两种方法，并且依靠一些数据可视化来指导我们去决定。

合理值修补

# Passengers 62 and 830 are missing Embarkment
full[c(62, 830), 'Embarked']

## [1] "" ""

好的，62号乘客与830号乘客的登船地点的数据缺失了。那么，现在我们来推断出他们登船地点的缺失值。因为，基于目前的数据，我们可以想象乘客类别（Pclass）与票价（Fare）这两个变量具有某种关联性。

# Passenger 62 paid
full[c(62, 830), 'Fare'][[1]][1]
# Passenger 830 paid
full[c(62, 830), 'Fare'][[1]][2]

## [1] 80
## [2] NA

# Passenger 62 Pclass
full[c(62, 830), 'Pclass'][[1]][1]
# Passenger 830 Pclass
full[c(62, 830), 'Pclass'][[1]][2]

## [1] 1
## [2] NA

我们可以看到，他们分别支付了 $80 与 $NA，他们的类别分别是 1 与 NA。那么，他们到底是在哪里登船的呢？

# Get rid of our missing passenger IDs
embark_fare <- full %>%
  filter(PassengerId != 62 & PassengerId != 830)

# Use ggplot2 to visualize 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()

就是这样！一个头等舱的乘客离开 Charbourg 的票价中值（C）正好对上了支付了 $80 缺缺少登船地的乘客。因此，可以安全的将 NA 值替换成 C 。

# Since their fare was $80 for 1st class, they most likely embarked from 'C'
full$Embarked[c(62, 830)] <- 'C'

我们已经修复了这里的 NA 值。继续看，第1044排的乘客是否具有缺失值。

# Show row 1044
full[1044, ]

##      PassengerId Survived Pclass               Name  Sex  Age SibSp Parch
## 1044        1044       NA      3 Storey, Mr. Thomas male 60.5     0     0
##      Ticket Fare Cabin Embarked Title Surname Fsize   Family    FsizeD
## 1044   3701   NA              S    Mr  Storey     1 Storey_1 singleton
##      Deck
## 1044 <NA>

这是从 Southampton 出发的三等舱乘客，他们的票价Fare其实隐藏着他们的舱位价格和开船地（n = 494）。

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()

通过这个数据可视化图表，我们可以很合理地将 NA 值的票价替代为 Fare的中值 $8.05。

# Replace missing fare value with median fare for class/embarkment
full$Fare[1044] <- median(full[full$Pclass == '3' & full$Embarked == 'S', ]$Fare, na.rm = TRUE)

预测值修补

最后，正如我们前面提到的，在我们的数据中有很多缺失的年龄值（Age）。我们要在填充缺失值时融入更多的幻想。为什么？因为我们可以基于其他已有变量的基础上来创建一个预测年龄的模型。

# Show number of missing Age values
sum(is.na(full$Age))

## [1] 263

我们可以使用 rpart（递归分割回归）来预测缺失的年龄数据，但现在我们准备使用的 mice 包是一个不同的东西。首先，我们先分解因子变量，然后使用mice修补数据。

# Make variables factors into factors
factor_vars <- c('PassengerId','Pclass','Sex','Embarked',
                 'Title','Surname','Family','FsizeD')

full[factor_vars] <- lapply(full[factor_vars], function(x) as.factor(x))

# Set a random seed
set.seed(129)

# Perform mice imputation, excluding certain less-than-useful variables:
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

# Save the complete output 
mice_output <- complete(mice_mod)

让我们比较我们mice后得到乘客年龄数据分布与原来年龄数据是否完全一致，以确保没有错误。

# Plot age distributions
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后的结果替换掉原始数据。

# Replace Age variable from the mice model.
full$Age <- mice_output$Age

# Show new number of missing Age values
sum(is.na(full$Age))

## [1] 0

其实现在，我们已经得到了一个完整的Age变量。接下来我会做一些润色工作，对Age做一些特征工程。

二轮特征工程

现在我们知道了每个人的年龄，我们可以创建一些新的年龄相关变量：Child 和 Mother。一个 child 只会是18岁以下的人，一个 mother 则必须满足 4 个条件：1）女性，2）超过18岁，3）至少有1个孩子，4）名字中没有冠以 Miss 的称呼。

# First we'll look at the relationship between age & survival
ggplot(full[1:891,], aes(Age, fill = factor(Survived))) + 
  geom_histogram() + 
  # I include Sex since we know (a priori) it's a significant predictor
  facet_grid(.~Sex) + 
  theme_few()

# Create the column child, and indicate whether child or adult
full$Child[full$Age < 18] <- 'Child'
full$Child[full$Age >= 18] <- 'Adult'

# Show counts
table(full$Child, full$Survived)

##        
##           0   1
##   Adult 484 274
##   Child  65  68

看起来像孩子的幸存率确实高一点（达到了50%），成人都有义务来保证孩子们不受到伤害。现在，我们将通过创建 Mother 变量来完成我们的特性工程。也许我们可以希望母亲们更有可能在泰坦尼克号上幸存下来。

# Adding Mother variable
full$Mother <- 'Not Mother'
full$Mother[full$Sex == 'female' & full$Parch > 0 & full$Age > 18 & full$Title != 'Miss'] <- 'Mother'

# Show counts
table(full$Mother, full$Survived)

##             
##                0   1
##   Mother      16  39
##   Not Mother 533 303

# Finish by factorizing our two new factor variables
full$Child  <- factor(full$Child)
full$Mother <- factor(full$Mother)

所有变量都应该被我们照顾到，所以不应该再有缺失值。因此，现在再加倍检查一次：

md.pattern(full)

## Warning in data.matrix(x): NAs introduced by coercion

## Warning in data.matrix(x): NAs introduced by coercion

## Warning in data.matrix(x): NAs introduced by coercion

##     PassengerId Pclass Sex Age SibSp Parch Fare Embarked Title Surname
## 150           1      1   1   1     1     1    1        1     1       1
##  61           1      1   1   1     1     1    1        1     1       1
##  54           1      1   1   1     1     1    1        1     1       1
## 511           1      1   1   1     1     1    1        1     1       1
##  30           1      1   1   1     1     1    1        1     1       1
## 235           1      1   1   1     1     1    1        1     1       1
## 176           1      1   1   1     1     1    1        1     1       1
##  92           1      1   1   1     1     1    1        1     1       1
##               0      0   0   0     0     0    0        0     0       0
##     Fsize Family FsizeD Child Mother Ticket Survived Deck Name Cabin     
## 150     1      1      1     1      1      1        1    1    0     0    2
##  61     1      1      1     1      1      1        0    1    0     0    3
##  54     1      1      1     1      1      0        1    1    0     0    3
## 511     1      1      1     1      1      1        1    0    0     0    3
##  30     1      1      1     1      1      0        0    1    0     0    4
## 235     1      1      1     1      1      1        0    0    0     0    4
## 176     1      1      1     1      1      0        1    0    0     0    4
##  92     1      1      1     1      1      0        0    0    0     0    5
##         0      0      0     0      0    352      418 1014 1309  1309 4402

哇塞！我们终于完成了修补缺失值的工作，泰坦尼克号的数据终于完整了！我们还成功地创建了一些新的变量，这将有助于我们建立一个模型，可靠地预测生存。

预测

最后，在我们精心处理过的泰坦尼克号数据集之上，已经准备好了幸存预测。为此，我们将使用随机森林分类算法。

分割训练数据与测试数据

第一步，把原始数据集分割成训练数据和测试数据两部分。

# Split the data back into a train set and a test set
train <- full[1:891,]
test <- full[892:1309,]

建模

我们将在训练数据集上使用randomForest（随机森林算法）建立模型。

# Set a random seed
set.seed(754)

# Build the model (note: not all possible variables are used)
rf_model <- randomForest(factor(Survived) ~ Pclass + Sex + Age + SibSp + Parch + 
                                            Fare + Embarked + Title + 
                                            FsizeD + Child + Mother,
                                            data = train)

# Show model error
plot(rf_model, ylim=c(0,0.36))
legend('topright', colnames(rf_model$err.rate), col=1:3, fill=1:3)

黑线表明我们数据整体的错误率低于20%。红线和绿线则表明了“死亡”和“幸存”的错误率。我们可以看到，相对于预测幸存，现在我们可以更准确预测死亡。我想知道这对我来说意味着什么？

变量的相对重要性

让我们看下变量的相对重要性。

# Get importance
importance    <- importance(rf_model)
varImportance <- data.frame(Variables = row.names(importance), 
                            Importance = round(importance[ ,'MeanDecreaseGini'],2))

# Create a rank variable based on importance
rankImportance <- varImportance %>%
  mutate(Rank = paste0('#',dense_rank(desc(Importance))))

# Use ggplot2 to visualize the relative importance of variables
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() + 
  theme_few()

哇，很高兴看到我们新建的Title变量具有最高的相对重要性。但更惊奇的是，Pclass变量只排到第5，但这也许只是作为一个孩子因为看电影泰坦尼克号的次数太多而产生的偏见。

预测

至此，我们已经准备好了最后一步——做出预测！但我们仍可以遍历前面的步骤进行调整以去适应数据采用不同的模型或使用变量的不同组合，来达到更好地预测。

# Predict using the test set
prediction <- predict(rf_model, test)

# Save the solution to a dataframe with two columns: PassengerId and Survived (prediction)
solution <- data.frame(PassengerID = test$PassengerId, Survived = prediction)

# Write the solution to file
write.csv(solution, file = 'rf_mod_Solution.csv', row.names = F)

注：本文由 Airing 翻译自 Megan L. Risdal. Exploring the Titanic Dataset https://www.kaggle.com/mrisdal/titanic/exploring-survival-on-the-titanic

本文由 Airing 创作，采用知识共享署名4.0 国际许可协议进行许可
本站文章除注明转载/出处外，均为本站原创或翻译，转载前请务必署名
最后编辑时间为: Dec 8, 2018 at 02:45 am