Day 067 — 特徴量重要度と SHAP — GBDTの予測根拠を解釈する

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

なぜ解釈性が必要か

LightGBM / XGBoost は高精度だが「なぜその予測をしたか」が分からない。
特徴量重要度と SHAP を使うことで次の3つが可能になる：

• 特徴量エンジニアリングの評価: 新しく作った特徴量が有効かを定量確認
• データリークの検出: 本来重要でないはずの特徴量が上位に来たら危険信号
• コンペ考察の投稿: 「なぜこの特徴量が重要か」を説明してノートブックを充実させる

📊 特徴量重要度の3種類

🔮 SHAPの概念

予測値の分解式

feature_importances_ vs SHAP

📈 特徴量重要度の典型的な分布（Titanic 風データ）

gain ベースの重要度（上位10特徴量）

gain vs split の差異（典型例）

🌊 Waterfall プロットの読み方（1サンプル例）

例: テストデータの1人目「女性・1等・31歳・運賃高め」→ 生存確率 0.92

🗂️ データスキーマ（Titanic 風 n=5,000）

📝 問題

セットアップコード（最初に実行）

📋 コピーimport numpy as np
import pandas as pd
import lightgbm as lgb
import shap
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import roc_auc_score
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)
n = 5000
data = pd.DataFrame({
    'Pclass':      np.random.choice([1, 2, 3], n, p=[0.24, 0.21, 0.55]),
    'Sex':         np.random.choice(['male', 'female'], n, p=[0.65, 0.35]),
    'Age':         np.random.normal(29.7, 14.5, n).clip(1, 80),
    'SibSp':       np.random.choice(range(9), n, p=[0.68,0.23,0.04,0.02,0.01,0.01,0.005,0.005,0.01]),
    'Parch':       np.random.choice(range(7), n, p=[0.76,0.13,0.09,0.005,0.005,0.005,0.005]),
    'Fare':        np.abs(np.random.normal(32, 50, n)),
    'Embarked':    np.random.choice(['S', 'C', 'Q'], n, p=[0.72, 0.19, 0.09]),
    'Cabin_type':  np.random.choice(['A','B','C','D','E','None'], n,
                                    p=[0.05,0.08,0.12,0.1,0.07,0.58]),
    'FamilySize':  None,
    'IsAlone':     None,
})
data['FamilySize'] = data['SibSp'] + data['Parch'] + 1
data['IsAlone'] = (data['FamilySize'] == 1).astype(int)
data['Survived'] = (
    (data['Sex'] == 'female') * 0.55 + (data['Pclass'] == 1) * 0.25 +
    (data['Age'] < 16) * 0.2 + np.random.normal(0, 0.12, n)
) > 0.4
data['Survived'] = data['Survived'].astype(int)

for col in ['Sex', 'Embarked', 'Cabin_type']:
    data[col] = LabelEncoder().fit_transform(data[col])

features = ['Pclass', 'Sex', 'Age', 'SibSp', 'Parch', 'Fare',
            'Embarked', 'Cabin_type', 'FamilySize', 'IsAlone']
X = data[features]
y = data['Survived']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,
                                                      random_state=42, stratify=y)
print("train:", X_train.shape, "test:", X_test.shape)

問1 — LightGBM を学習し、gain と split の特徴量重要度を DataFrame にまとめよ

📋 コピーmodel_gain = lgb.LGBMClassifier(
    n_estimators=500, num_leaves=63, learning_rate=0.05,
    min_child_samples=20, random_state=42, verbose=-1, importance_type='gain'
)
# ここを実装: fit → gain_imp, split_imp を取得し imp_df を作る
# imp_df の列: feature, gain, split
# gain で降順ソート

問2 — gain ベースの横棒グラフを描け（上位10特徴量、降順）

📋 コピー# ここを実装
# plt.barh で横棒グラフ
# タイトル: 'Top 10 Feature Importances (LightGBM, gain)'

問3 — SHAP 値を計算し、beeswarm プロットを描け

📋 コピー# ここを実装
# shap.TreeExplainer(model_gain) → explainer
# explainer(X_test) → shap_values
# shap.plots.beeswarm(shap_values, max_display=10)

