データ分析プロセス・可視化備忘録 : Home Credit Default Riskコンペのデータセットを用いて

はじめに

とある企業より「次回選考ではウチが保有するテーブルデータを用いて分析してみて！」と言われたので
データ分析のフロー、可視化方法など確認用の投稿です。

今回使用するデータセットはなんでもよかったのですが、
就職試験ではテーブルデータを用いて分類のモデルを作成して欲しいようなので
分類問題であるHome Credit Default Riskコンペのデータセットを使用します。

(*注)あくまで自己学習用の投稿なのでKaggleのカーネルを眺める方がよっぽど有意義であることを保証します。

目次

• 欠損値の補完
• 可視化
• パラメータの調整①(GridSearch)
• パラメータの調整②（Optuna）
• 特徴量エンジニアリング（Featuretools）
• おわりに

欠損値の補完

欠損とは文字通りデータが失われている状態です。
欠損値の発生要因は

- MCAR (Missing Completely At Random) / 欠損値が完全にランダムに発生
- MAR (Missing At Random) / 欠損は変数自体の値には依存していないが、他の変数には依存関係を持っている。
- MNAR (Missing Not At Random) / 分析における他の変数及び、変数自体に対して依存関係を持っている。

上記3タイプに分けられるよう。
そして上記3つのどの原因をもってして欠損しているのか、それによって必要な処理（補完or削除）が必要なのかが異なるよう。

[参考]
欠損データ分析 (missing data analysis)
データ解析における欠損値対応フロー

とはいえ、就職試験も近いので、先を急ぐ(笑)
今回は下記の特徴量を用いる。

['SK_ID_CURR','EXT_SOURCE_1','EXT_SOURCE_2','EXT_SOURCE_3','DAYS_BIRTH',
'AMT_REQ_CREDIT_BUREAU_HOUR','AMT_REQ_CREDIT_BUREAU_DAY','AMT_REQ_CREDIT_BUREAU_WEEK',
'AMT_REQ_CREDIT_BUREAU_MON','AMT_REQ_CREDIT_BUREAU_QRT','AMT_REQ_CREDIT_BUREAU_YEAR']

また、今回はあくまで手法や流れの確認なので、dfをサイズダウン。

# shape = (10000, 10)
df = df.iloc[0:10000,:]

欠損がちらほら...欠損を埋める関数を用意。
この関数は、

• 欠損率が任意の割合を超える特徴量の削除
• type objectの特徴量は最頻値で補完
• type int / float はランダムフォレストで補完
• one-hot encoding

ランダムフォレストは決定木を複数用いて、決定木単体よりも精度をあげるもの。
肝心な決定木はというとある条件(ジニ不純度, エントロピー, 分類誤差)に従ってデータを分類していくモデル。
欠損値補完においてランダムフォレストを用いることは確率的に同分類とされる他IDのデータを参考に欠損を補完をするという理解（適当w）

from missingpy import MissForest
from tqdm import tqdm

def drop_fill(df, thredhold):
    """
    note : 欠損率がthreshold以上の列は削除
            dtypesがobjectのものはmodeで欠損値を埋め
            dtypesがint / floatはRandomForestで埋める
    ---------------------
    df : df
    thredhold : float
    ---------------------
    attribute : df
    """
    
    # 欠損率が任意の数値未満の列を抽出
    object_columns = []
    for i in tqdm(df.columns):
        if df[i].dtypes == 'O' and (df[i].isnull().sum() / len(df[i])) < thredhold:
            object_columns.append(i)
    
    # 欠損値が任意の数値未満の列を抽出
    float_int_columns = []
    for i in tqdm(df.columns):
        if df[i].dtypes == 'float64' or 'int64' and (df[i].isnull().sum() / len(df[i])) < thredhold:
            float_int_columns.append(i)
    
    # objectは最頻値で補完
    df_object = pd.DataFrame(df[object_columns]).fillna(df[object_columns].mode().iloc[0])
    df_int_float = pd.DataFrame(df[float_int_columns])
    df_finalized = pd.concat([df_object, df_int_float], axis = 1)
    
    # onehot encoding
    df_onehot = pd.get_dummies(df_finalized[df_finalized.columns])
    
    # float, intはrandom forestで補完
    imputer = MissForest()
    df_imputed = imputer.fit_transform(df_onehot)

    return pd.DataFrame(df_imputed, columns=df_onehot.columns)

可視化

ここはIrisデータを元に可視化のコードを羅列するスタイル (手抜きw)

import seaborn as sns

# < ヒストグラム付散布図 >
# 通常の散布図に比べて、ヒストグラムを載せることで、最頻値や分布の状況を簡単に把握することができる
sns.jointplot(x='sepal length (cm)', y='petal length (cm)', data=df)

# < 箱ひげ図 >
# データの四分位点などを同時に表現できるグラフ
sns.boxplot(x="species", y="sepal length (cm)", data=df)

# < バイオリンプロット >
# 中央値、四分位点のほかに、データの分布密度も同時に確認できるプロット
sns.violinplot(x="species", y="sepal length (cm)", data=df)

# < 相関図 >
# 散布図行列、二変数間の相関を示す。
sns.pairplot(df, hue="Species", size=2.5)

# < ヒートマップ >
# method='spearman'でスピアマン、'kendall'でケンドールも指定可能
corr_mat = df.corr(method='pearson')
sns.heatmap(corr_mat,
            vmin=-1.0,
            vmax=1.0,
            center=0,
            annot=True, # True:格子の中に値を表示
            fmt='.1f',
            xticklabels=corr_mat.columns.values,
            yticklabels=corr_mat.columns.values
           )

