【SHAP】スタッキング(stacking)で特徴量の解釈はできるのか|pythonアンサンブル学習

import itertools
import numpy as np
import pandas as pd
import shap
from sklearn.datasets import load_breast_cancer
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import GridSearchCV, train_test_split
from sklearn.ensemble import StackingClassifier
from sklearn.preprocessing import PolynomialFeatures

# データのロード
data = load_breast_cancer()
X, y = data.data, data.target
feature_names = data.feature_names

#シード値の設定
seed=42

#交差項の作成
#２次元まで
poly = PolynomialFeatures(degree=2, interaction_only=True, include_bias=False)
X_poly = poly.fit_transform(X)
feature_names_poly = poly.get_feature_names(input_features=feature_names)

#特徴量選択
#分散分析
selector = SelectKBest(f_classif, k=30)
X_selected = selector.fit_transform(X_poly, y)
selected_feature_indices = selector.get_support(indices=True)
selected_feature_names = [feature_names_poly[i] for i in selected_feature_indices]

#教師データと評価データの分割
X_train, X_test, y_train, y_test = train_test_split(X_selected, y, test_size=0.2, random_state=seed)

#基本モデルとメタモデルの設定
base_models = [
    ('rf', RandomForestClassifier(random_state=seed)),
    ('gb', GradientBoostingClassifier(random_state=seed))
]
meta_model = LogisticRegression(random_state=42)

#ハイパーパラメータのチューニング
#グリッドサーチ
param_grid_rf = {
    'min_samples_leaf': [3, 5, 7],
    'max_depth': [5, 8, 10]
}
param_grid_gb = {
    'max_depth': [3, 5, 8],
    'learning_rate': [0.01, 0.1, 0.2]
}
grid_search_rf = GridSearchCV(base_models[0][1], param_grid_rf, scoring='roc_auc', cv=5)
grid_search_gb = GridSearchCV(base_models[1][1], param_grid_gb, scoring='roc_auc', cv=5)
grid_search_rf.fit(X_train, y_train)
grid_search_gb.fit(X_train, y_train)

#基本モデルごとに特徴量重要度を算出
rf_feature_importances = grid_search_rf.best_estimator_.feature_importances_
gb_feature_importances = grid_search_gb.best_estimator_.feature_importances_

#特徴量重要度を正規化
rf_feature_importances = rf_feature_importances / np.sum(rf_feature_importances)
gb_feature_importances = gb_feature_importances / np.sum(gb_feature_importances)

#教師データから予測を行う
rf_test_predictions = grid_search_rf.best_estimator_.predict_proba(X_test)[:, 1]
gb_test_predictions = grid_search_gb.best_estimator_.predict_proba(X_test)[:, 1]

#各基本モデルの予測と、教師データの相関
rf_corr = np.corrcoef(rf_test_predictions, y_test)[0, 1]
gb_corr = np.corrcoef(gb_test_predictions, y_test)[0, 1]

#基本モデルの重みを計算
rf_weight = rf_corr / (rf_corr + gb_corr)
gb_weight = gb_corr / (rf_corr + gb_corr)

#勾配ブースティング（XGBoost, LightGBM, CatBoostなど）で作成したモデルを読み込み、Shap値を導くためのインスタンス
rf_explainer = shap.TreeExplainer(grid_search_rf.best_estimator_)
gb_explainer = shap.TreeExplainer(grid_search_gb.best_estimator_)

#各基本モデルのShap値を計算
rf_shap_values = np.array(rf_explainer.shap_values(X_test))[1]
gb_shap_values = np.array(gb_explainer.shap_values(X_test))[1]

#重みをつけたShap値を計算
stacking_shap_values = rf_weight * rf_shap_values + gb_weight * gb_shap_values

#絶対値をとり、平均化
stacking_feature_importances_shap = np.mean(np.abs(stacking_shap_values), axis=0)

#特徴量重要度を降順にする
sorted_indices = np.argsort(stacking_feature_importances_shap)[::-1]
sorted_feature_names = [selected_feature_names[i] for i in sorted_indices]
sorted_feature_importances_shap = stacking_feature_importances_shap[sorted_indices]

#視覚化
plt.figure(figsize=(12, 6))
plt.bar(selected_feature_names, sorted_feature_importances_shap)
plt.xticks(rotation=90)
plt.xlabel("Features")
plt.ylabel("SHAP Importance")
plt.title("SHAP-based feature importances for the stacking model")
plt.show()