Day 107 — 特徴量選択の高度手法 — SHAP・Boruta・PermutationImportance

📚 背景知識（読んでから問題へ）

🧩 「本物の特徴量」を見極める3つの検証手法

🎯 問題

9個の特徴量（本物の情報を持つ信号5個・純粋なノイズ3個・値の種類が非常に多い疑似ID列1個）の中から「本当に効いている特徴量」を4つの視点で見極める。

import numpy as np
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

RANDOM_STATE = 42
rng = np.random.RandomState(RANDOM_STATE)

# 5個の「本物」信号特徴量
X_base, y = make_classification(
    n_samples=3000, n_features=5, n_informative=5, n_redundant=0,
    n_clusters_per_class=2, class_sep=1.2, random_state=RANDOM_STATE
)
feat_names = [f"Signal{i+1}" for i in range(5)]
df = pd.DataFrame(X_base, columns=feat_names)

# 3個の純粋なノイズ特徴量（yと無関係な標準正規乱数）
for i in range(3):
    df[f"Noise{i+1}"] = rng.normal(0, 1, size=len(df))

# 1個の高カーディナリティな疑似ID列
df["HighCardID"] = rng.randint(0, 100000, size=len(df))

FEATURES = list(df.columns)
X = df[FEATURES].values

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, stratify=y, random_state=RANDOM_STATE
)

タスク

📐 組み込み重要度 vs Permutation Importance（各系列の最大値を100%として正規化）

実際に計算した値をもとに作図。オレンジ=組み込み重要度、青=Permutation Importance。ノイズ系4特徴量（右側4組）で、オレンジは薄く残るのに青はほぼ消えることに注目。

📊 実行結果の比較（4手法）

実際にコードを実行して得た値（乱数シード固定・環境により多少前後する）。

Test AUC

組み込み重要度 vs Permutation Importance（下位4特徴量に注目）

Boruta風シャドウ特徴量テスト（20試行中の勝利回数）

SHAP: mean(|SHAP値|)

SHAPのランキングはPermutation Importanceとほぼ同じ並びになった。Signal群（0.05〜0.11）とノイズ群（0.0016〜0.0041）の間に明確な断絶があり、判断がしやすい。

💡 ヒント

from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
from sklearn.metrics import roc_auc_score
import shap

# 1. ベースライン
rf = RandomForestClassifier(n_estimators=300, max_depth=6, random_state=42, n_jobs=-1)
rf.fit(X_train, y_train)
auc = roc_auc_score(y_test, rf.predict_proba(X_test)[:, 1])

# 2. 組み込み重要度
imp_builtin = pd.Series(rf.feature_importances_, index=FEATURES).sort_values(ascending=False)

# 3. Permutation Importance（testデータ）
perm = permutation_importance(rf, X_test, y_test, n_repeats=30, random_state=42, scoring="roc_auc")
imp_perm = pd.Series(perm.importances_mean, index=FEATURES).sort_values(ascending=False)

# 4. Boruta風シャドウ特徴量テスト
N_TRIALS = 20
hit_counts = pd.Series(0, index=FEATURES)
for trial in range(N_TRIALS):
    trial_rng = np.random.RandomState(1000 + trial)
    X_shadow = np.column_stack([trial_rng.permutation(X_train[:, j]) for j in range(X_train.shape[1])])
    X_combined = np.hstack([X_train, X_shadow])
    combined_names = FEATURES + [f"shadow_{f}" for f in FEATURES]
    rf_shadow = RandomForestClassifier(n_estimators=150, max_depth=6, random_state=trial_rng, n_jobs=-1)
    rf_shadow.fit(X_combined, y_train)
    imp_combined = pd.Series(rf_shadow.feature_importances_, index=combined_names)
    max_shadow = imp_combined[[f"shadow_{f}" for f in FEATURES]].max()
    for f in FEATURES:
        if imp_combined[f] > max_shadow:
            hit_counts[f] += 1

# 5. SHAP
explainer = shap.TreeExplainer(rf)
shap_values = explainer.shap_values(X_test)
sv = shap_values[1] if isinstance(shap_values, list) else shap_values[:, :, 1]
mean_abs_shap = pd.Series(np.abs(sv).mean(axis=0), index=FEATURES).sort_values(ascending=False)

✅ 模範解答

import numpy as np
import pandas as pd
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
from sklearn.metrics import roc_auc_score
import shap

