今回は特徴量選択 (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 で分類して、モデルに組み込みの特徴量の可視化してみよう。
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 ,
'shuffle' : False ,
}
X, y = make_classification(**args)
clf = RandomForestClassifier(n_estimators=100 ,
random_state=42 )
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 次元が上位にきていることが確認できる。
上記の重要度を棒グラフにプロットしたものは以下のようになる。
RandomForest から得られる特徴量の重要度
この特徴量の重要度をそのまま使って特徴量選択をすることもできる。
たとえば上位 N 件を取り出した場合で交差検証のスコアを比較することが考えられる。
しかし、スコアの増加が小さな特徴量は単なるノイズかどうかの判断が難しい。
たとえば上記であれば「2」の特徴量は事前知識がなければノイズかどうか怪しいところ。
特徴量の重要度を Null Importance と比較する
そこで、特徴量の重要度をより細かく検証するために Null Importance と比較する。
前述したとおり Null Importance は目的変数をシャッフルして学習させたモデルから計算した特徴量の重要度となる。
Null Importance を基準として、元の特徴量の重要度がどれくらい大きいか比べることでノイズかどうかの判断がしやすくなる。
以下のサンプルコードでは Null Importance と元の特徴量の重要度を上位 10 件についてヒストグラムとしてプロットしている。
Null Importance は何回も計算するのが比較的容易なので、このようにヒストグラムで比較できる。
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 )
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]
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)):
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
すると、次のようなヒストグラムが得られる。
Null Importance のヒストグラムとの比較
先頭の 5 次元については Null Importance の分布から大きく離れている。
一方で、それ以降の特徴量の重要度は Null Importance の分布の中に入ってしまっていることがわかる。
ここから、先頭の 5 次元以降はノイズに近い特徴量であると判断する材料になる。
Null Importance を使って特徴量を選択してみる
それでは、実際に特徴量の選択までやってみよう。
以下のサンプルコードでは、Null Importance と元の特徴量の重要度の比率からスコアを計算している。
そして、スコアの大きさを元に上位 N% の特徴量を取り出して交差検証のスコアを確認している。
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)
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 ]
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
すると、以下のようなグラフが得られる。
Null Importance を用いて特徴量選択の可視化
今回は先頭 5 次元が重要とわかっている。
上記のグラフをみても上位 5% (= 5 次元) が選択されるまではスコアが順調に伸びている様子が確認できる。
しかし上位 1% (= 1 次元) までいくと推論に使える特徴量まで削ってしまいスコアが落ちていることがわかる。
ちなみに、実際には上記のように比較するのではなくスコアに一定の閾値を設けたり上位 N 件を決め打ちで選択してしまう場合も多いようだ。
LightGBM + Pandas を使った場合のサンプル
続いては、よくある例としてモデルに LightGBM を使って、データが Pandas のデータフレームの場合も確認しておこう。
サンプルコードは次のとおり。
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_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)
col_names = [f'col{i}' for i in range (n_cols)]
train_x = pd.DataFrame(X, columns=col_names)
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 ]
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% までいくとスコアが落ちていることがわかる。
データに Pandas モデルに LightGBM を使った場合
めでたしめでたし。
![ディジタル画像処理[改訂第二版]](https://m.media-amazon.com/images/I/41FCF736uWL._SL500_.jpg "ディジタル画像処理[改訂第二版]")
