Kaggle机器学习实战（5）——再探泰坦尼克

机器学习
, ,

之前已经看过一篇Megan Risdal的kernel，她使用R语言、随机森林进行了建模。那时我是个纯新手，连依葫芦画瓢都不会。现在通过周志华老师的《机器学习》、吴恩达老师的coursera课程、查阅相关文档掌握的调包侠知识，总算是可以把这个经典入门案例拿过来练练手。

关于EDA部分，由于和其他人的步骤大同小异，也没有新的发现，就略过不谈。在模型方面，首先是调用sklearn包中一些常见的分类器查看其分类效果，其次是使用Stacking方法进行集成，最后的结果会提交到kaggle上。

本文参考了4篇点赞数量较高的kernel：

读取数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# 调包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import re
import warnings
warnings.filterwarnings("ignore")

from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LogisticRegressionCV
from sklearn.neural_network import MLPClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.svm import SVC
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import BernoulliNB
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.ensemble import AdaBoostClassifier
from sklearn.ensemble import GradientBoostingClassifier
from xgboost import XGBClassifier
from lightgbm import LGBMClassifier

对数据进行简单预处理。

1
2
3
4
5
6
7
8
9
10
train = pd.read_csv("train.csv")
test = pd.read_csv("test.csv")

y = train["Survived"]
train = train.drop(["Survived"], axis=1)
all_data = pd.concat([train, test])
PassengerId = test["PassengerId"]
all_data = all_data.drop(["PassengerId"], axis=1)
print('训练集和测试集共有{0}个样本和{1}个特征。'.format(all_data.shape[0], all_data.shape[1]))
all_data.head()

缺失值处理

查看缺失值

使用isnull()方法查看缺失值，并按降序排列。

1
all_data.isnull().sum().sort_values(ascending=False)

可以看到：Cabin具有1014个缺失值；Age具有263个缺失值；Embarked具有2个缺失值；Fare具有1个缺失值。我们将分别用不同的手段处理它们。

处理’Cabin’的缺失值

Cabin（客舱号）缺失值达到77%，我们可以考虑舍弃掉这个变量。仔细分析，可以理解成：不为空的表示该乘客有客舱，为空的表示没有客舱。因此，把这个变量转换为一个二值变量，即“有客舱”，就实现了对数据的利用。或许提取客舱号第一位的字母能更进一步的利用此条数据，但从他人的kernel中看出似乎没有太大的效果。

1
all_data["Cabin"] = all_data["Cabin"].apply(lambda x: 0 if type(x) == float else 1)

处理’Age’的缺失值

Age（年龄）缺失值达到了20%，已经超过可以舍弃的阈值了。但年龄绝对是一个非常重要的变量，因此要好好地去填充缺失值。一种方法是使用随机数填充，随机数的上下限是（均值±标准差）。另一种是使用其他变量（如票价、船舱等级）来“猜”年龄。

这里参照Anisotropic和Yassine Ghouzam，结合了两种方法。Yassine发现年龄与船舱等级、亲属数量有较强的相关性，因此找缺失年龄的样本是否具有与其同样Pclass、SibSp和Parch的样本，有的话就取这些样本年龄的中值，否则取所有样本年龄的中值。我在处理没有类似的样本时，取（均值±标准差）中的随机整数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 获取空值的index
index_NaN_age = list(all_data["Age"][all_data["Age"].isnull()].index)

# 默认用随机数填充
age_mean = all_data["Age"].mean()
age_std = all_data["Age"].std()
age_random_list = np.random.randint(age_mean-age_std, age_mean+age_std, size=len(index_NaN_age))
all_data["Age"][all_data["Age"].isnull()] = age_random_list

# 当存在类似情况时，用类似样本的中位数填充空值
for i in index_NaN_age:
    same_samples_age = all_data["Age"][((all_data["SibSp"] == all_data["SibSp"].iloc[i]) & (all_data["Parch"] == all_data["Parch"].iloc[i]) & (all_data["Pclass"] == all_data["Pclass"].iloc[i]))].median()
    if not np.isnan(same_samples_age):
        all_data["Age"].iloc[i] = same_samples_age

# 检查是否完全填充
assert all_data["Age"].isnull().sum() == 0

处理’Embarked’和’Fare’的缺失值