RANDOM_STATE = 42
rng = np.random.RandomState(RANDOM_STATE)

# ====== データ準備 ======
X_base, y = make_classification(
    n_samples=3000, n_features=5, n_informative=5, n_redundant=0,
    n_clusters_per_class=2, class_sep=1.2, random_state=RANDOM_STATE
)
feat_names = [f"Signal{i+1}" for i in range(5)]
df = pd.DataFrame(X_base, columns=feat_names)
for i in range(3):
    df[f"Noise{i+1}"] = rng.normal(0, 1, size=len(df))
df["HighCardID"] = rng.randint(0, 100000, size=len(df))

FEATURES = list(df.columns)
X = df[FEATURES].values
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, stratify=y, random_state=RANDOM_STATE
)

# ====== 1. ベースラインモデル ======
rf = RandomForestClassifier(n_estimators=300, max_depth=6, random_state=RANDOM_STATE, n_jobs=-1)
rf.fit(X_train, y_train)
auc = roc_auc_score(y_test, rf.predict_proba(X_test)[:, 1])
print(f"Test AUC: {auc:.4f}")

# ====== 2. 組み込み重要度 ======
imp_builtin = pd.Series(rf.feature_importances_, index=FEATURES).sort_values(ascending=False)

# ====== 3. Permutation Importance（testデータ、n_repeats=30） ======
perm = permutation_importance(rf, X_test, y_test, n_repeats=30, random_state=RANDOM_STATE, scoring="roc_auc", n_jobs=-1)
imp_perm = pd.Series(perm.importances_mean, index=FEATURES).sort_values(ascending=False)

# ====== 4. Boruta風シャドウ特徴量テスト（20試行） ======
N_TRIALS = 20
hit_counts = pd.Series(0, index=FEATURES)
for trial in range(N_TRIALS):
    trial_rng = np.random.RandomState(1000 + trial)
    X_shadow = np.column_stack([trial_rng.permutation(X_train[:, j]) for j in range(X_train.shape[1])])
    X_combined = np.hstack([X_train, X_shadow])
    combined_names = FEATURES + [f"shadow_{f}" for f in FEATURES]
    rf_shadow = RandomForestClassifier(n_estimators=150, max_depth=6, random_state=trial_rng, n_jobs=-1)
    rf_shadow.fit(X_combined, y_train)
    imp_combined = pd.Series(rf_shadow.feature_importances_, index=combined_names)
    max_shadow = imp_combined[[f"shadow_{f}" for f in FEATURES]].max()
    for f in FEATURES:
        if imp_combined[f] > max_shadow:
            hit_counts[f] += 1

# ====== 5. SHAP ======
explainer = shap.TreeExplainer(rf)
shap_values = explainer.shap_values(X_test)
sv = shap_values[1] if isinstance(shap_values, list) else shap_values[:, :, 1]
mean_abs_shap = pd.Series(np.abs(sv).mean(axis=0), index=FEATURES).sort_values(ascending=False)

# ====== 6. 4手法の比較表 ======
comparison = pd.DataFrame({
    "builtin_rank": imp_builtin.rank(ascending=False).astype(int),
    "builtin_value": imp_builtin,
    "perm_rank": imp_perm.rank(ascending=False).astype(int),
    "perm_value": imp_perm,
    "boruta_hits/20": hit_counts,
    "shap_rank": mean_abs_shap.rank(ascending=False).astype(int),
    "shap_value": mean_abs_shap,
}).loc[FEATURES].sort_values("perm_rank")
print(comparison.round(4).to_string())

実際に実行して得た出力（RANDOM_STATE=42固定・環境により多少前後する）:

Test AUC: 0.9487

            builtin_rank  builtin_value  perm_rank  perm_value  boruta_hits/20  shap_rank  shap_value
Signal2                1         0.2159          1      0.1268              20          1      0.1140
Signal3                3         0.1929          2      0.1123              20          2      0.0922
Signal1                5         0.1530          3      0.0938              20          5      0.0520
Signal4                4         0.1826          4      0.0680              20          4      0.0713
Signal5                2         0.2054          5      0.0649              20          3      0.0853
HighCardID              7         0.0124          6      0.0003               6          7      0.0036
Noise1                  9         0.0120          7      0.0002               0          9      0.0016
Noise2                  8         0.0123          8      0.0000               1          8      0.0023
Noise3                  6         0.0136          9     -0.0001              10          6      0.0041