問4 — テストデータ最初のサンプルの予測を SHAP waterfall で説明せよ

📋 コピーidx = 0
print(f"予測確率: {model_gain.predict_proba(X_test.iloc[[idx]])[:, 1][0]:.4f}")
print(f"実際の値: {y_test.iloc[idx]}")
# ここを実装: shap.plots.waterfall(shap_values[idx])

問5 — 考察

1. gain と split で重要度の順位が大きく違う特徴量はあったか？　あった場合、なぜそうなるか説明せよ
2. SHAP 値から「予測に最も影響する特徴量」は何か？　タイタニックのドメイン知識と一致するか？

💡 ヒント

📋 コピー# gain と split を別々のモデルで取得
model_gain = lgb.LGBMClassifier(..., importance_type='gain')
model_gain.fit(X_train, y_train)
gain_imp = model_gain.feature_importances_

model_split = lgb.LGBMClassifier(..., importance_type='split')
model_split.fit(X_train, y_train)
split_imp = model_split.feature_importances_

imp_df = pd.DataFrame({
    'feature': features,
    'gain': gain_imp,
    'split': split_imp
}).sort_values('gain', ascending=False).reset_index(drop=True)

📋 コピー# 問2: 棒グラフ
top10 = imp_df.head(10).sort_values('gain')  # ascending=True で下から上に
fig, ax = plt.subplots(figsize=(8, 5))
ax.barh(top10['feature'], top10['gain'], color='#22c55e')
ax.set_xlabel('Feature Importance (gain)')
ax.set_title('Top 10 Feature Importances (LightGBM, gain)')
plt.tight_layout(); plt.show()

# 問3: SHAP
explainer = shap.TreeExplainer(model_gain)
shap_values = explainer(X_test)
shap.plots.beeswarm(shap_values, max_display=10)

# 問4: waterfall
shap.plots.waterfall(shap_values[0])

✅ 模範解答

📋 コピーimport numpy as np
import pandas as pd
import lightgbm as lgb
import shap
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import roc_auc_score
import warnings
warnings.filterwarnings('ignore')

# ── データ準備（セットアップと同じ） ──
np.random.seed(42)
n = 5000
data = pd.DataFrame({
    'Pclass':      np.random.choice([1, 2, 3], n, p=[0.24, 0.21, 0.55]),
    'Sex':         np.random.choice(['male', 'female'], n, p=[0.65, 0.35]),
    'Age':         np.random.normal(29.7, 14.5, n).clip(1, 80),
    'SibSp':       np.random.choice(range(9), n, p=[0.68,0.23,0.04,0.02,0.01,0.01,0.005,0.005,0.01]),
    'Parch':       np.random.choice(range(7), n, p=[0.76,0.13,0.09,0.005,0.005,0.005,0.005]),
    'Fare':        np.abs(np.random.normal(32, 50, n)),
    'Embarked':    np.random.choice(['S', 'C', 'Q'], n, p=[0.72, 0.19, 0.09]),
    'Cabin_type':  np.random.choice(['A','B','C','D','E','None'], n,
                                    p=[0.05,0.08,0.12,0.1,0.07,0.58]),
    'FamilySize':  None,
    'IsAlone':     None,
})
data['FamilySize'] = data['SibSp'] + data['Parch'] + 1
data['IsAlone'] = (data['FamilySize'] == 1).astype(int)
data['Survived'] = (
    (data['Sex'] == 'female') * 0.55 + (data['Pclass'] == 1) * 0.25 +
    (data['Age'] < 16) * 0.2 + np.random.normal(0, 0.12, n)
) > 0.4
data['Survived'] = data['Survived'].astype(int)
for col in ['Sex', 'Embarked', 'Cabin_type']:
    data[col] = LabelEncoder().fit_transform(data[col])

features = ['Pclass', 'Sex', 'Age', 'SibSp', 'Parch', 'Fare',
            'Embarked', 'Cabin_type', 'FamilySize', 'IsAlone']
X = data[features]
y = data['Survived']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,
                                                      random_state=42, stratify=y)