这两项缺失值很少，可以逐个填充。

先来看Embarked（登船港口），只有两个缺失值。将他们挑出：

1
2
# 筛出两条缺失值
all_data[all_data["Embarked"].isnull()]

大多数人使用众数S填充，因为从S登船的占对大多数；而Megan Risdal则通过观察，选择了用C填充。

可以发现：她们票号一致，所以必然是同一个港口上船的。观察同为头等舱的三个港口上船乘客的票价的均值、中位数：

1
all_data[all_data["Pclass"]==1].groupby("Embarked")["Fare"].describe()

其实无论是用C还是S填充都是合理的。80刀的票价，在C港口登船的头等舱的所有票价中位于中位数附近，且与在S港口登船的头等舱票价均值相近（接近75%四分位点）。由于S是众数（无论实在样本里，还是在众多解答者的选择中），最终决定用S填充。

1
all_data["Embarked"] = all_data["Embarked"].fillna("S")

对于存在唯一缺失值的Fare，我先找了找存不存在与其Ticket一样的样本，没有找到。于是就找同为三等舱、S港口登船的样本的中位数来填充。

1
2
3
4
5
fare_med = all_data["Fare"][(all_data["Embarked"]=='S') & (all_data["Pclass"]==3)].median()
all_data["Fare"] = all_data["Fare"].fillna(fare_med)

# 确认所有缺失值处理完毕
all_data.isnull().sum()

特征工程

处理’Name’

毫不奇怪地，大多数高分kernel都提取了姓名中的称谓，比如Mr、Miss等，其他一些稀少的称谓归为Rare一类，将称谓作为一个新的变量使用，并抛弃原“姓名”这个没有用的变量。这里需要使用到正则表达式，详细用法参见python re的官方文档。

数据的Name形如Braund, Mr. Owen Harris，我们要提取的Mr.这些称谓冠词的共同点是末尾有一个缩写点，所以用re.search()查找点号前的子串。

1
2
3
4
5
6
7
8
9
10
def getTitle(name):
    title_search = re.search('([A-Za-z]+)\.', name)
    if title_search:
        return title_search.group(1)
    return ""

all_data["Title"] = all_data["Name"].apply(getTitle)

# 查看各种称谓出现的次数
all_data["Title"].value_counts()

我打算只保留”Mr”、”Miss”（同义词”Mlle”、”Ms” ）、”Mrs”（同义词”Mme”）、”Master”这四个比较常见的称谓，其他归入”Rare”。

1
2
3
4
5
6
7
8
rare_title = ["Dr","Rev","Col","Major","Don","Jonkheer","Sir","Lady","Countess","Capt","Dona"]

all_data["Title"].replace(rare_title, "Rare", inplace=True)
all_data["Title"].replace(["Mlle","Ms"], "Miss", inplace=True)
all_data["Title"].replace(["Mme"], "Mrs", inplace=True)

# 检查是否完成合并
all_data["Title"].value_counts()

对于Name变量的操作就到此为止。

处理’SibSp’和’Parch’

我所阅读的3篇kernel都是把这两项相加（还要把自己算进去），得到familySize（家庭人数）这一新特征，其中如果该新特征等于1，就表示是他是孤身一人登船，他遇难的概率会大大增加。这就是一个对预测生还率有用的变量。

参考Anisotropic和Manav Sehgal的解答，我会创建一个名为“孤身一人”的二值变量，其它变量均舍去。Yassine Ghouzam还考虑了人数大于4的大家庭，以后也可以考虑一下。

1
2
3
4
5
6
all_data["FamilySize"] = all_data["SibSp"] + all_data["Parch"] + 1
all_data["IsAlone"] = 0
all_data.loc[all_data["FamilySize"] == 1, "IsAlone"] = 1

# 查看是否为一人的人数
all_data["IsAlone"].value_counts()

数据转化

通过阅读了这几篇kernel，我发现部分解答对数据的处理存在如下共同点：

  • Age、Fare分段：这两项是连续的数值型数据，由于我们采用的集成模型大多数是基于决策树的，决策树模型本质上是将它当离散值处理的，因此手动将其分为我们想要的几段，可能会获得更好的性能。
  • 将字符型变量标签编码，以及虚拟变量，由于模型在处理字符型变量时会自动转换为哑变量，这一步似乎不做也无妨？
  • 由于最后选择的变量全是代表类别的离散型，无需做特征缩放。