🪜 Step-by-Step 解説

1ベースラインモデルを学習する

rf = RandomForestClassifier(n_estimators=300, max_depth=6, random_state=42, n_jobs=-1)
rf.fit(X_train, y_train)
# Test AUC: 0.9487

2組み込み重要度を確認する — 一見もっともらしいが、罠がある

Signal2 0.2159, Signal5 0.2054, ... Noise3 0.0136, HighCardID 0.0124, Noise2 0.0123, Noise1 0.0120

3Permutation Importanceで「本当に効いているか」を検証する

HighCardID  : mean=0.0003
Noise1      : mean=0.0002
Noise2      : mean=0.0000
Noise3      : mean=-0.0001

4Boruta風シャドウ特徴量テストで統計的に判定する

Signal1〜5: 20/20回すべてでシャドウの最大値を上回った
Noise3: 10/20回（五分五分）
HighCardID: 6/20回
Noise2: 1/20回
Noise1: 0/20回

5SHAPでモデルの判断根拠を特徴量に分配する

Signal2 0.1140 … Noise3 0.0041, HighCardID 0.0036, Noise2 0.0023, Noise1 0.0016

64手法をまとめて「残すべき特徴量」を決める

🧮 数学・統計の補足（文系向け）

数式を見たら

🏆 Kaggleでの実践的な使い方

よく使われるコンペカテゴリ: ☑ 表形式データ（Tabular） / ☐ 自然言語処理（NLP） / ☐ 画像認識（CV） / ☑ 時系列（Time Series）

• 特徴量エンジニアリング後の絞り込み: 集計特徴量・Rolling/Lag特徴量（Day 099・100）を大量に生成した後、Permutation ImportanceかSHAPのmean(|SHAP値|)が低い特徴量を機械的に足切りし、上位N個だけで再学習することでCVスコアと学習速度の両方が改善することが多い
• Boruta（boruta_py）: 過去の上位解法ノートブックで「Borutaで数百の候補特徴量から本物の特徴量だけを抽出した」という報告は頻出。計算コストが高い（多数のモデルを繰り返し学習する）ため、GBDT系の高速版と組み合わせて使うのが実務的
• SHAPは「特徴量選択」と「モデル解釈」の二役: Kaggleの最終レポート（Solution write-up）ではSHAP summary plotで「どの特徴量がなぜ効いたか」を説明するのが定番。Public/Privateのシェイクアップ対策として、特定の1〜2個の特徴量にSHAP依存度が偏りすぎていないか確認するのにも使われる
• 相関の強い特徴量への注意: Permutation Importanceは、ある特徴量をシャッフルしても相関する別の特徴量が同じ情報を持っていれば「性能があまり落ちない＝重要度が低く見える」という弱点がある。強相関のグループはまとめてシャッフルする「grouped permutation importance」で対処するのが上位陣のテクニック

⚠️ よくある誤解・ミス

🔗 関連する重要な概念

• feature_importances_（Phase 2, Day 059）: 今日比較した「組み込み重要度」そのもの。今日学んだ弱点を踏まえて使う必要がある
• 不均衡データ対策（Day 106）: 評価指標にAUCやPR-AUCを使う判断は共通。今日のPermutation Importanceもscoring="roc_auc"を指定した
• SHAP値の解釈（Phase 4, Day 108予定）: 今日はグローバルなmean(|SHAP値|)のみ扱った。次回は1件ごとの予測（ローカル解釈）とSHAP summary plot・dependence plotを学ぶ

🚀 次のステップ

• 発展: Signal1〜5同士に人工的な相関を持たせたデータ（例: Signal1_dup = Signal1 + 少量のノイズ）を追加し、Permutation Importanceで元のSignal1とSignal1_dupの重要度がどう分散するかを確認してみる。相関特徴量が重要度を「薄める」現象を体験できる
• 次回予告: Day 108「SHAP値の解釈（グローバル/ローカル解釈・Kaggleでの活用）」 — 今日はモデル全体の傾向（グローバル）だけを見たが、次回は1件の予測ごとに「なぜこの予測になったか」を説明するローカル解釈（force plot・waterfall plot）とKaggleでの具体的な活用法に進む

📝 自己評価（解いた後に記入）