# ── 問1: 特徴量重要度 ──
model_gain = lgb.LGBMClassifier(
    n_estimators=500, num_leaves=63, learning_rate=0.05,
    min_child_samples=20, random_state=42, verbose=-1, importance_type='gain'
)
model_gain.fit(X_train, y_train)
gain_imp = model_gain.feature_importances_

model_split = lgb.LGBMClassifier(
    n_estimators=500, num_leaves=63, learning_rate=0.05,
    min_child_samples=20, random_state=42, verbose=-1, importance_type='split'
)
model_split.fit(X_train, y_train)
split_imp = model_split.feature_importances_

imp_df = pd.DataFrame({
    'feature': features,
    'gain':    gain_imp,
    'split':   split_imp,
}).sort_values('gain', ascending=False).reset_index(drop=True)

print("=" * 60)
print(f"{'Feature':<14} {'Gain':>10} {'Split':>8}")
print("-" * 60)
for _, row in imp_df.iterrows():
    print(f"  {row['feature']:<12} {row['gain']:>10.2f} {row['split']:>8.0f}")

auc = roc_auc_score(y_test, model_gain.predict_proba(X_test)[:, 1])
print(f"\nTest AUC: {auc:.4f}")

# ── 問2: 棒グラフ ──
top10 = imp_df.head(10).sort_values('gain')
fig, ax = plt.subplots(figsize=(8, 5))
ax.barh(top10['feature'], top10['gain'], color='#22c55e', alpha=0.85)
ax.set_xlabel('Feature Importance (gain)', fontsize=12)
ax.set_title('Top 10 Feature Importances (LightGBM, gain)', fontsize=13)
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig('feature_importance.png', dpi=120)
plt.show()

# ── 問3: SHAP beeswarm ──
explainer = shap.TreeExplainer(model_gain)
shap_values = explainer(X_test)
shap.plots.beeswarm(shap_values, max_display=10)

# ── 問4: waterfall（1サンプル） ──
idx = 0
prob   = model_gain.predict_proba(X_test.iloc[[idx]])[:, 1][0]
actual = y_test.iloc[idx]
print(f"\nサンプル #{idx}")
print(f"  予測確率: {prob:.4f}（{'生存予測' if prob > 0.5 else '死亡予測'}）")
print(f"  実際の値: {'生存(1)' if actual == 1 else '死亡(0)'}")
shap.plots.waterfall(shap_values[idx])

# ── ボーナス: SHAP ベースの特徴量選択 ──
mean_abs_shap = np.abs(shap_values.values).mean(axis=0)
shap_imp_df = pd.DataFrame({
    'feature': features,
    'mean_abs_shap': mean_abs_shap
}).sort_values('mean_abs_shap', ascending=False)

print("\nSHAP 重要度（mean |SHAP|）")
for _, row in shap_imp_df.iterrows():
    bar = "█" * int(row['mean_abs_shap'] * 100)
    print(f"  {row['feature']:<14}: {row['mean_abs_shap']:.4f}  {bar}")

# 低重要度特徴量（IsAlone, SibSp, Parch など）を削除して再学習
low_imp_features = shap_imp_df[shap_imp_df['mean_abs_shap'] < 0.01]['feature'].tolist()
if low_imp_features:
    top_features = [f for f in features if f not in low_imp_features]
    model_reduced = lgb.LGBMClassifier(n_estimators=500, random_state=42, verbose=-1)
    model_reduced.fit(X_train[top_features], y_train)
    auc_reduced = roc_auc_score(y_test, model_reduced.predict_proba(X_test[top_features])[:, 1])
    print(f"\n全特徴量 AUC:        {auc:.4f}")
    print(f"低重要度除外 AUC:    {auc_reduced:.4f}")
    print(f"除外した特徴量: {low_imp_features}")

🪜 Step-by-Step 解説

1なぜ gain を使うべきか（split の罠）

📋 コピー# Fare は連続値でカーディナリティが高い
# → 決定木が「1000円で分割、1001円で分割、…」と細かく分岐しがち
# → split（分岐回数）が多くなり、重要度が水増しされる

