Kaggle / Titanic - Machine Learning from Disaster / 独学7日間で上位5.6%を達成した話

## 目次
[1.なぜまとめようと思ったのか](#なぜまとめようと思ったのか)
[2.前提](#前提)
[3.実際に独学でやったこと](#実際に独学でやったこと)
[4.感想](#感想)

## なぜまとめようと思ったのか
[Kaggle](https://www.kaggle.com/)で入門用として有名なお題である"[Titanic - Machine Learning from Disaster](https://www.kaggle.com/competitions/titanic)"をやってみて7日間(1日最長2時間くらい)で上位5.6%(Bronze Medal相当)を達成しました。

Day1〜7までの具体的な取り組みについてのメモは[はてなブログ](https://blog.hatena.ne.jp/hideoyaji/hideoyaji.hatenablog.com/)に書きましたが、自分用でまとめておきたいのと、これからデータサイエンスでpythonを勉強する方にKaggleのコンペは敷居が高そうに見えるけど実際には楽しくできるよ、ってことを伝えたいのでまとめていきたいと思います。

## 前提
### 僕のpythonについての能力・学習機会
コーチやメンターがいる訳ではなく独学です。基本はProgateの有料課金コースを終えて、その後はWebサイトで収集できる情報をベースに環境構築をして、写経をしたり、Djangoフレームワークを使って趣味のロードバイクやトレランのレース情報をスクレイピングして一覧表示するサイトをつくったり([Cardio Freaks](http://cardiofreaks.pythonanywhere.com/))できるようになりましたが、駆け出しです。

また、これを書いている2023年3月2日時点で[株式会社JDSC](https://jdsc.ai/)の営業を僕の会社として受託して支援しており、その中でデータサイエンティストやエンジニアの仕事を横目で見ることができ、かつコミュニケーションする機会があります。

python x DSという意味では業務での実績がないことに変わりはありません。

pythonではないですがプログラミングでの仕事での実績としては、僕が新卒で入社して一人前の戦略コンサルタントに育ててくれた[株式会社野村総合研究所](https://www.nri.com/jp)の先輩が立ち上げた[D4DR株式会社](https://www.d4dr.jp/)のプロジェクトにジョインして、Google Spreadsheet上でGPT-3を使って企業情報を収集したくらいです。

### 今回のチャレンジにおける環境
+ 当初：Macbook Air(M1, 8GB/512GB)
+ 途中から：Macbooc Air(M2, 24GB/512GB)
+ Jupyter Lab

ちなみに、MBAのスペックが途中から上がったのはメモリ8GBのMBAはメモリプレッシャーがすぐに黄色になり、動作がもっさりしたためです。pythonが影響しているというよりは、Jupyter Labを動かしつつ、ブラウザのタブやウィンドウを多数開いたり、普段のコンサルで使うアプリ(Slack, Teams, Messenger等)を多数立ち上げていたことが原因だと思っています。24GBにして快適になりました(FYI)。

### 独学をやる前にKaggleでやったこと
1. KaggleのWebサイトをザッと眺めた
1. Kaggle x 初級者 等で検索して日本語サイトを読んだ
1. [Titanicテーマの動画](https://youtu.be/8yZMXCaFshs)を見た
1. 同、テーマサイトの中身を読んだ(英語は時にDeepLを活用)
1. [MANAV SHEGAL氏(Gold Medalist)のNotebook](https://www.kaggle.com/code/startupsci/titanic-data-science-solutions)を見ながら写経

特に5.は重要で、Notebookを公開されている方で詳しく解説をつけているものを見ることは有効でした。

Kaggleでのお題についての理解の仕方、分析の仕方、データの整形、特徴量生成、モデル選択、学習、評価、そしてKaggleへのsubmition(成果ファイルの提出)までの一通りを経験することができます。もちろん理解が追いついていないのですが、独学を終えた今ではそれでもとても有効だと思っています。

## 実際に独学でやったこと
### Day1, 2：データの理解、欠損(null)値の補完
#### データの理解
写経の時と同じようにまずはデータを可視化するところから始めました。
よく使ったのが下記です。
```python:jupyter lab
# pandasのDataFrameの概要を見る
train_df.info()
train_df.shape()

# pandasのDataFrameを眺める
train_df.head()
train_df.head(10)
train_df.tail(5)
train_df.sample(10)

# pandasのDataFrameの要約統計量(max/min, 四分位とか)を眺める
train_df.describe()
train_df.describe(include='all')
```
![はてなブログ](https://cdn-ak.f.st-hatena.com/images/fotolife/h/hideoyaji/20230221/20230221214725.png)
![はてなブログ](https://cdn-ak.f.st-hatena.com/images/fotolife/h/hideoyaji/20230221/20230221214826.png)
object型の変数があること、いくつかの変数にnull値があることが分かります。機械学習で使用するモデルの大半は(僕の今の理解では)object型、すなわち文字列に近いものは扱えませんし、null値があると処理ができません。ですので、変換、補完が必要です。

#### 変数の型変換
##### カテゴリー系のobject変数のint型への変換
次にobject型の変数について数値int型に変換しました。
```python:jupyter lab
# 学習、テスト用データ双方をfor文で扱うためにcombineを作成
combine = [[train_df], [test_df]]

# 'Sex'についてmale=0, female=1とする
for dataset in combine:
    dataset['Sex'] = dataset['Sex'].map( {'female':1, 'male':0} ).astype(int)

# 'Embarked'はtrain_dfで2つnull値があるので最頻値freq_portで埋める
freq_port = train_df.Embarked.dropna().mode()[0] # .mode()で最頻値をSeriesで取得
for dataset in combine:
    dataset['Embarked'] = dataset['Embarked'].fillna(freq_port)
# 'Embarked'を文字列型から0/1/2の数値カテゴリ型に変換
for dataset in combine:
    dataset['Embarked'] = dataset['Embarked'].map( {'S':0, 'C':1, 'Q':2} ).astype(int)
```

![はてなブログ](https://cdn-ak.f.st-hatena.com/images/fotolife/h/hideoyaji/20230221/20230221214906.png)

##### 文字列型objectである'Name'の変換
![スクリーンショット 2023-03-02 21.52.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/6d211410-e7e5-a8ea-660e-3e648642c4ee.png)
'Name'は文字列なのでそのままでは扱えません。dropさせても良いのですがここは写経の通りにやってみました。まずは正規化表現を用いて肩書部分('Mr'とか)を取り出します。
```python:Jupyter Lab
for dataset in combine:
    dataset['Title'] = dataset.Name.str.extract('([A-Za-z]+)\.', expand=False)

pd.crosstab(train_df_drpd['Title'], train_df_drpd['Sex'])
```
![スクリーンショット 2023-03-02 21.56.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/cc1bcba7-cc69-733f-6c92-6e0413b329f3.png)

次に肩書を一度くくります。ここは写経です。正直英語圏で生活していないと分からないと思います。
```python:Jupyter Lab
or dataset in combine:
    dataset['Title'] = dataset['Title'].replace(['Lady', 'Countess', 'Capt', 'Col', 'Don', 'Dr', 'Major', 'Rev', 'Sir', 'Jonkheer', 'Dona'], 'Rare')
    dataset['Title'] = dataset['Title'].replace('Mlle', 'Miss')
    dataset['Title'] = dataset['Title'].replace('Ms', 'Miss')
    dataset['Title'] = dataset['Title'].replace('Mme', 'Mrs')

train_df_drpd = train_df_drpd.drop(['Name'], axis=1)
test_df_drpd = test_df_drpd.drop(['Name'], axis=1)
combine = [train_df_drpd, test_df_drpd]
```
同時に'Name'はdorpしてしまいます。
![スクリーンショット 2023-03-02 21.58.03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/17f2a0df-1694-5c56-2262-f4aeeead1a91.png)
更に数値カテゴリー型にします。
```python:Jupyter Lab
title_mapping = {'Mr': 1, 'Miss': 2, 'Mrs': 3, 'Master': 4, 'Rare': 5}

for dataset in combine:
    dataset['TitleId'] = dataset['Title'].map(title_mapping)
    dataset['TitleId'] = dataset['TitleId'].fillna(0)

train_df_drpd = train_df_drpd.drop(['Title'], axis=1)
test_df_drpd = test_df_drpd.drop(['Title'], axis=1)
combine = [train_df_drpd, test_df_drpd]

train_df_drpd = train_df_drpd.drop(['Ticket'], axis=1)
test_df_drpd = test_df_drpd.drop(['Ticket'], axis=1)
combine = [train_df_drpd, test_df_drpd]
```
![スクリーンショット 2023-03-02 22.01.11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/a6ab52c2-e1be-f3ae-ec83-e4dd9cf13219.png)
上記では'Ticket'もdropしています。


#### null値の補完
##### 'Age'の補完
まずnull値の多い'Age'を補完します。傾向が見たいのでseabornのFacegrid関数を使って3(Pclass)x2(Survived)のグリッドを作成、グリッドの中にmatplotlibのhist関数で'Age'のヒストグラムを表示させてみます。
```python:Jupyter Lab
%matplotlib inline
grid = sns.FacetGrid(train_df_drpd, row='Pclass', col='Survived', hue='Survived') # col=''で指定した属性で分けてグラフを表示
grid.map(plt.hist, 'Age', bins=20) # Ageを軸として20分割(棒の数)でヒストグラム表示(この場合80歳/20なので4歳刻み)
```
ちなみにseabornとかmatplotlibを使う際に`%matplotlib inline`をいちいち書かないと2度目の実行からはグラフが出力されないということがありました。ですので、おまじないのように書いてます。
![スクリーンショット 2023-03-02 21.41.54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/0993f355-d7e5-2dab-67a1-5ec99c2f7a4a.png)

同様に'Survived'と'Sex'を入れ替えて描画。
```python:Jupyter Lab
%matplotlib inline
grid = sns.FacetGrid(train_df_drpd, row='Pclass', col='Sex', hue='Sex') # col=''で指定した属性で分けてグラフを表示
grid.map(plt.hist, 'Age', bins=20) # Ageを軸として20分割(棒の数)でヒストグラム表示(この場合80歳/20なので4歳刻み)
```
![スクリーンショット 2023-03-02 21.43.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/4d6f5b35-a75f-10f3-32da-87d7d2c57e11.png)

'Age'は'Pclass'とも相関していますし、'Sex'によって中央値・平均値が微妙にずれています。ですので'Pclass'と'Sex'のマトリクスを作って中央値Median()で補完することとしました。

```python:Jupyter Lab
# まず'Sex', 'Pclass'で行列を準備
guess_ages = np.zeros((2, 3))

#  行列'guess_ages'に'Sex', 'Pclass'における'Age'の中央値を代入
for dataset in combine:
    for i in range(0, 2):
        for j in range(0, 3):
            # 'Sex'がiで'Pclass'がj+1のレコードをDataFrameとして抽出
            guess_df = dataset[(dataset['Sex'] == i) &
                               (dataset['Pclass'] == j+1)]['Age'].dropna()
            # 中央値を算出
            age_guess = guess_df.median()
            # 行列(i, j)に中央値を代入
            guess_ages[i, j] = int(age_guess/0.5 + 0.5) * 0.5  # ここは分からない…職人技？

    for i in range(0, 2):
        for j in range(0, 3):
            dataset.loc[(dataset.Age.isnull()) & (dataset.Sex == i) &
                        (dataset.Pclass == j+1), 'Age'] = guess_ages[i, j]

    dataset['Age'] = dataset['Age'].astype(int)
```
![スクリーンショット 2023-03-02 21.46.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/1d63a1f0-52fd-27e5-7a3a-527440189fb5.png)

##### 'Fare'の補完
テスト用データであるpandas DataFrameのtest_df_drpdの'Fare'に1つだけnull値があるのでそれを中央値で補完しました。
```python:Jupyter Lab
test_df_drpd['Fare'].fillna(test_df_drpd['Fare'].dropna().median(), inplace=True)
```
![スクリーンショット 2023-03-02 22.05.12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/619e3640-edbd-1e96-42c5-b41be0812fc6.png)
これにて全てがint64 / float64という数値型の変数、特徴量となりました。データのクレンジングが終わったということですね。ここまででDay1,2でした。

### 目的変数と説明変数、説明変数同士の相関関係の把握
#### 変数同士の相関係数をマトリクスで可視化
```python:Jupyter Lab
train_df_drpd_corr = train_df_drpd.corr()
print(train_df_drpd_corr)
```
![スクリーンショット 2023-03-02 22.07.47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/fb1dd4ae-54eb-1dc0-effe-f9c2b8f51b1d.png)
相関係数の表はよく使いました。これを見ながら目的変数である'Survived'との相関度が高い変数をピックアップしました。

+ 'Survived'と相関が高いのは'Sex', 'TitleId', 'Pclass', 'Fare'
+ 'Sex'と'TitleId'は相関度が高い
+ 'Pclass'と’Fare', 'Age'も相関度が高い
+ 'SibSp'と’Parch'は相関度が高い
+ 'Embarked'は特に何も相関していない

あたりが分かったので、説明変数同士の相関(多重共線性)を排除することを意識して、単純に’Sex'と'Pclass'のみでモデリングしてみることにしました。

#### モデリングと評価
##### Day3：1st Submission
機械学習モデルについては写経の時にRandom Forestが一番効いているのが分かったので、今回の独学ではRandom Forest一択としています。
```python:Jupyter Lab
from sklearn.ensemble import RandomForestClassifier

t = train_df_drpd['Survived'].values
x = train_df_drpd[['Sex', 'Pclass']].values
x_test = test_df_drpd[['Sex', 'Pclass']].values

rf = RandomForestClassifier()
rf.fit(x, t)
pred_rf = rf.predict(x_test)
rf.score(x, t)

submission = pd.DataFrame({
    'PassengerId':test_df_drpd['PassengerId'],
    'Survived':pred_rf
})
```
+ モデルスコア：0.7867564534231201
+ テストスコア：0.77511

参加者の中で下位30%くらいだったかと思います。ここまでがDay3です。

その後、上記の仮説で出した特徴量を色々と組み合わせてテストしますが、この0.77511をなかなか超えられず一度目の壁にぶつかります。

##### Day4,5：'IsInfant'をモデルに組み込み上位7.2％を達成
下記の特徴量を作ってみました。
+ FamilySize：家族構成人数で'SibSp'と'Parch'を加算したもの
+ IsSmallFamily:'FamilySize'<=3の人
+ IsSingle：'FamilySize'=0の人
+ IsInfant：'Age'が7歳未満の人
```python:Jupyter Lab
# FamilySize
for dataset in combine:
    dataset['FamilySize'] = dataset['SibSp'] + dataset['Parch']
train_df_drpd[['FamilySize', 'Survived']].groupby(['FamilySize'], as_index=False).mean().sort_values(by='FamilySize', ascending=True)

# IsSmallFamily
for dataset in combine:
    dataset.loc[(dataset['FamilySize'] >= 1) & (dataset['FamilySize'] <= 3), 'IsSmallFamily'] = 1
    dataset.loc[(dataset['FamilySize'] < 1) | (dataset['FamilySize'] >= 4), 'IsSmallFamily'] = 0

train_df_drpd['IsSmallFamily'] = train_df_drpd['IsSmallFamily'].astype(int)
test_df_drpd['IsSmallFamily'] = test_df_drpd['IsSmallFamily'].astype(int)

combine = [train_df_drpd, test_df_drpd]

# IsSingle
for dataset in combine:
    dataset.loc[dataset['FamilySize'] == 0, 'IsSingle'] = 1
    dataset.loc[dataset['FamilySize'] != 0, 'IsSingle'] = 0

train_df_drpd['IsSingle'] = train_df_drpd['IsSingle'].astype(int)
test_df_drpd['IsSingle'] = test_df_drpd['IsSingle'].astype(int)

combine = [train_df_drpd, test_df_drpd]

# IsInfant
for dataset in combine:
    dataset.loc[dataset['Age'] < 7, 'IsInfant'] = 1
    dataset.loc[dataset['Age'] >= 7, 'IsInfant'] = 0
    
train_df_drpd['IsInfant'] = train_df_drpd['IsInfant'].astype(int)
test_df_drpd['IsInfant'] = test_df_drpd['IsInfant'].astype(int)
combine = [train_df_drpd, test_df_drpd]
```
![スクリーンショット 2023-03-02 22.22.40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/3276e7dd-d718-bc1e-159a-855c7a62eb09.png)

それで色々と特徴量を組み合わせてみました。
```python:Jupyter Lab
t = train_df_drpd['Survived'].values
x = train_df_drpd[['Sex', 'Pclass', 'IsInfant']].values
x_test = test_df_drpd[['Sex', 'Pclass', 'IsInfant']].values

rf.fit(x, t)
pred_rf = rf.predict(x_test)
rf.score(x, t)
```
+ モデルスコア：0.8058361391694725
+ テストスコア：0.77751

ベストスコアを更新。順位も3,745/14,205と上位26.4％と大躍進しました。
ここまでがDay4です。

Day5ではDay4で作成したモデルに更に’TitleId’を追加して学習させました。
```python:Jupyter Lab
t = train_df_drpd['Survived'].values
x = train_df_drpd[['Sex', 'Pclass', 'IsInfant', 'TitleId']].values
x_test = test_df_drpd[['Sex', 'Pclass', 'IsInfant', 'TitleId']].values

rf.fit(x, t)
pred_rf = rf.predict(x_test)
rf.score(x, t)
```
+ モデルスコア：0.8080808080808081
+ テストスコア：0.78947
![20230226174037.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/65d0b9c7-3153-99d7-d01f-c2df18a2915a.png)

結果'Sex', 'Pclass', 'IsInfant', 'TitleId'の組み合わせとなりましたが、これでベストスコアを更新し1,038 / 14,269で上位7.2%を記録することができました。

##### Day6,7：'Fare'に関する特徴量を組み込み、上位5.6%を達成
その後伸び悩み、一度基本に立ち返るということで学習用データを眺めてみました。
![20230301201821.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/633510d8-ced5-f38f-92ce-98bed86a67e9.png)
Day5でのBest Scoreを叩き出したモデルの特徴量は'Sex', 'Pclass', 'TitleId', 'IsInfant'の4つ。性別、階級、家族的な要素が含まれたものになっていることに気づきました。

そして’Fare'については特徴量化していない。’Fare'も’Survived'と相関係数はまあまあ高いもののこれもモデルに入れると過学習が起こり、テストスコアが下がることも分かっていました。

更にもう一度テーブルを眺めてみて、'Fare'を'FamiliSize(+1)'で割って'UnitFare'としてみたらどうだろうかと思いつきました。
```python:Jupyter Lab
train_df_drpd['UnitFareBand'] = pd.cut(train_df_drpd['UnitFare'], 10)
train_df_drpd[['UnitFareBand', 'Survived']].groupby(['UnitFareBand'], as_index=False).mean().sort_values(by='UnitFareBand', ascending=True)
```
![20230301202114.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/f3384ad4-d9b3-f578-00bb-411d12109793.png)
pandasの.cut()の変数を色々と変えて'Survived'の平均値を見てみたところ、'UnitFare'<50あたりで生存率が35%前後と極めて低くなることが分かりました。

要するに1人あたりのチケット金額が50を超えている人は生き残る確率が高いということですね。ということで’IsLowUnitPrice'という特徴量を生成しました。
```python:Jupyter Lab
for dataset in combine:
    dataset.loc[dataset['UnitFare'] < 50, 'IsLowUnitFare'] = 0
    dataset.loc[dataset['UnitFare'] >= 50, 'IsLowUnitFare'] = 1

train_df_drpd['IsLowUnitFare'] = train_df_drpd['IsLowUnitFare'].astype(int)
test_df_drpd['IsLowUnitFare'] = test_df_drpd['IsLowUnitFare'].astype(int)

combine = [train_df_drpd, test_df_drpd]
```
+ モデルスコア：0.8114478114478114
+ テストスコア：0.79186
![20230301202554.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224055/8e186e5d-ffb0-8c58-6330-da551566e35a.png)
これで799位 / 14,182人で上位5.6％を達成しました。独学を始めてDay7のことでした。

# 感想
+ Day1,2：データクレンジング
+ Day3：モデリング開始と最初のサブミッションで下位30%からスタート
+ Day4：特徴量生成とモデリングにより上位26.4%
+ Day5：特徴量の組み合わせを変更して上位7.2%
+ Day6：ベストモデルに変数、特徴量を追加して低迷
+ Day7：触っていなかった変数から新特徴量を生成して上位5.6%を達成

こんな感じになりました。

ここからの学びとしては、
+ 変数、特徴量で目的変数への相関係数0.3以上をモデルに組み込むことが基本
+ 上記に達していないものでも独立性の高い変数、特徴量は精度向上に寄与する
+ 上記でスタックした場合には使っていない独立性の高い変数から特徴量を生成しモデルに組み込むことでブレイクスルーできる可能性がある

というあたりです。再現性のないやり方をしている初心者なのでもっと効率よく特徴量エンジニアリング、モデリングをしたいです。そのためにmatplotlibやseaborn等を活用してデータのビジュアライズを瞬時にできるように熟練する、というのが目下の課題です。

長文、乱文にお付き合いいただき、ありがとうございました。