对连续值分段，pandas有cut和qcut两种方式可用，前者依据值本身的大小来确定切割点，后者依据出现的频率来确定切割点。参考Anisotropic和Manav Sehgal的解答，对Age采用cut划分为5个分段，对Fare采用qcut划分为4个分段。这一步帮助我们确定分割点，下一步是手动创建分割后的标签变量。

1
2
3
4
5
all_data["AgeCut"] = pd.cut(all_data["Age"], 5)
all_data["AgeCut"].value_counts()

all_data["FareCut"] = pd.qcut(all_data["Fare"], 4)
all_data["FareCut"].value_counts()

年龄以16岁为一个阶梯，票价的分割点则是[7.896, 14.454, 31.275]。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 年龄分段
all_data.loc[all_data["Age"]<16, "AgeGroup"] = 0
all_data.loc[(all_data["Age"]>=16) & (all_data["Age"]<32), "AgeGroup"] = 1
all_data.loc[(all_data["Age"]>=32) & (all_data["Age"]<48), "AgeGroup"] = 2
all_data.loc[(all_data["Age"]>=48) & (all_data["Age"]<64), "AgeGroup"] = 3
all_data.loc[all_data["Age"]>=64, "AgeGroup"] = 4

# 票价分段
all_data.loc[all_data["Fare"]<7.896, "FareGroup"] = 0
all_data.loc[(all_data["Fare"]>=7.896) & (all_data["Fare"]<14.454), "FareGroup"] = 1
all_data.loc[(all_data["Fare"]>=14.454) & (all_data["Fare"]<31.275), "FareGroup"] = 2
all_data.loc[all_data["Fare"]>=31.275, "FareGroup"] = 3

# 取整
all_data["AgeGroup"] = all_data["AgeGroup"].astype("int")
all_data["FareGroup"] = all_data["FareGroup"].astype("int")

接着进行特征选择，把没有用的特征抛弃掉。

1
2
3
4
5
drop_features = ["Name","Age","SibSp","Parch","Ticket","Fare","FamilySize","AgeCut","FareCut"]
all_data.drop(drop_features, axis=1, inplace=True)

# 查看结果
all_data.head()

剩下来的变量均为类别型变量，且不存在明确的顺序关系，因此将其虚拟化。由于get_dummies()仅将字符型变量虚拟化，因此一些变量需要手动转换。

1
2
3
4
5
6
7
8
9
10
11
12
13
# 性别转换为二值变量
all_data["Sex"] = all_data["Sex"].map({"female": 0, "male": 1}).astype(int)

# 整型转换为字符型
all_data["Pclass"] = all_data["Pclass"].astype("str")
all_data["AgeGroup"] = all_data["AgeGroup"].astype("str")
all_data["FareGroup"] = all_data["FareGroup"].astype("str")

# 虚拟变量
all_data = pd.get_dummies(all_data)

# 查看虚拟化后的变量名
all_data.columns

至此我们已经完成了特征工程，我们把训练集和测试集分开，准备进行下一步。

1
2
3
# 分离训练集和测试集
X_train = all_data[:y.shape[0]]
X_test = all_data[y.shape[0]:]

运用单个模型

作为一次机器学习训练，自然是把所有能用的分类模型都拿过来试一试。除了scikit-learn自带的几个模型之外，还使用了XGBoost和LightGBM这两个比赛中常用的集成模型。把每个模型都调用来比较一下，找到最优的分类器，看看其精度如何，最后和Stacking模型的精度比一比。

尝试使用的模型：

  • 逻辑回归
  • 神经网络
  • CART
  • SVM
  • KNN
  • 朴素贝叶斯
  • 随机森林
  • 极端随机树
  • AdaBoost
  • GBDT
  • XGBoost
  • LightGBM

