最新アンサンブル学習SklearnStackingの性能調査(LBGM, RGF, ET, RF, LR, KNNモデルをHeamyとSklearnで比較する)



# TL;DR
sklearnのスタッキング、使ってみたらすごい優秀な子だったので標準になるかもしれない

| library                    | layer1                        | layer2  | public score | speed |
| -------------------------- | ----------------------------- | ------- | ------------ | ----- |
| lgbm model                 | lgb                           | -       | 0.74162      |    -   |
| heamy single stacking      | lbg, rgf, et, rf, lr, knn     | lr      | 0.76076      |    64s   |
| heamy multiple stacking    | lgb, rgf, et, rf, lr, knn     | lgb, lr | 0.74162      |    76s (64 + 12)   |
| heamy single stacking   nn | lgb, rgf, et, rf, lr, knn, nn | lgb, nn | 0.76076      |    119s   |
| heamy multiple stacking nn | lgb, rgf, et, rf, lr, knn, nn | lgb, nn | 0.75598      |    139s (119 + 20)   |
| sklearn single stacking    | lgb, rgf, et, rf, lr, knn     | lr      | 0.77511      |    30s   |
| sklearn multiple  stacking | lgb, rgf, et, rf, lr, knn     | lgb, lr | 0.76555      |    28s   |

# Abstract