# gain の場合:
# → 各分岐の「実際の不純度改善量」を計算して平均
# → 分岐回数に左右されない信頼性の高い指標

# 結果の違い例（Titanic風データ）
# gain 順位: Sex > Pclass > Age > Fare > FamilySize
# split 順位: Sex > Pclass > Age > Fare > FamilySize  ←大体同じだが
#             Fareの split だけ Ageより高くなることがある

2SHAP の TreeExplainer とは

📋 コピー# TreeExplainer は決定木系に特化した高速なSHAP実装
# 全サンプルのSHAP値を O(TLD²) で計算（T=ツリー数, L=葉, D=深さ）

explainer = shap.TreeExplainer(model_gain)
shap_values = explainer(X_test)  # Explanation オブジェクト

# .values: SHAP 値の行列 (n_samples, n_features)
# .base_values: 基準値（全体の平均予測値）
# .data: 特徴量の実際の値

print(f"SHAP values shape: {shap_values.values.shape}")    # (1000, 10)
print(f"base_value（基準値）: {shap_values.base_values[0]:.4f}")  # ≈ 平均生存率
print(f"SHAP合計 ≈ 予測値: {shap_values.values[0].sum() + shap_values.base_values[0]:.4f}")

3beeswarm プロットの読み方

📋 コピー# beeswarm プロットの各要素
# - Y軸: 特徴量（|SHAP|の平均で降順）
# - X軸: SHAP値（正 = 生存確率を押し上げ、負 = 押し下げ）
# - 色: 特徴量の値（赤 = 高い値、青 = 低い値）
# - 点の密度: 多くのサンプルが同じ SHAP値を持つ = 点が重なる

# 読み方の例
# Sex の赤い点（female=1）が右（+方向）に集中している
# → 女性であることが生存確率を大きく押し上げる

# Pclass の青い点（1等=1）が右に、赤い点（3等=3）が左にある
# → 1等 → 生存UP、3等 → 生存DOWN

4SHAP を使った特徴量選択の手順

📋 コピー# Step A: 全特徴量で学習
# Step B: mean |SHAP| で重要度をランキング
mean_abs_shap = np.abs(shap_values.values).mean(axis=0)
# Step C: 重要度が低い特徴量を候補としてリスト
# Step D: 候補を除外して再学習 → AUC を比較
# Step E: AUC が維持されていれば除外を決定

# Kaggle の実践: IsAlone は FamilySize と情報重複
# → IsAlone を削除しても AUC がほぼ変わらないことが多い
# → シンプルなモデルの方が汎化しやすい

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

Shapley値とは（ゲーム理論から）

不純度（gini）とは

gain = 分岐後の Gini 不純度の改善量。大きいほど「効果的な分岐」を作れる特徴量。

🏆 Kaggleでの実践的な使い方

📋 コピー# Kaggle コンペでよく使うパターン
# 特徴量選択のループ
n_top = 7
top_features = shap_imp_df.head(n_top)['feature'].tolist()
model_top = lgb.LGBMClassifier(n_estimators=500, random_state=42, verbose=-1)
model_top.fit(X_train[top_features], y_train)
auc_top = roc_auc_score(y_test, model_top.predict_proba(X_test[top_features])[:, 1])
print(f"Top {n_top} features AUC: {auc_top:.4f} (全特徴量: {auc:.4f})")

# ランク平均 SHAP（複数モデルのSHAP重要度を平均してより安定した特徴量選択）
# shap1 = shap.TreeExplainer(lgb_model)(X_test).values
# shap2 = shap.TreeExplainer(xgb_model)(X_test).values
# mean_shap = (np.abs(shap1).mean(0) + np.abs(shap2).mean(0)) / 2

🚫 よくある誤解・ミス

🚀 次のステップ

• 発展: SHAP Interaction Values — 2つの特徴量の「掛け算効果」を測る。shap.plots.scatter(shap_values[:, 'Age'], color=shap_values[:, 'Sex']) で視覚化。
• 次回予告（Day 068）: Phase 2 総復習 — House Prices コンペ End-to-End パイプライン（提出まで）

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