先使用默认参数都跑一遍，用交叉验证查看准确率的均值和标准差。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
models = {}
models['LogisticRegression'] = LogisticRegressionCV()
models['NeuralNetwork'] = MLPClassifier()
models['CART'] = DecisionTreeClassifier()
models['SVM'] = SVC()
models['KNN'] = KNeighborsClassifier()
models['NaiveBayes'] = BernoulliNB()
models['RandomForest'] = RandomForestClassifier()
models['ExtraTree'] = ExtraTreesClassifier()
models['AdaBoost'] = AdaBoostClassifier()
models['GBDT'] = GradientBoostingClassifier()
models['XGBoost'] = XGBClassifier()
models['LightGBM'] = LGBMClassifier()

kf = KFold(10)

for model in models:
    cv_result = cross_val_score(models[model], X_train, y, cv=kf, scoring='accuracy')
    print('%s模型的交叉验证得分平均值%.2f%%，标准差%.2f%%。' % (model, cv_result.mean()*100, cv_result.std()*100))

这里没有设置随机数种子，因此每次跑模型，结果都是不一样的。

1
2
3
4
5
6
7
8
9
10
11
12
LogisticRegression模型的交叉验证得分平均值80.47%，标准差3.86%。
NeuralNetwork模型的交叉验证得分平均值80.81%，标准差4.46%。
CART模型的交叉验证得分平均值80.36%，标准差3.30%。
SVM模型的交叉验证得分平均值79.58%，标准差4.06%。
KNN模型的交叉验证得分平均值79.13%，标准差4.39%。
NaiveBayes模型的交叉验证得分平均值74.98%，标准差5.07%。
RandomForest模型的交叉验证得分平均值80.93%，标准差3.55%。
ExtraTree模型的交叉验证得分平均值80.59%，标准差3.24%。
AdaBoost模型的交叉验证得分平均值79.01%，标准差4.40%。
GBDT模型的交叉验证得分平均值81.71%，标准差3.35%。
XGBoost模型的交叉验证得分平均值82.27%，标准差3.63%。
LightGBM模型的交叉验证得分平均值81.15%，标准差3.81%。

可以看出大多数模型的准确率在80%附近，XGBoost效果最好，可以进一步调参优化；朴素贝叶斯较差，可以放弃。

下面使用网格搜索法对随机森林、极端随机树、AdaBoost、GBDT、XGBoost、LightGBM进行调参。主要调试n_estimators，对于Boosting类的模型，还会调试learnig_rate这个参数，其余保持默认设置。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
bagging_models = {'RandomForest': RandomForestClassifier(), 'ExtraTree': ExtraTreesClassifier()}
boosting_models = {'AdaBoost': AdaBoostClassifier(), 'GBDT': GradientBoostingClassifier(), 'XGBoost': XGBClassifier(), 'LightGBM': LGBMClassifier()}

bagging_params = {'n_estimators': [10, 50, 100, 200, 500, 800]}
boosting_params = {'n_estimators': [10, 50, 100, 200, 500, 800], 'learning_rate':[0.005, 0.01, 0.1, 1]}

kf = KFold(10)

for model in bagging_models:
    grid = GridSearchCV(estimator=bagging_models[model], param_grid=bagging_params, cv=kf, scoring='accuracy')
    grid_result = grid.fit(X_train, y)
    print('%s模型的最优参数是%s，得分%.2f%%。' % (model, grid_result.best_params_, grid_result.best_score_*100))

for model in boosting_models:
    grid = GridSearchCV(estimator=boosting_models[model], param_grid=boosting_params, cv=kf, scoring='accuracy')
    grid_result = grid.fit(X_train, y)
    print('%s模型的最优参数是%s，得分%.2f%%。' % (model, grid_result.best_params_, grid_result.best_score_*100))

一次跑下来的结果为：

1
2
3
4
5
6
RandomForest模型的最优参数是{'n_estimators': 50}，得分81.37%。
ExtraTree模型的最优参数是{'n_estimators': 200}，得分80.81%。
AdaBoost模型的最优参数是{'learning_rate': 0.01, 'n_estimators': 800}，得分80.25%。
GBDT模型的最优参数是{'learning_rate': 0.1, 'n_estimators': 50}，得分82.15%。
XGBoost模型的最优参数是{'learning_rate': 0.005, 'n_estimators': 800}，得分82.49%。
LightGBM模型的最优参数是{'learning_rate': 0.01, 'n_estimators': 500}，得分81.82%。

如果使用单一模型，则XGBoost效果最好。接下来，尝试使用Stacking模型。

运用Stacking模型