[white, inc](http://bit.ly/33fuxi1) の ソフトウェアエンジニア [r2en](http://bit.ly/2qGILex) です。 
自社では新規事業を中心としたコンサルタント業務を行なっており、
普段エンジニアは、[新規事業を開発する無料のクラウド型ツール](http://bit.ly/33iJeRe) を開発したり、
新規事業のコンサルティングからPoC開発まで携わります

[![image](https://user-images.githubusercontent.com/17031124/70364775-4b5f7500-18d1-11ea-91bc-857813d89762.png)](http://bit.ly/33iJeRe)

今回は、機械学習の技術調査を行なったので記事で共有させていただきます
以下から文章が長くなりますので、口語で記述させていただきます

----------------------

scikit-learn 0.22で新しく、アンサンブル学習のStackingを分類と回帰それぞれに使用できるようになったため、自分が使っているHeamyと使用感を比較する

KaggleのTitanicデータセットを使い、性能や精度、速度を検証する

アンサンブルに使用する機械学習モデルは、lightgbm, regularized greedy forest, extremely randomized trees, random forest, logistic regression, K Nearest Neighbor, 3layer nural network になる

# Introduction

そもそもアンサンブル学習とは、複数の機械学習モデルを組み合わせてモデルを作り、予測することを指す。

単一のモデルよりもアンサンブルしたモデルは高精度になることが多く、分析コンペでは多く取り入れられている手法である

アンサンブル学習にも、平均や加重平均、StackingやBlendingなど様々な手法がある

現在では、自作stacking、[pystacknet](https://github.com/h2oai/pystacknet)、[heamy](https://github.com/rushter/heamy)などが主流で使われているように思う

ここでは実装が量を多く含めてしまうため、詳細なことは割愛させていただき、よりわかりやすい説明が載っている以下の文献を参考にしていただきたい

[Kaggle Ensembling Guide](https://mlwave.com/kaggle-ensembling-guide/)
[Kaggleでかつデータ分析の技術](https://www.amazon.co.jp/gp/product/4297108437/ref=ppx_yo_dt_b_asin_title_o04_s00?ie=UTF8&psc=1)
[アンサンブル手法のStackingを実装と図で理解する](https://blog.ikedaosushi.com/entry/2018/10/21/204842)

# Environment

#### 検証マシン

```
OS: macOS HighSierra 10.13.6(Retina, Early 2015)
CPU: 3.1GHz Intel Core i7
MEM: 16GB 1867MHz DDR3
GPU: Intel Iris Graphics 6100 1536MB
```

#### 検証環境

テスト等は行なっていない
dockerfileを使用している
jupyternotebookを使用している
もし追試する場合は、以下リンクから僕のリポジトリをクローンして欲しい
実行の手順書も.mdに記載している

https://github.com/r2en/research_sklearn_stacking

#### インストール

自身でインストールされる場合の方法も載せておく

```shell
$ pip install -U scikit-learn==0.22.0
```

```python
import sklearn
sklearn.__version__
'0.22.0'
```

#### 速度計測

```python
import time
from contextlib import contextmanager

@contextmanager
def timer(name):
    t0 = time.time()
    yield
    print(f'[{name}] done in {time.time() - t0:.0f} s')
```

#### 検証データ


[Titanic: Machine Learning from Disaster](https://www.kaggle.com/c/titanic/data)

通称タイタニックコンペ

わかりやすい、読まれやすい、スタック部分だけ集中して見てもらいやすい、追試してもらいやすい、
時系列データではない、学習とテストの分布が似ているため最低限のスタッキングできる要件は満たしている
という部分で採用した

ただ、Stackingするにはデータ数少なすぎる、評価指標が正答率でブレやすいなどのデメリットもあるため、
あくまでも実装方法の参考や、性能指標もなんとなくn=1のデータセットにはこういう精度なんだという気持ちで見て欲しい

このコンペは、タイタニック号に乗船した各乗客のデータを元に、タイタニック号が氷山に衝突し沈没した際生存したかどうか(Survived)を予測する

データの説明変数と目的変数は以下

| 変数名 | 特徴 |
| -------- | -------- |
| PassengerId　|乗客識別ユニークID|
|Survived|生死|
|Pclass|チケットクラス|
|Name|乗客の名前|
|Sex|性別|
|Age|年齢|
|SibSp|タイタニックに同乗している兄弟/配偶者の数|
|Parch|タイタニックに同乗している親/子供の数|
|Ticket|チケット番号|
|Cabin|客室番号|
|Embarked|出港地（タイタニックへ乗った港）|

#### 前処理

アンサンブル同士の性能比較であり、タイタニックコンペでの高い順位を求めているわけではない為、最低限どの機械学習モデルでも学習予測できる最低限の前処理しか行わない

正規化をしているのはNuralNetworkモデルが学習予測できるようにする為

「他のモデルに影響を及ばさないのか」については、[決定木は本当に変換に依存しないのか？](https://amalog.hateblo.jp/entry/decision-tree-scaling)の記事の説明がわかりやすい

もし、影響を多少及ぼしても今回は無視するものとする


```python
import re
import numpy
import pandas
from sklearn.preprocessing import StandardScaler
```

```python
def first_dataset():
    
    train = pandas.read_csv('train.csv')
    test = pandas.read_csv('test.csv')
    
    datasets = [train, test]
    
    def get_title(name):
        if re.search(' ([A-Za-z]+)\.', name):
            return re.search(' ([A-Za-z]+)\.', name).group(1)
        return ""

    
    for dataset in datasets:

        dataset['Cabin'] = dataset['Cabin'].apply(lambda x: 1 if type(x) == str else 0)
        
        dataset['Age'] = dataset['Age'].fillna(-1).astype(int)
        
        dataset['Fare'] = dataset['Fare'].fillna(-1).astype(int)

        dataset['Sex'] = dataset['Sex'].map( {'female': 0, 'male': 1} ).astype(int)
        
        dataset['Title'] = dataset['Name'].apply(get_title)
        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')
        title_mapping = {"Mr": 1, "Miss": 2, "Mrs": 3, "Master": 4, "Rare": 5}
        dataset['Title'] = dataset['Title'].map(title_mapping)
        dataset['Title'] = dataset['Title'].fillna(-1)

        dataset['Embarked'] = dataset['Embarked'].fillna('S')
        dataset['Embarked'] = dataset['Embarked'].map( {'S': 0, 'C': 1, 'Q': 2} ).astype(int)
    
        dataset.drop(['PassengerId', 'Ticket', 'Name'], axis=1, inplace=True)
        
    X_train = train.drop(['Survived'], axis=1)
    y_train = train['Survived']
    X_test = test
    
    std = StandardScaler()
    std.fit(X_train)
    X_train = std.transform(X_train).astype(numpy.float32)
    X_test = std.transform(X_test).astype(numpy.float32)

    return {'X_train': X_train, 'X_test': X_test, 'y_train': y_train}
    
df = first_dataset()
```

# Method

## Single LightGBM

スタッキングすること自体が、有用かどうかを検証するため
単体モデルの予測性能を測る


```python
import lightgbm as lgbm
from sklearn.model_selection import train_test_split

X_train, X_valid, y_train, y_valid = train_test_split(df['X_train'], df['y_train'], test_size=0.2, random_state=0)

train_dataset = lgbm.Dataset(data=X_train, label=y_train, free_raw_data=False)
test_dataset = lgbm.Dataset(data=X_valid, label=y_valid, free_raw_data=False)
final_train_dataset = lgbm.Dataset(data=df['X_train'], label=df['y_train'], free_raw_data=False)

lgbm_params = {
    'boosting': 'dart', 
    'application': 'binary',
    'learning_rate': 0.05,
    'min_data_in_leaf': 20,
    'feature_fraction': 0.7,
    'num_leaves': 41,
    'metric': 'binary_logloss',
    'drop_rate': 0.15
}

evaluation_results = {}
clf = lgbm.train(train_set=train_dataset,
                 params=lgbm_params,
                 valid_sets=[train_dataset, test_dataset], 
                 valid_names=['Train', 'Test'],
                 evals_result=evaluation_results,
                 num_boost_round=500,
                 early_stopping_rounds=100,
                 verbose_eval=20
                )
                
clf_final = lgbm.train(train_set=final_train_dataset,
                      params=lgbm_params,
                      num_boost_round=500,
                      verbose_eval=0
                      )

y_pred = numpy.round(clf_final.predict(df['X_test'])).astype(int)

passengerId = pandas.read_csv('test.csv')['PassengerId']
dataframe = pandas.DataFrame({'PassengerId': passengerId, 'Survived': y_pred})

dataframe.to_csv('submission_single_lgbm_model.csv', index=False)
```

## Heamy single Stacking


```python
from heamy.dataset import Dataset
from heamy.estimator import Regressor, Classifier
from heamy.pipeline import ModelsPipeline

from rgf.sklearn import RGFClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier
```


```python
ds = Dataset(preprocessor=first_dataset, use_cache=False)
```

```python
et_params = {'n_estimators': 100, 'max_features': 0.5, 'max_depth': 18, 'min_samples_leaf': 4, 'n_jobs': -1}
rf_params = {'n_estimators': 125, 'max_features': 0.2, 'max_depth': 25, 'min_samples_leaf': 4, 'n_jobs': -1}
rgf_params = {'algorithm': 'RGF_Sib', 'loss': 'Log'}
```

```python
from keras.layers import Dense
from keras.models import Sequential

def NuralNetClassifier(X_train, y_train, X_test, y_test=None):
    input_dim = X_train.shape[1]
    
    model = Sequential()
    model.add(Dense(12, input_dim=input_dim, activation='relu'))
    model.add(Dense(6, activation='relu'))
    model.add(Dense(1, activation='sigmoid'))
    model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])
    model.fit(X_train, y_train, epochs=30, batch_size=10, verbose=0)
    y_pred = numpy.ravel(model.predict(X_test))
    
    return y_pred
```


```python    
def LightGBMClassifier(X_train, y_train, X_test, y_test=None):
    lgbm_params = {
        'boosting': 'dart', 
        'application': 'binary',
        'learning_rate': 0.05,
        'min_data_in_leaf': 20,
        'feature_fraction': 0.7,
        'num_leaves': 41,
        'metric': 'binary_logloss',
        'drop_rate': 0.15
    }
    
    X_train, X_valid, y_train, y_valid = train_test_split(X_train, y_train, test_size=0.2, random_state=0)
    train_dataset = lgbm.Dataset(data=X_train, label=y_train, free_raw_data=False)
    test_dataset = lgbm.Dataset(data=X_valid, label=y_valid, free_raw_data=False)
    
    final_train_dataset = lgbm.Dataset(data=X_train, label=y_train, free_raw_data=False)
    
    evaluation_results = {}
    
    clf = lgbm.train(train_set=train_dataset,
                     params=lgbm_params,
                     valid_sets=[train_dataset, test_dataset], 
                     valid_names=['Train', 'Test'],
                     evals_result=evaluation_results,
                     num_boost_round=500,
                     early_stopping_rounds=100,
                     verbose_eval=0
                    )
    
    clf_final = lgbm.train(train_set=final_train_dataset,
                          params=lgbm_params,
                          num_boost_round=500,
                          verbose_eval=0
                          )

    y_pred = clf_final.predict(X_test)

    
    return y_pred
```

```python
pipeline = ModelsPipeline(
    Classifier(estimator=LightGBMClassifier, dataset=ds, use_cache=False),
    Classifier(estimator=NuralNetClassifier, dataset=ds, use_cache=False),
    Classifier(estimator=RGFClassifier, dataset=ds, use_cache=False, parameters=rgf_params),
    Classifier(estimator=ExtraTreesClassifier, dataset=ds, use_cache=False, parameters=et_params),
    Classifier(estimator=RandomForestClassifier, dataset=ds, use_cache=False, parameters=rf_params),
    Classifier(estimator=LogisticRegression, dataset=ds, use_cache=False),
    Classifier(estimator=KNeighborsClassifier, dataset=ds, use_cache=False)
)
```

```python
stack_ds = pipeline.stack(k=10, seed=0, add_diff=False, full_test=True)
stacker = Classifier(dataset=stack_ds, estimator=LogisticRegression, use_cache=False)
y_pred = stacker.predict()
```

```python
dataframe = pandas.DataFrame({'PassengerId': pandas.read_csv('test.csv')['PassengerId'], 'Survived': numpy.round(y_pred).astype(int)})
dataframe.to_csv('submission_heamy_single_stacking_model.csv', index=False)
```


## Heamy multiple Stacking


```python
from sklearn.metrics import log_loss
```

```python
pipeline2 = ModelsPipeline(
    Classifier(estimator=LightGBMClassifier, dataset=stack_ds, use_cache=False),
    Classifier(estimator=NuralNetClassifier, dataset=stack_ds, use_cache=False)
)
```

```python
weights = pipeline2.find_weights(log_loss)
predictions = pipeline2.weight(weights).execute()
```

```
Best Score (log_loss): 0.3849248890947457
Best Weights: [0.50000511 0.49999489]
```

```python
dataframe = pandas.DataFrame({'PassengerId': pandas.read_csv('test.csv')['PassengerId'], 'Survived': numpy.round(predictions).astype(int)})
dataframe.to_csv('submission_heamy_multiple_stacking_model.csv', index=False)
```


## sklearn single Stacking


```python
from sklearn.ensemble import StackingClassifier
from keras.wrappers.scikit_learn import KerasClassifier
```

```python
lgbm_params = {
        'boosting': 'dart', 
        'application': 'binary',
        'learning_rate': 0.05,
        'min_data_in_leaf': 20,
        'feature_fraction': 0.7,
        'num_leaves': 41,
        'metric': 'binary_logloss',
        'drop_rate': 0.15
}
keras_params = {'epochs': 10, 'batch_size': 10}
```

```python
def build_fn():
    clf = Sequential()
    clf.add(Dense(12, input_dim=9, activation='relu'))
    clf.add(Dense(6, activation='relu'))
    clf.add(Dense(1, activation='sigmoid'))
    clf.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])
    return clf
```


今回本当に申し訳ないが、KerasClassifierというsklearn準拠のラッパーおよび、自作でsklearn準拠のestimatorの中にkeras nnを差し込んだが、エラーを起こして動かなかった

僕の実装力不足ですね

なので、sklearnではnnを使用したアンサンブル学習ができてない

しかし、heamyにも記載はしないが同様にNNを抜いたアンサンブルの環境での実験は行なっている為安心していただきたい

下記が当該エラーになる

```
python3.7/dist-packages/sklearn/ensemble/_base.py in _validate_estimators(self)
    249                 raise ValueError(
    250                     "The estimator {} should be a {}.".format(
--> 251                         est.__class__.__name__, is_estimator_type.__name__[3:]
    252                     )
    253                 )

ValueError: The estimator KerasClassifier should be a classifier.
```

```python
estimators = [
    ('lgb', lgbm.LGBMClassifier(**lgbm_params)),
    #('nn', KerasClassifier(build_fn=build_fn, **keras_params)),
    ('rgf', RGFClassifier(**rgf_params)),
    ('et', ExtraTreesClassifier(**et_params)),
    ('rf', RandomForestClassifier(**rf_params)),
    ('lr', LogisticRegression()),
    ('knn', KNeighborsClassifier())
]
```

```python
clf = StackingClassifier(estimators=estimators, final_estimator=LogisticRegression())
```

```python
clf.fit(df['X_train'], df['y_train'])
predictions = clf.predict(df['X_test'])
dataframe = pandas.DataFrame({'PassengerId': pandas.read_csv('test.csv')['PassengerId'], 'Survived': numpy.round(predictions).astype(int)})
dataframe.to_csv('submission_sklearn_single_stacking_model.csv', index=False)
```

## sklearn multiple Stacking


```python
final_estimator = StackingClassifier(
    estimators= [
        ('lgb', lgbm.LGBMClassifier(**lgbm_params)),
        ('lr', LogisticRegression())
    ],
    final_estimator=LogisticRegression()
)
```

```python
clf = StackingClassifier(
    estimators= [
        ('lgb', lgbm.LGBMClassifier(**lgbm_params)),
        #('nn', KerasClassifier(build_fn=build_fn, **keras_params)),
        ('rgf', RGFClassifier(**rgf_params)),
        ('et', ExtraTreesClassifier(**et_params)),
        ('rf', RandomForestClassifier(**rf_params)),
        ('lr', LogisticRegression()),
        ('knn', KNeighborsClassifier())
    ],
    final_estimator=final_estimator
)
```

```python
clf.fit(df['X_train'], df['y_train'])
predictions = clf.predict(df['X_test'])
dataframe = pandas.DataFrame({'PassengerId': pandas.read_csv('test.csv')['PassengerId'], 'Survived': numpy.round(predictions).astype(int)})
dataframe.to_csv('submission_sklearn_multiple_stacking_model.csv', index=False)
```


# Result

sklearn stacking で NuralNetworkが使えないことで、比較対象を多くせざるおえなくなったが、以下の通りの結果になった

| library                    | layer1                        | layer2  | public score | speed |
| -------------------------- | ----------------------------- | ------- | ------------ | ----- |
| lgbm model                 | lgb                           | -       | 0.74162      |    -   |
| heamy single stacking      | lbg, rgf, et, rf, lr, knn     | lr      | 0.76076      |    64s   |
| heamy multiple stacking    | lgb, rgf, et, rf, lr, knn     | lgb, lr | 0.74162      |    76s (64 + 12)   |
| heamy single stacking   nn | lgb, rgf, et, rf, lr, knn, nn | lgb, nn | 0.76076      |    119s   |
| heamy multiple stacking nn | lgb, rgf, et, rf, lr, knn, nn | lgb, nn | 0.75598      |    139s (119 + 20)   |
| sklearn single stacking    | lgb, rgf, et, rf, lr, knn     | lr      | 0.77511      |    30s   |
| sklearn multiple  stacking | lgb, rgf, et, rf, lr, knn     | lgb, lr | 0.76555      |    28s   |

# Discussion

- 一番高精度だったのは、sklearnという結果になった、次点で、NNありのHeamyという結果になった
- multiple stackingよりもsingle stackingの方が性能が良いという結果がでた <- データ数が少ない為、多段スタッキングよりも一段スタッキングの方が各ライブラリとも良い値がでたことが理由と思われる
- heamyよりもskleanの方のスタッキングの方が良い精度でる
- お手軽さ・わかりやすさ・メンテナンス性は断然sklearn
- 中身を複雑に変えたい場合は、heamyの方が良い気がする
- ただ、sklearnのスタッキングにはpipelineも取り込めるので少しだけ複雑なことはできる
- 速度が圧倒的にsklearnスタッキングの方が上となった
- sklearnには各モデルの重みを設定できるところが見当たらなかった

NuralNetworkが扱えるようになれば、sklearnのスタッキングは大変良いのでは？

# Reference

- [Release Highlights for scikit-learn 0.22](https://scikit-learn.org/stable/auto_examples/release_highlights/plot_release_highlights_0_22_0.html)
- [1.11.8. Stacked generalization¶](https://scikit-learn.org/stable/modules/ensemble.html#stacking)
- [heamy documentation](https://heamy.readthedocs.io/en/latest/usage.html)
- [github: heamy](https://github.com/rushter/heamy)
- [Kaggle: Introduction to Ensembling/Stacking in Python](https://www.kaggle.com/arthurtok/introduction-to-ensembling-stacking-in-python)
- [Kaggle: Stacking using heamy](https://www.kaggle.com/rushter/stacking-using-heamy)
- [Forest Cover -Ensembling and Stacking with heamy](https://www.kaggle.com/justfor/ensembling-and-stacking-with-heamy#Forest-Cover--Ensembling-and-Stacking-with-heamy)
- [アンサンブル手法のStackingを実装と図で理解する](https://blog.ikedaosushi.com/entry/2018/10/21/204842)
- [KAGGLE ENSEMBLING GUIDE](https://mlwave.com/kaggle-ensembling-guide/)
- [Kaggleでよく使われるStacking/Blendingをheamy、Stacknetをpystacknetで高速に実装する](https://blog.ikedaosushi.com/entry/2018/12/15/212508)
- [Qiita: Pythonでアンサンブル(スタッキング)学習 & 機械学習チュートリアル in Kaggle](https://qiita.com/hokuto_HIRANO/items/2c35a81fbc95f0e4b7c1#%E6%9C%80%E5%BE%8C%E3%81%AB)
- [Amazon: Kaggleでかつデータ分析の技術](https://www.amazon.co.jp/gp/product/4297108437/ref=ppx_yo_dt_b_asin_title_o04_s00?ie=UTF8&psc=1)
- [github: pystacknet](https://github.com/h2oai/pystacknet)
- [決定木は本当に変換に依存しないのか？](https://amalog.hateblo.jp/entry/decision-tree-scaling)
- [Keras with scikit-learnのメモ](http://nobunaga.hatenablog.jp/entry/2017/10/17/071849)

※ 記事中に引用した文献や、Referenceで取り上げさせていただいた文献から
は、いずれも引用に許可を取っていないため、指摘された場合は文面から削除させていただきます
どれも大変参考にさせていただいております。ありがとうございます