今回は特徴量選択 (Feature Selection) の手法のひとつとして使われることのある Null Importance を試してみる。 Null Importance というのは、目的変数をシャッフルして意味がなくなった状態で学習させたモデルから得られる特徴量の重要度を指す。 では、それを使ってどのように特徴量選択をするかというと、シャッフルしなかったときの重要度との比率をスコアとして計算する。 もし、シャッフルしたときの重要度が元となった重要度よりも小さくなっていれば、スコアは大きくなって特徴量に意味があるとみなせる。 一方で、シャッフルしたときの重要度が元とさほど変わらなければ、スコアは小さくなってその特徴量は単なるノイズに近い存在と判断できる。 あとはスコアに一定の閾値を設けたり、上位 N 件を取り出すことで特徴量選択ができるようだ。
今回使った環境は次のとおり。
$ sw_vers ProductName: Mac OS X ProductVersion: 10.15.6 BuildVersion: 19G73 $ python -V Python 3.8.3
下準備
下準備として、あらかじめ今回使うパッケージをインストールしておく。
$ pip install scikit-learn lightgbm matplotlib pandas tqdm
題材とするデータについて
題材としては scikit-learn が生成するダミーデータを用いる。
sklearn.datasets.make_classification() 関数を使うと、データの次元数や推論の難易度などを調整してダミーデータが作れる。
以下のサンプルコードでは二値分類のダミーデータを生成している。 このダミーデータは全体で 100 次元あるものの、先頭の 5 次元だけが推論する上で意味のあるデータとなっている。 試しにダミーデータを RandomForest で分類して、モデルに組み込みの特徴量の可視化してみよう。
#!/usr/bin/env python # -*- coding: utf-8 -*- import numpy as np from sklearn.datasets import make_classification from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import StratifiedKFold from sklearn.model_selection import cross_validate from matplotlib import pyplot as plt def main(): # 疑似的な教師信号を作るためのパラメータ args = { # データ点数 'n_samples': 1_000, # 次元数 'n_features': 100, # その中で意味のあるもの 'n_informative': 5, # 重複や繰り返しはなし 'n_redundant': 0, 'n_repeated': 0, # タスクの難易度 'class_sep': 0.65, # 二値分類問題 'n_classes': 2, # 生成に用いる乱数 'random_state': 42, # 特徴の順序をシャッフルしない (先頭の次元が informative になる) 'shuffle': False, } # ノイズを含んだ教師データを作る X, y = make_classification(**args) # 分類器にランダムフォレストを使う clf = RandomForestClassifier(n_estimators=100, random_state=42) # Stratified 5-Fold CV + ROC AUC でモデルを検証する folds = StratifiedKFold(n_splits=5, shuffle=True, random_state=42, ) # すべての次元を使った場合 cv_all_result = cross_validate(clf, X, y, cv=folds, return_estimator=True, scoring='roc_auc', ) # 意味のある特徴量だけを使った場合 cv_ideal_result = cross_validate(clf, # 使う次元を先頭だけに絞り込む X[:, :args['n_informative']], y, return_estimator=False, scoring='roc_auc', ) # それぞれの状況でのメトリックのスコア print('All used AUC:', np.mean(cv_all_result['test_score'])) print('Ideal AUC:', np.mean(cv_ideal_result['test_score'])) # 学習済みモデルを取り出す clfs = cv_all_result['estimator'] # 特徴量の重要度を取り出す importances = [clf.feature_importances_ for clf in clfs] mean_importances = np.mean(importances, axis=0) std_importances = np.std(importances, axis=0) sorted_indices = np.argsort(mean_importances)[::-1] # 重要度の高い特徴量を表示する MAX_TOP_N = 10 rank_n = min(X.shape[1], MAX_TOP_N) print('Feature importance ranking (TOP {rank_n})'.format(rank_n=rank_n)) for rank, idx in enumerate(sorted_indices[:rank_n], start=1): params = { 'rank': rank, 'idx': idx, 'importance': mean_importances[idx] } print('{rank}. feature {idx:02d}: {importance}'.format(**params)) # 特徴量の重要度を可視化する plt.figure(figsize=(6, 8)) plt.barh(range(rank_n), mean_importances[sorted_indices][:rank_n][::-1], color='g', ecolor='r', yerr=std_importances[sorted_indices][:rank_n][::-1], align='center') plt.yticks(range(rank_n), sorted_indices[:rank_n][::-1]) plt.ylabel('Features') plt.xlabel('Importance') plt.grid() plt.show() if __name__ == '__main__': main()
上記を実行してみよう。 はじめに、すべての次元を使って学習した場合と、実際に意味のある次元 (先頭 5 次元) だけを使った場合の ROC AUC を出力している。 そして、それに続いて値の大きさで降順ソートした特徴量の重要度が表示される。
$ python rfimp.py All used AUC: 0.8748692786445511 Ideal AUC: 0.947211653938503 Feature importance ranking (TOP 10) 1. feature 00: 0.0843898282531266 2. feature 03: 0.05920082172020726 3. feature 01: 0.05713024480668112 4. feature 04: 0.05027373664878948 5. feature 02: 0.017385737849956347 6. feature 55: 0.010379957110544775 7. feature 52: 0.009663802073627043 8. feature 87: 0.00907575944742235 9. feature 34: 0.008999985108968961 10. feature 82: 0.008999816679827738
先頭の 5 次元だけを使った場合の方が ROC AUC のスコアが高くなっていることがわかる。 また、RandomForest のモデルに組み込まれている特徴量の重要度を見ても先頭の 5 次元が上位にきていることが確認できる。 上記の重要度を棒グラフにプロットしたものは以下のようになる。
この特徴量の重要度をそのまま使って特徴量選択をすることもできる。 たとえば上位 N 件を取り出した場合で交差検証のスコアを比較することが考えられる。 しかし、スコアの増加が小さな特徴量は単なるノイズかどうかの判断が難しい。 たとえば上記であれば「2」の特徴量は事前知識がなければノイズかどうか怪しいところ。
特徴量の重要度を Null Importance と比較する
そこで、特徴量の重要度をより細かく検証するために Null Importance と比較する。 前述したとおり Null Importance は目的変数をシャッフルして学習させたモデルから計算した特徴量の重要度となる。 Null Importance を基準として、元の特徴量の重要度がどれくらい大きいか比べることでノイズかどうかの判断がしやすくなる。
以下のサンプルコードでは Null Importance と元の特徴量の重要度を上位 10 件についてヒストグラムとしてプロットしている。 Null Importance は何回も計算するのが比較的容易なので、このようにヒストグラムで比較できる。
#!/usr/bin/env python # -*- coding: utf-8 -*- import logging import sys from itertools import chain from tqdm import tqdm import numpy as np from sklearn.datasets import make_classification from sklearn.ensemble import RandomForestClassifier from matplotlib import pyplot as plt from sklearn.model_selection import StratifiedKFold from sklearn.model_selection import cross_validate LOGGER = logging.getLogger(__name__) def feature_importance(X, y): """特徴量の重要度を交差検証で計算する""" clf = RandomForestClassifier(n_estimators=100, random_state=42) folds = StratifiedKFold(n_splits=5, shuffle=True, random_state=42, ) cv_result = cross_validate(clf, X, y, cv=folds, return_estimator=True, scoring='roc_auc', n_jobs=-1, ) clfs = cv_result['estimator'] importances = [clf.feature_importances_ for clf in clfs] mean_importances = np.mean(importances, axis=0) return mean_importances def main(): # 結構時間がかかるのでログを出す logging.basicConfig(level=logging.INFO, stream=sys.stderr) args = { 'n_samples': 1_000, 'n_features': 100, 'n_informative': 5, 'n_redundant': 0, 'n_repeated': 0, 'class_sep': 0.65, 'n_classes': 2, 'random_state': 42, 'shuffle': False, } X, y = make_classification(**args) # ベースとなる特徴量の重要度を計算する LOGGER.info('Starting base importance calculation') base_importance = feature_importance(X, y) # データのシャッフルに再現性をもたせる np.random.seed(42) # Null Importance を何度か計算する LOGGER.info('Starting null importance calculation') TRIALS_N = 20 null_importances = [] for _ in tqdm(range(TRIALS_N)): # 目的変数をシャッフルする y_permuted = np.random.permutation(y) # シャッフルした状態で特徴量の重要度を計算する null_importance = feature_importance(X, y_permuted) null_importances.append(null_importance) null_importances = np.array(null_importances) # 列と行を入れ替える transposed_null_importances = null_importances.transpose(1, 0) sorted_indices = np.argsort(base_importance)[::-1] sorted_base_importance = base_importance[sorted_indices] sorted_null_importance = transposed_null_importances[sorted_indices] # ベースとなる特徴量の重要度の上位と Null Importance を可視化する LOGGER.info('Starting visualization') HIST_ROWS = 3 HIST_COLS = 3 fig, axes = plt.subplots(HIST_ROWS, HIST_COLS, figsize=(8, 18)) for index, ax in enumerate(chain.from_iterable(axes)): # Null Importance をヒストグラムにする col_null_importance = sorted_null_importance[index] ax.hist(col_null_importance, label='Null Importance', color='b') # ベースとなる特徴量の重要度の場所に縦線を引く col_base_importance = sorted_base_importance[index] ax.axvline(col_base_importance, label='Base Importance', color='r') # グラフの体裁 ax.set_xlabel('Importance') ax.set_ylabel('Frequency') ax.set_title(f'Feature: {sorted_indices[index]}') ax.legend() plt.show() if __name__ == '__main__': main()
上記を実行してみよう。
$ python nullimphist.py
すると、次のようなヒストグラムが得られる。
先頭の 5 次元については Null Importance の分布から大きく離れている。 一方で、それ以降の特徴量の重要度は Null Importance の分布の中に入ってしまっていることがわかる。 ここから、先頭の 5 次元以降はノイズに近い特徴量であると判断する材料になる。
Null Importance を使って特徴量を選択してみる
それでは、実際に特徴量の選択までやってみよう。 以下のサンプルコードでは、Null Importance と元の特徴量の重要度の比率からスコアを計算している。 そして、スコアの大きさを元に上位 N% の特徴量を取り出して交差検証のスコアを確認している。
#!/usr/bin/env python # -*- coding: utf-8 -*- import logging import sys from tqdm import tqdm import numpy as np from sklearn.datasets import make_classification from sklearn.ensemble import RandomForestClassifier from matplotlib import pyplot as plt from sklearn.model_selection import StratifiedKFold from sklearn.model_selection import cross_validate LOGGER = logging.getLogger(__name__) def _cross_validate(X, y): clf = RandomForestClassifier(n_estimators=100, random_state=42) folds = StratifiedKFold(n_splits=5, shuffle=True, random_state=42, ) cv_result = cross_validate(clf, X, y, cv=folds, return_estimator=True, scoring='roc_auc', n_jobs=-1, ) return cv_result def cv_mean_feature_importance(X, y): """特徴量の重要度を交差検証で計算する""" cv_result = _cross_validate(X, y) clfs = cv_result['estimator'] importances = [clf.feature_importances_ for clf in clfs] mean_importances = np.mean(importances, axis=0) return mean_importances def cv_mean_test_score(X, y): """OOF Prediction のメトリックを交差検証で計算する""" cv_result = _cross_validate(X, y) mean_test_score = np.mean(cv_result['test_score']) return mean_test_score def main(): logging.basicConfig(level=logging.INFO, stream=sys.stderr) n_cols = 100 args = { 'n_samples': 1_000, 'n_features': n_cols, 'n_informative': 5, 'n_redundant': 0, 'n_repeated': 0, 'class_sep': 0.65, 'n_classes': 2, 'random_state': 42, 'shuffle': False, } X, y = make_classification(**args) columns = np.arange(n_cols) LOGGER.info('Starting base importance calculation') base_importance = cv_mean_feature_importance(X, y) np.random.seed(42) LOGGER.info('Starting null importance calculation') TRIALS_N = 20 null_importances = [] for _ in tqdm(range(TRIALS_N)): y_permuted = np.random.permutation(y) null_importance = cv_mean_feature_importance(X, y_permuted) null_importances.append(null_importance) null_importances = np.array(null_importances) # Null Importance を用いた特徴量選択の一例 # Base Importance と Null Importance の値の比をスコアとして計算する criterion_percentile = 50 # 基準となるパーセンタイル (ここでは中央値) percentile_null_imp = np.percentile(null_importances, criterion_percentile, axis=0) null_imp_score = base_importance / (percentile_null_imp + 1e-6) # スコアの大きさで降順ソートする sorted_indices = np.argsort(null_imp_score)[::-1] # 上位 N% の特徴量を使って性能を比較してみる use_feature_importance_top_percentages = [100, 75, 50, 25, 10, 5, 1] mean_test_scores = [] selected_cols_len = [] for percentage in use_feature_importance_top_percentages: # スコアが上位のカラムを取り出す sorted_columns = columns[sorted_indices] num_of_features = int(n_cols * percentage / 100) selected_cols = sorted_columns[:num_of_features] X_selected = X[:, selected_cols] LOGGER.info(f'Null Importance score TOP {percentage}%') LOGGER.info(f'selected features: {selected_cols}') LOGGER.info(f'selected feature len: {len(selected_cols)}') selected_cols_len.append(len(selected_cols)) mean_test_score = cv_mean_test_score(X_selected, y) LOGGER.info(f'mean test score: {mean_test_score}') mean_test_scores.append(mean_test_score) # 結果を可視化する _, ax1 = plt.subplots(figsize=(8, 4)) ax1.plot(mean_test_scores, color='b', label='mean test score') ax1.set_xlabel('percentile') ax1.set_ylabel('mean test score') ax1.legend() ax2 = ax1.twinx() ax2.plot(selected_cols_len, color='r', label='selected features len') ax2.set_ylabel('selected features len') ax2.legend() plt.xticks(range(len(use_feature_importance_top_percentages)), use_feature_importance_top_percentages) plt.show() if __name__ == '__main__': main()
上記を実行してみよう。
$ python nullimpselect.py INFO:__main__:Starting base importance calculation INFO:__main__:Starting null importance calculation 100%|███████████████████████████████████████████| 20/20 [00:56<00:00, 2.82s/it] INFO:__main__:Null Importance score TOP 100% INFO:__main__:selected features: [ 0 3 1 4 2 55 13 52 34 87 67 18 29 56 82 94 23 75 58 89 54 47 97 61 85 93 64 44 43 45 63 7 37 38 27 41 26 12 72 10 51 49 96 50 22 53 81 66 78 69 60 79 16 90 19 80 24 17 36 71 15 8 39 20 98 40 11 86 84 6 74 73 9 65 31 35 92 91 68 48 25 99 88 76 42 28 5 70 30 95 33 83 14 59 57 32 21 46 77 62] INFO:__main__:selected feature len: 100 INFO:__main__:mean test score: 0.8658391949639143 INFO:__main__:Null Importance score TOP 75% INFO:__main__:selected features: [ 0 3 1 4 2 55 13 52 34 87 67 18 29 56 82 94 23 75 58 89 54 47 97 61 85 93 64 44 43 45 63 7 37 38 27 41 26 12 72 10 51 49 96 50 22 53 81 66 78 69 60 79 16 90 19 80 24 17 36 71 15 8 39 20 98 40 11 86 84 6 74 73 9 65 31] INFO:__main__:selected feature len: 75 INFO:__main__:mean test score: 0.8774203740902301 INFO:__main__:Null Importance score TOP 50% INFO:__main__:selected features: [ 0 3 1 4 2 55 13 52 34 87 67 18 29 56 82 94 23 75 58 89 54 47 97 61 85 93 64 44 43 45 63 7 37 38 27 41 26 12 72 10 51 49 96 50 22 53 81 66 78 69] INFO:__main__:selected feature len: 50 INFO:__main__:mean test score: 0.909876116663287 INFO:__main__:Null Importance score TOP 25% INFO:__main__:selected features: [ 0 3 1 4 2 55 13 52 34 87 67 18 29 56 82 94 23 75 58 89 54 47 97 61 85] INFO:__main__:selected feature len: 25 INFO:__main__:mean test score: 0.932729764573096 INFO:__main__:Null Importance score TOP 10% INFO:__main__:selected features: [ 0 3 1 4 2 55 13 52 34 87] INFO:__main__:selected feature len: 10 INFO:__main__:mean test score: 0.9513031915436441 INFO:__main__:Null Importance score TOP 5% INFO:__main__:selected features: [0 3 1 4 2] INFO:__main__:selected feature len: 5 INFO:__main__:mean test score: 0.9629145987828075 INFO:__main__:Null Importance score TOP 1% INFO:__main__:selected features: [0] INFO:__main__:selected feature len: 1 INFO:__main__:mean test score: 0.6537074505769904
すると、以下のようなグラフが得られる。
今回は先頭 5 次元が重要とわかっている。 上記のグラフをみても上位 5% (= 5 次元) が選択されるまではスコアが順調に伸びている様子が確認できる。 しかし上位 1% (= 1 次元) までいくと推論に使える特徴量まで削ってしまいスコアが落ちていることがわかる。 ちなみに、実際には上記のように比較するのではなくスコアに一定の閾値を設けたり上位 N 件を決め打ちで選択してしまう場合も多いようだ。
LightGBM + Pandas を使った場合のサンプル
続いては、よくある例としてモデルに LightGBM を使って、データが Pandas のデータフレームの場合も確認しておこう。 サンプルコードは次のとおり。
#!/usr/bin/env python3 # -*- coding: utf-8 -*- import logging import sys import numpy as np import pandas as pd import lightgbm as lgb from sklearn.metrics import roc_auc_score from tqdm import tqdm from sklearn.datasets import make_classification from sklearn.model_selection import StratifiedKFold from matplotlib import pyplot as plt LOGGER = logging.getLogger(__name__) class ModelExtractionCallback(object): """see: https://blog.amedama.jp/entry/lightgbm-cv-model""" def __init__(self): self._model = None def __call__(self, env): self._model = env.model def _assert_called_cb(self): if self._model is None: raise RuntimeError('callback has not called yet') @property def cvbooster(self): self._assert_called_cb() return self._model def _cross_validate(train_x, train_y, folds): """LightGBM で交差検証する関数""" lgb_train = lgb.Dataset(train_x, train_y) model_extraction_cb = ModelExtractionCallback() callbacks = [ model_extraction_cb, ] lgbm_params = { 'objective': 'binary', 'metric': 'auc', 'first_metric_only': True, 'verbose': -1, } lgb.cv(lgbm_params, lgb_train, folds=folds, num_boost_round=1_000, early_stopping_rounds=100, callbacks=callbacks, verbose_eval=20, ) return model_extraction_cb.cvbooster def _predict_oof(cv_booster, train_x, train_y, folds): """学習済みモデルから Out-of-Fold Prediction を求める""" split = folds.split(train_x, train_y) fold_mappings = zip(split, cv_booster.boosters) oof_y_pred = np.zeros_like(train_y, dtype=float) for (_, val_index), booster in fold_mappings: val_train_x = train_x.iloc[val_index] y_pred = booster.predict(val_train_x, num_iteration=cv_booster.best_iteration) oof_y_pred[val_index] = y_pred return oof_y_pred def cv_mean_feature_importance(train_x, train_y, folds): """交差検証したモデルを使って特徴量の重要度を計算する""" cv_booster = _cross_validate(train_x, train_y, folds) importances = cv_booster.feature_importance(importance_type='gain') mean_importance = np.mean(importances, axis=0) return mean_importance def cv_mean_test_score(train_x, train_y, folds): """交差検証で OOF Prediction の平均スコアを求める""" cv_booster = _cross_validate(train_x, train_y, folds) # OOF Pred を取得する oof_y_pred = _predict_oof(cv_booster, train_x, train_y, folds) test_score = roc_auc_score(train_y, oof_y_pred) return test_score def main(): logging.basicConfig(level=logging.INFO, stream=sys.stderr) n_cols = 100 args = { 'n_samples': 1_000, 'n_features': n_cols, 'n_informative': 5, 'n_redundant': 0, 'n_repeated': 0, 'class_sep': 0.65, 'n_classes': 2, 'random_state': 42, 'shuffle': False, } X, y = make_classification(**args) # Pandas のデータフレームにする col_names = [f'col{i}' for i in range(n_cols)] train_x = pd.DataFrame(X, columns=col_names) # インデックスが 0 から始まる連番とは限らないのでこういうチェックを入れた方が良い train_x.index = train_x.index * 10 + 10 train_y = pd.Series(y, name='target') folds = StratifiedKFold(n_splits=5, shuffle=True, random_state=42, ) LOGGER.info('Starting base importance calculation') base_importance = cv_mean_feature_importance(train_x, train_y, folds) LOGGER.info('Starting null importance calculation') TRIALS_N = 20 null_importances = [] for _ in tqdm(range(TRIALS_N)): train_y_permuted = np.random.permutation(train_y) null_importance = cv_mean_feature_importance(train_x, train_y_permuted, folds) null_importances.append(null_importance) null_importances = np.array(null_importances) criterion_percentile = 50 percentile_null_imp = np.percentile(null_importances, criterion_percentile, axis=0) null_imp_score = base_importance / (percentile_null_imp + 1e-6) sorted_indices = np.argsort(null_imp_score)[::-1] # 上位 N% の特徴量を使って性能を比較してみる use_feature_importance_top_percentages = [100, 75, 50, 25, 10, 5, 1] mean_test_scores = [] percentile_selected_cols = [] for percentage in use_feature_importance_top_percentages: sorted_columns = train_x.columns[sorted_indices] num_of_features = int(n_cols * percentage / 100) selected_cols = sorted_columns[:num_of_features] selected_train_x = train_x[selected_cols] LOGGER.info(f'Null Importance score TOP {percentage}%') LOGGER.info(f'selected features: {list(selected_cols)}') LOGGER.info(f'selected feature len: {len(selected_cols)}') percentile_selected_cols.append(selected_cols) mean_test_score = cv_mean_test_score(selected_train_x, train_y, folds) LOGGER.info(f'mean test_score: {mean_test_score}') mean_test_scores.append(mean_test_score) _, ax1 = plt.subplots(figsize=(8, 4)) ax1.plot(mean_test_scores, color='b', label='mean test score') ax1.set_xlabel('Importance TOP n%') ax1.set_ylabel('mean test score') ax1.legend() ax2 = ax1.twinx() ax2.plot([len(cols) for cols in percentile_selected_cols], color='r', label='selected features len') ax2.set_ylabel('selected features len') ax2.legend() plt.xticks(range(len(use_feature_importance_top_percentages)), use_feature_importance_top_percentages) plt.show() if __name__ == '__main__': main()
上記を実行してみよう。
$ python lgbmnullimp.py
すると、次のようなヒストグラムが得られる。 こちらも、上位 5% を選択するまで徐々にスコアが上がって、上位 1% までいくとスコアが落ちていることがわかる。
めでたしめでたし。
- 作者:もみじあめ
- 発売日: 2020/02/29
- メディア: Kindle版
![ディジタル画像処理[改訂第二版]](https://m.media-amazon.com/images/I/41FCF736uWL._SL500_.jpg "ディジタル画像処理[改訂第二版]")