之前在房价预测中，接触到了Stacking这一kaggle神器，这次来亲手尝试一下，构造一个Stacking类，虽然代码是直接抄过来的。

原作者kernel地址：《Stacked Regressions : Top 4% on LeaderBoard》。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
from sklearn.base import BaseEstimator, TransformerMixin, ClassifierMixin, clone

class StackingAveragedModels(BaseEstimator, ClassifierMixin, TransformerMixin):
    def __init__(self, base_models, meta_model, n_folds=5):
        self.base_models = base_models
        self.meta_model = meta_model
        self.n_folds = n_folds

    # 定义fit方法
    def fit(self, X, y):
        self.base_models_ = [list() for x in self.base_models]
        self.meta_model_ = clone(self.meta_model)
        kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=1)
        out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))
        # 用K折后的数据训练每一个初级分类器
        for i, model in enumerate(self.base_models):
            for train_index, holdout_index in kfold.split(X, y):
                instance = clone(model)
                self.base_models_[i].append(instance)
                instance.fit(X.iloc[train_index], y.iloc[train_index])
                y_pred = instance.predict(X.iloc[holdout_index])
                out_of_fold_predictions[holdout_index, i] = y_pred
        # 使用次级分类器拟合初级分类器预测的结果
        self.meta_model_.fit(out_of_fold_predictions, y)
        return self

    # 定义predict方法
    def predict(self, X):
        meta_features = np.column_stack([np.column_stack([model.predict(X) for model in base_models]).mean(axis=1) for base_models in self.base_models_])
        return self.meta_model_.predict(meta_features)

Anisotropic使用的Stacking是以随机森林、极端随机树、AdaBoost、GBDT、SVM作为元分类器，XGBoost作为次级分类器。于是我这样做：随机森林、极端随机树、AdaBoost、GBDT、LightGBM作为元分类器，XGBoost作为次级分类器。

1
2
3
4
5
6
7
8
9
10
rf = RandomForestClassifier(n_estimators=50)
et = ExtraTreesClassifier(n_estimators=200)
ab = AdaBoostClassifier(n_estimators=800, learning_rate=0.01)
gb = GradientBoostingClassifier(n_estimators=50, learning_rate=0.1)
xg = XGBClassifier(n_estimators=800, learning_rate=0.005)
lg = LGBMClassifier(n_estimators=500, learning_rate=0.01)

stacked_averaged_models = StackingAveragedModels(base_models=(rf, et, ab, gb, lg), meta_model=xg)
score = cross_val_score(stacked_averaged_models, X_train, y, cv=kf, scoring='accuracy')
print('Stacking模型的交叉验证得分平均值%.2f%%，标准差%.2f%%。' % (score.mean()*100, score.std()*100))
1
Stacking模型的交叉验证得分平均值81.37%，标准差3.92%。

结果令人失望：验证集精度不如单一模型。无论如何，分别提交XGBoost和Stacking的结果看看测试集精度如何。

1
2
3
4
5
6
7
8
9
10
11
12
# XGBoost
xg = XGBClassifier(n_estimators=800, learning_rate=0.005)
xg.fit(X_train, y)
xg_pred = xg.predict(X_test)
XGBoostSubmission = pd.DataFrame({'PassengerId': PassengerId, 'Survived': xg_pred})
XGBoostSubmission.to_csv("XGBoostSubmission.csv", index=False)

# Stacking
stacked_averaged_models.fit(X_train, y)
stacking_pred = stacked_averaged_models.predict(X_test)
StackingSubmission = pd.DataFrame({'PassengerId': PassengerId, 'Survived': stacking_pred})
StackingSubmission.to_csv("StackingSubmission.csv", index=False)

Stacking模型在测试集上的得分为0.77990；XGBoost模型在测试集上的得分为0.79904，位列前14.9%。

总结

泰坦尼克号生还预测是一个非常适合入门的机器学习项目。通过这一次实战训练，有以下收获：

  • 掌握pandas的数据清洗和特征工程方法。
  • 熟悉sklearn的各种分类模型及其调参。
  • 了解Stacking，虽然这个任务中性能似乎不如XGBoost，但可能适合别的任务。
  • 在Kaggle阅读了越来越多的kernel，了解机器学习的前沿方法。