[参考]
pythonでデータを可視化したいならseabornを使おう！
seaborn の heatmap で美しく可視化
seaborn: statistical data visualization

パラメータの調整①(GridSearch)

モデルを作成するにあたり、アルゴリズムのパラメーターを設定する必要が。
グリッドサーチはシンプルで全てのパラメータを試す、総当たり方式なので当然探索に時間がかかる。

from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.ensemble import RandomForestClassifier

def rf_pipeline(X_train, y_train):
    
    steps = [('SC', StandardScaler()),
             ('PCA', PCA()),
             ('RF', RandomForestClassifier())
            ]

    params = {'PCA__whiten':['True', 'False'],
              'RF__max_depth':[3, 10],
              'RF__n_estimators':[50, 100],
              'RF__min_samples_split':[2, 10],
              'RF__min_samples_leaf':[1, 10],
              'RF__bootstrap':[True, False],
              'RF__criterion':["gini", "entropy"]
              
             }

    pipeline = Pipeline(steps)
    gs = GridSearchCV(pipeline, params, cv=10)
    gs.fit(X_train, y_train)

    return gs.best_params_, gs.best_estimator_, gs.best_score_

CPU times: user 39min 5s, sys: 6.99 s, total: 39min 12s
Wall time: 25min 57s

[参考]
Python: 無名数化によるデータの前処理
sklearn.ensemble.RandomForestClassifier¶

パラメータの調整②（Optuna）

GridSearchではパラメータ探索に25分も費やしてしまった。
就職試験で与えられる時間は3時間のみなので、これでは命取りである。

__Optuna__を使用してみる。

Optuna はハイパーパラメータの最適化を自動化するためのソフトウェアフレームワーク。

Optuna は次の試行で試すべきハイパーパラメータの値を決めるために、完了している試行の履歴を用いている。

完了している試行の履歴に基づき、有望そうな領域を推定し、その領域の値を実際に試すということを繰り返す。

Tree-structured Parzen Estimator というベイズ最適化アルゴリズムの一種を用いている。

[参考]
機械学習のためのベイズ最適化入門
hyperoptって何してんの？
GridSearchCVはもう古い！Optunaでサポートベクターマシンもラクラクチューニング

まずOptunaが最適化に使うための関数を定義する必要がある。

import optuna
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_validate
 

def objective(train_X, test_X, train_y, test_y, trial):

    max_depth = trial.suggest_int('max_depth', 3, 15) 
    n_estimators = trial.suggest_int('n_estimators', 0,300)
    min_samples_split = trial.suggest_int('min_samples_split', 2,10)
    min_samples_leaf = trial.suggest_int('min_samples_leaf', 1,10) 
    bootstrap = trial.suggest_categorical('bootstrap', [True, False])
    criterion = trial.suggest_categorical('criterion', ["gini", "entropy"])

    model = RandomForestClassifier(max_depth=max_depth,
                                   n_estimators=n_estimators,
                                   min_samples_split=min_samples_split,
                                   min_samples_leaf=min_samples_leaf,
                                   bootstrap=bootstrap,
                                   criterion=criterion
                                  )
 
    model.fit(train_X, train_y)
    y_pred = model.predict(test_X)
    
    kf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
    scores = cross_validate(model, X=train_X, y=train_y, cv=kf)
    
    return 1.0 - scores['test_score'].mean()

from sklearn.model_selection import train_test_split
from functools import partial
train_X, test_X, train_y, test_y = train_test_split(df_X_new, df_y)
 
f = partial(objective, train_X, test_X, train_y, test_y)
study = optuna.create_study()
optuna.logging.disable_default_handler()
study.optimize(f, n_trials=10)

探索範囲（パラメータ）が若干異なるとはいえ、以下の通り驚愕の時短。

CPU times: user 1min 19s, sys: 632 ms, total: 1min 19s
Wall time: 1min 20s

特徴量エンジニアリング（Featuretools）

現在のところ精度は

0.700845876776329

精度向上のために特徴量を生成します。
featuretools先生の出番。

Featuretools is a framework to perform automated feature engineering. It excels at transforming transactional and relational datasets into feature matrices for machine learning.

自動で特徴量いっぱい作ってくれる。w
Will Koehrsen氏のカーネルは大変勉強になった。

[参考]
Feature-Engineeringのリンク集めてみた

import featuretools as ft

es = ft.EntitySet(id = 'homecredit')
es.entity_from_dataframe(entity_id = 'homecredit', dataframe = df_X_new, index = 'SK_ID_CURR')

agg_primitives = ["mean", "sum", "mode"]

feature_matrix, feature_defs = ft.dfs(entityset = es, target_entity = 'homecredit' ,
                                      agg_primitives=agg_primitives, trans_primitives = ['add_numeric', 'multiply_numeric'])

feature_train_X, feature_test_X, train_y_, test_y_= train_test_split(feature_matrix, df_y)
# feature_matrix.shape = (10000, 1540)

特徴量生成したdfでパラメータ探索を行い、精度算定したところ精度が上がった。

0.7327325962618401

おわりに

今回は就職試験対策にデータ分析フローを再度確認した。
3時間という限られた時間の中で結果を出すために、特徴量生成、パラメータ探索を高速で行えるよう、
featuretools, optunaを使用した。

精度が上がったのでひとまずめでたしめでたし。