初心者がKaggle Titanicで上位1.5%(0.83732)以内に入るアプローチ解説_3

[前回](https://qiita.com/shiroino11111/items/fb6aa6b7dba2ddc3ce04)に続き、[Kaggle Titanic](https://www.kaggle.com/c/titanic/submissions)で上位1.5%(0.83732)へのアプローチを解説していきます。
使用するコードは[Github](https://github.com/shiroino11111/taitanic-0.83732-)のtitanic(0.83732)_3です。
[前回](https://qiita.com/shiroino11111/items/fb6aa6b7dba2ddc3ce04)に出した提出スコアから、0.83732まで向上させる解説を行います。
###1.必要なライブラリをインポートし、CSVを読み込む。

```py
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.ensemble import RandomForestRegressor
from sklearn.pipeline import Pipeline,make_pipeline
from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import SelectKBest
from sklearn import model_selection
from sklearn.model_selection import GridSearchCV
import warnings
warnings.filterwarnings('ignore')
```
```py
# CSVを読み込む
train= pd.read_csv("train.csv")
test= pd.read_csv("test.csv")

# データの統合
dataset = pd.concat([train, test], ignore_index = True)

# 提出用に
PassengerId = test['PassengerId']
```

```py
# 客室階層による生存率比較
dataset['Cabin'] = dataset['Cabin'].fillna('Unknown') # 客室データが欠損している場合はUnknownを代入
dataset['Deck']= dataset['Cabin'].str.get(0) #Cabin（部屋番号）の頭文字（０番目の文字）取得

# チケットの文字数による生存率比較
Ticket_Count = dict(dataset['Ticket'].value_counts()) # チケットの文字数でグループ分け
dataset['TicketGroup'] = dataset['Ticket'].apply(lambda x:Ticket_Count[x]) # グループの振り分け

# チケットの文字数での生存率が高いグループと低いグループの２つに分ける。
# 高ければ２，低ければ１を代入
def Ticket_Label(s):
    if (s >= 2) & (s <= 4): # 文字数での生存率が高いグループ
        return 2
    elif ((s > 4) & (s <= 8)) | (s == 1): # 文字数での生存率が低いグループ
        return 1
    elif (s > 8):
        return 0

dataset['TicketGroup'] = dataset['TicketGroup'].apply(Ticket_Label)
```
###2.敬称を利用
Kaggle上位者のコードを見ると、名前の敬称を利用することが高スコアのカギだと分かります。
敬称とは、Nameの中間当たりに含まれているMr、Mrs、Missなどのことです。Mrを使わずDr(医師)やRev(聖職者や牧師)といった職業が記載されていることもあります。
これらの情報を抽出してグルーピングします。

```py
# 'Honorifics'(敬称)による特徴別に分ける
dataset['Honorifics'] = dataset['Name'].apply(lambda x:x.split(',')[1].split('.')[0].strip()) #敬称(','と'.'の間の単語)を抽出

# 敬称をグループ分け
# 例：'Capt', 'Col', 'Major', 'Dr', 'Rev'は'Officer'とする
Honorifics_Dict = {}
Honorifics_Dict.update(dict.fromkeys(['Capt', 'Col', 'Major', 'Dr', 'Rev'], 'Officer'))
Honorifics_Dict.update(dict.fromkeys(['Don', 'Sir', 'the Countess', 'Dona', 'Lady'], 'Royalty'))
Honorifics_Dict.update(dict.fromkeys(['Mme', 'Ms', 'Mrs'], 'Mrs'))
Honorifics_Dict.update(dict.fromkeys(['Mlle', 'Miss'], 'Miss'))
Honorifics_Dict.update(dict.fromkeys(['Mr'], 'Mr'))
Honorifics_Dict.update(dict.fromkeys(['Master','Jonkheer'], 'Master'))
dataset['Honorifics'] = dataset['Honorifics'].map(Honorifics_Dict)
sns.barplot(x="Honorifics", y="Survived", data=dataset, palette='Set3')

"""敬称一覧
Mr：男 , Master：男の子, Jonkheer：オランダ貴族(男),
Mlle：マドモワゼル (フランス未婚女性), Miss：未婚女性、女の子, Mme：マダム(フランス既婚女性), Ms：女性(未婚・既婚問わず), Mrs：既婚女性, 
Don：男(スペイン), Sir：男(イギリス), the Countess：伯爵夫人, Dona：既婚女性(スペイン), Lady：既婚女性(イギリス),
Capt：船長, Col：大佐, Major：軍人, Dr：医者, Rev：聖職者や牧師
"""
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/16428cc8-80a6-1771-10f1-635704295ca1.png)
やはり成人男性は生存率が低く、女性や子供が高いです。が、今回、貴族などのRoyaltyグループが子供以上に高いことが発見できます。
この時代の貴族が優先して助かっているのは分かる気がします。
生存率に貴族であるかどうかをデータとして使えるのは武器になりそうです。

###3.欠損値の代入を見直し
[前回](https://qiita.com/shiroino11111/items/fb6aa6b7dba2ddc3ce04)、[前々回](https://qiita.com/shiroino11111/items/bc3889fa38ff32d46c13)は、とりあえず欠損値に中央値を入れました。
予測の精度を上げるにはこれらを見直します。
###3.1 'Age'(年齢)の欠損値の見直し
機械学習で予測したものを年齢の欠損値に代入
予測には先程の敬称(職業)データも使えそうです。（Drで５才といった予測を抑えられる)

```py
## 年齢の欠損値を予測して代入
# 年齢予測のために使う項目を抽出し、ダミー変数を作成
age = dataset[['Age','Pclass','Sex','Honorifics']]
age_dummies = pd.get_dummies(age)
age_dummies.head(3)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/2b07fc6f-0bea-7a94-5eb8-2617a1a9e0cf.png)

```py
# ageが分かるものと欠損しているものに分ける
known_age = age_dummies[age_dummies.Age.notnull()].as_matrix()
null_age = age_dummies[age_dummies.Age.isnull()].as_matrix()

# 特徴量と正解データに分ける
age_X = known_age[:, 1:]
age_y = known_age[:, 0]

# 年齢予測モデルを作成し、予測した値を代入
rf = RandomForestRegressor()
rf.fit(age_X, age_y)
pred_Age = rf.predict(null_age[:, 1:])
dataset.loc[(dataset.Age.isnull()),'Age'] = pred_Age
```
###3.2 'Embarked'(出港地)の欠損値の見直し
次に'Embarked'(出港地)の欠損値を埋めるために、欠損しているデータを確認します。

```py
# 'Embarked'(出港地)が欠損しているデータを表示
dataset[dataset['Embarked'].isnull()]
```
![2020-01-10 (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/88b8d4f7-20d3-40e5-d3f5-713c0aa3d9f5.png)
どちらも'Pclass'(チケットクラス)が１で、'Fare'(料金)が８０です。
'Pclass'(チケットクラス)が１の、'Embarked'(出港地)ごとの'Fare'(料金)中央値を比較すると、Cが最も近いです。２つの欠損値にはCを代入します。

```py
# 'Pclass'(チケットクラス)が１の、'Embarked'(出港地)ごとの'Fare'(料金)中央値を表示
C = dataset[(dataset['Embarked']=='C') & (dataset['Pclass'] == 1)]['Fare'].median()
print("Cの中央値", C)
S = dataset[(dataset['Embarked']=='S') & (dataset['Pclass'] == 1)]['Fare'].median()
print("Sの中央値", S)
Q = dataset[(dataset['Embarked']=='Q') & (dataset['Pclass'] == 1)]['Fare'].median()
print("Qの中央値", Q)

# 'Embarked'の欠損値にCを代入
dataset['Embarked'] = dataset['Embarked'].fillna('C')
```
Cの中央値 76.7292
Sの中央値 52.0
Qの中央値 90.0
###3.3 'Fare'(料金)の欠損値の見直し
データを確認すると、'Pclass'(チケットクラス)が3で、'Embarked'(出港地)が'S'だと分かります。
そのため、この欠損値には'Pclass'(チケットクラス)が3で、'Embarked'(出港地)が'S'の中央値を代入します。
これでAge,Embarked, Fareの欠損値を埋められたので確認します。

```py
# 'Fare'(料金)が欠損しているデータを表示
dataset[dataset['Fare'].isnull()]

# 'Pclass'(チケットクラス)が3で、'Embarked'(出港地)が'S'の中央値を代入
fare_median=dataset[(dataset['Embarked'] == "S") & (dataset['Pclass'] == 3)].Fare.median()
dataset['Fare']=dataset['Fare'].fillna(fare_median)

# 全体の欠損データの個数を確認
dataset_null = dataset.fillna(np.nan)
dataset_null.isnull().sum()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/aed3fb5a-c1ef-363b-c718-b3e8b505f3a8.png)
Age              0
Cabin            0
Embarked         0
Fare             0
Name             0
Parch            0
PassengerId      0
Pclass           0
Sex              0
SibSp            0
Survived       418
Ticket           0
Deck             0
TicketGroup      0
Honorifics       0
dtype: int64

欠損値がなくなっています。

###4.家族人数
[前々回](https://qiita.com/shiroino11111/items/bc3889fa38ff32d46c13)上手く扱えなかった同乗している兄弟/配偶者の数、親/子供の数を使えるデータに加工します
同乗する家族をまとめ、同乗家族の人数による生存率でグループ分けします。

```py
# 同乗している兄弟/配偶者の数による生存率比較
sns.barplot(x="SibSp", y="Survived", data=train, palette='Set3')

# 同乗している親/子供の数による生存率比較
sns.barplot(x="Parch", y="Survived", data=train, palette='Set3')

# 同乗している家族の数
dataset['FamilySize']=dataset['SibSp']+dataset['Parch']+1
sns.barplot(x="FamilySize", y="Survived", data=dataset, palette='Set3')

# 家族数による生存率でグループ分け
def Family_label(s):
    if (s >= 2) & (s <= 4):
        return 2
    elif ((s > 4) & (s <= 7)) | (s == 1):
        return 1
    elif (s > 7):
        return 0
dataset['FamilyLabel']=dataset['FamilySize'].apply(Family_label)
sns.barplot(x="FamilyLabel", y="Survived", data=dataset, palette='Set3')
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/54ee0915-6c66-1cca-94a6-eb953db12edd.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/3805c4b9-dbc4-5eb4-141b-364a9f1fa4f4.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/8d044ed9-df1b-c54d-1928-efc4ede01b71.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/91c0701f-4cc2-0334-cb83-369b4e589183.png)
きれいな差に分けることができました。

###5.苗字における生存率調整
上の'SibSp'と'Parch'では３等親以降の家族関係が不明なので、家族でなく苗字における生存率で調査します。
苗字における生存率で大きな差が確認できます。

```py
# 苗字による特徴を調べる
dataset['Surname'] = dataset['Name'].apply(lambda x:x.split(',')[0].strip()) # 苗字(名前の","の前の単語）を抽出
Surname_Count = dict(dataset['Surname'].value_counts()) # 苗字の個数を数える
dataset['Surname_Count'] = dataset['Surname'].apply(lambda x:Surname_Count[x]) # 苗字の個数を代入

# 苗字にダブりがある人を、女・子供のグループと、大人で男性のグループに分ける
Female_Child_Group=dataset.loc[(dataset['Surname_Count']>=2) & ((dataset['Age']<=12) | (dataset['Sex']=='female'))]
Male_Adult_Group=dataset.loc[(dataset['Surname_Count']>=2) & (dataset['Age']>12) & (dataset['Sex']=='male')]

# 女・子供グループにおける苗字ごとの生存率平均の個数を比較
Female_Child_mean = Female_Child_Group.groupby('Surname')['Survived'].mean() # 苗字ごとの生存率平均
Female_Child_mean_count = pd.DataFrame(Female_Child_mean.value_counts()) # 苗字ごとの生存率平均の数
Female_Child_mean_count.columns=['GroupCount']
Female_Child_mean_count
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/ecc8664b-99c8-2bc4-bc57-23bc3e319eaa.png)

```py
# 男（大人）グループにおける苗字ごとの生存率平均の個数を比較
Male_Adult_mean = Male_Adult_Group.groupby('Surname')['Survived'].mean() # 苗字ごとの生存率平均
Male_Adult_mean_count = pd.DataFrame(Male_Adult_mean.value_counts()) # 苗字ごとの生存率平均の数
Male_Adult_mean_count.columns=['GroupCount']
Male_Adult_mean_count
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/710c21d4-ad6a-ca51-49ce-9b9a78a5bf89.png)
どちらのグループも大抵１か０となり、グループで差が大きく存在することが分かります。
女子供（大人・男）を持つ家族と同じ苗字であれば、全員が生存（死亡）する、というルールでもあるのでしょうか？
これだけはっきりした特徴量は貴重です。
このルールと逆の結果となるものを外れ値として扱うことで、スコアの改善に貢献が期待できます。
何をするかというと、データの書き換えを行います。
女子供（大人・男）を持つ家族と同じ苗字であるのに、全員が死亡（生存）している苗字の人物は、全て逆側のルールに従わせたプロフィールのデータにします。

```py
# 各グループの例外を処理
# 各グループの例外となる苗字を抽出
# Dead_List：女・子供グループで全員死亡した苗字
# Survived_List：男（大人）グループで全員死亡した苗字
Dead_List = set(Female_Child_mean[Female_Child_mean.apply(lambda x:x==0)].index)
print("Dead_List", Dead_List, sep="\n")
Survived_List = set(Male_Adult_mean[Male_Adult_mean.apply(lambda x:x==1)].index)
print("Survived_List", Survived_List, sep="\n")
```
Dead_List
{'Danbom', 'Turpin', 'Zabour', 'Bourke', 'Olsson', 'Goodwin', 'Cacic', 'Robins', 'Canavan', 'Lobb', 'Palsson', 'Ilmakangas', 'Oreskovic', 'Lefebre', 'Sage', 'Johnston', 'Arnold-Franchi', 'Skoog', 'Attalah', 'Lahtinen', 'Jussila', 'Ford', 'Vander Planke', 'Rosblom', 'Boulos', 'Rice', 'Caram', 'Strom', 'Panula', 'Barbara', 'Van Impe'}
Survived_List
{'Chambers', 'Beane', 'Jonsson', 'Cardeza', 'Dick', 'Bradley', 'Duff Gordon', 'Greenfield', 'Daly', 'Nakid', 'Taylor', 'Frolicher-Stehli', 'Beckwith', 'Kimball', 'Jussila', 'Frauenthal', 'Harder', 'Bishop', 'Goldenberg', 'McCoy'}

```py
# テストデータを書き換える
# データをtrainとtestに分解 
train = dataset.loc[dataset['Survived'].notnull()]
test = dataset.loc[dataset['Survived'].isnull()]

# 女・子供グループで全員死亡した苗字の人→６０歳の男性、敬称はMrに。
# 男（大人）グループで全員生存した苗字の人→５才の女性、敬称はMissに。
test.loc[(test['Surname'].apply(lambda x:x in Dead_List)),'Sex'] = 'male'
test.loc[(test['Surname'].apply(lambda x:x in Dead_List)),'Age'] = 60
test.loc[(test['Surname'].apply(lambda x:x in Dead_List)),'Title'] = 'Mr'
test.loc[(test['Surname'].apply(lambda x:x in Survived_List)),'Sex'] = 'female'
test.loc[(test['Surname'].apply(lambda x:x in Survived_List)),'Age'] = 5
test.loc[(test['Surname'].apply(lambda x:x in Survived_List)),'Title'] = 'Miss'

# 再びデータを結合
dataset = pd.concat([train, test])
```

###6.再び予測を行う

```py
# 使用する変数を抽出
dataset6 = dataset[['Survived','Pclass','Sex','Age','Fare','Embarked','Honorifics','FamilyLabel','Deck','TicketGroup']]
# ダミー変数を作成
dataset_dummies = pd.get_dummies(dataset6)
dataset_dummies.head(3)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/479337/33a92082-aba0-6868-b738-e4d87d914042.png)

```py
# データをtrainとtestに分解 
#（ 'Survived'が存在するのがtrain, しないのがtest ）
train_set = dataset_dummies[dataset_dummies['Survived'].notnull()]
test_set = dataset_dummies[dataset_dummies['Survived'].isnull()]
del test_set["Survived"]

# trainデータを変数と正解に分離
X = train_set.as_matrix()[:, 1:] # Pclass以降の変数
y = train_set.as_matrix()[:, 0] # 正解データ

# 予測モデルの作成
pipe = Pipeline([('classify', RandomForestClassifier(random_state = 10, max_features = 'sqrt'))])

param_test = {'classify__n_estimators':list(range(20, 30, 1)), 
              'classify__max_depth':list(range(3, 10, 1))}
gsearch = GridSearchCV(estimator = pipe, param_grid = param_test, scoring='accuracy', cv=10)
gsearch.fit(X, y)
print(gsearch.best_params_, gsearch.best_score_)

# testデータの予測
predictions = gsearch.predict(test_set)

# Kaggle提出用csvファイルの作成
submission = pd.DataFrame({"PassengerId": PassengerId, "Survived": predictions.astype(np.int32)})
submission.to_csv("submission6.csv", index=False)
```
'classify__max_depth': 5, 'classify__n_estimators': 28}
0.8451178451178452
提出スコアは0.81818でした。

###7.特徴量を減らして予測
前回に比べて特徴量が２６個と大幅に増えたので、重要でない特徴量を除外します。

```py
pipe = Pipeline([('select',SelectKBest(k=20)),  # 予測に役立つ特徴量を２０個使ってモデルを作成
               ('classify', RandomForestClassifier(random_state = 10, max_features = 'sqrt'))])

param_test = {'classify__n_estimators':list(range(20, 30, 1)), 
              'classify__max_depth':list(range(3, 10, 1))}
gsearch = GridSearchCV(estimator = pipe, param_grid = param_test, scoring='accuracy', cv=10)
gsearch.fit(X, y)
print(gsearch.best_params_, gsearch.best_score_)
```
{'classify__max_depth': 6, 'classify__n_estimators': 26} 
0.8451178451178452

```py
select = SelectKBest(k = 20)
clf = RandomForestClassifier(random_state = 10, warm_start = True, 
                                  n_estimators = 26,
                                  max_depth = 6, 
                                  max_features = 'sqrt')
pipeline = make_pipeline(select, clf)
pipeline.fit(X, y)
```
先程のモデルとmax_depthとn_estimatorsが変化しました。
このmax_depthとn_estimators情報を用いて、再度特徴量を２０に絞って予測モデルを作成し、予測します。

```py
# 与えられたmax_depthとn_estimatorsを利用して、特徴量を２０に絞って再度予測モデルを作成し、予測
select = SelectKBest(k = 20)
clf = RandomForestClassifier(random_state = 10,
                             warm_start = True, 
                             n_estimators = 26,
                             max_depth = 6, 
                             max_features = 'sqrt')
pipeline = make_pipeline(select, clf)
pipeline.fit(X, y)

cv_score = model_selection.cross_val_score(pipeline, X, y, cv= 10)
print("CV Score : Mean - %.7g | Std - %.7g " % (np.mean(cv_score), np.std(cv_score)))

# testデータの予測
predictions = pipeline.predict(test_set)

# Kaggle提出用csvファイルの作成
submission = pd.DataFrame({"PassengerId": PassengerId, "Survived": predictions.astype(np.int32)})
submission.to_csv("submission7.csv", index=False)
```
CV Score : Mean - 0.8451402 | Std - 0.03276752 
これで提出スコアが0.83732となったはずです。
２０１９年時点での順位は２１７位です。これは上位1.5％に相当します。

###8.まとめ
欠損値を論理的に埋め、敬称など新たに特徴量を生成したり、テストデータの書き換えを行うことで、Kaggle Titanicの上位1.5％に相当する0.83732というスコアを出しました。
様々なデータ加工が出てきており、Titanicがデータ分析力のチュートリアルと扱われているのがよく分かります。

今回でTitanicは一旦終了です。
この記事を拝見した方のお役に立てれば幸いです。