「日拱一码」030 机器学习常用库——SHAP

## 基于模型类型的方法
# 1. TreeExplainer
import shap
import xgboost
from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import LabelEncoder
import pandas as pd

# 加载 Ames 房价数据集
housing = fetch_openml(name="house_prices", as_frame=True)
X = housing.data
y = housing.target

y = y.astype(float)

# 对分类变量进行标签编码
label_encoders = {}
for column in X.select_dtypes(include=['object']).columns:
    le = LabelEncoder()
    X[column] = le.fit_transform(X[column])
    label_encoders[column] = le

# 分割数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 训练 XGBoost 模型
model = xgboost.XGBRegressor(objective='reg:squarederror', random_state=42)
model.fit(X_train, y_train)

# 使用 SHAP 的 TreeExplainer
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# 绘制特征重要性摘要图
shap.summary_plot(shap_values, X_test, feature_names=X.columns)

KernelExplainer

通过对数据进行扰动来估计SHAP值，适用于任何模型，但计算成本较高，适合小型数据集

# 2. KernelExplainer
import shap
from sklearn.ensemble import RandomForestRegressor
from sklearn.datasets import fetch_openml
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split

# 加载 Ames 房价数据集
housing = fetch_openml(name="house_prices", as_frame=True)
X = housing.data
y = housing.target

# 将目标变量转换为数值类型（如果需要）
y = y.astype(float)

# 对分类变量进行标签编码
label_encoders = {}
for column in X.select_dtypes(include=['object']).columns:
    le = LabelEncoder()
    X[column] = le.fit_transform(X[column])
    label_encoders[column] = le
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 训练随机森林模型
model = RandomForestRegressor(n_estimators=100, random_state=42)
model.fit(X_train, y_train)

# 使用KernelExplainer
explainer = shap.KernelExplainer(model.predict, X_train)
shap_values = explainer.shap_values(X_test, nsamples=100)

# 绘制特征重要性摘要图
shap.summary_plot(shap_values, X_test)

DeepExplainer

通过近似模型的梯度来计算SHAP值，用于深度学习模型（如TensorFlow和Keras模型）

# 3. DeepExplainer
import numpy as np
import tensorflow as tf
from tensorflow.keras.datasets import mnist
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Conv2D, Flatten, MaxPooling2D
import shap

# 加载 MNIST 数据集
(X_train, y_train), (X_test, y_test) = mnist.load_data()

# 数据预处理
X_train = X_train.reshape(-1, 28, 28, 1).astype('float32') / 255.0
X_test = X_test.reshape(-1, 28, 28, 1).astype('float32') / 255.0

model = Sequential([
    Conv2D(32, kernel_size=(3, 3), activation='relu', input_shape=(28, 28, 1)),
    MaxPooling2D(pool_size=(2, 2)),
    Flatten(),
    Dense(128, activation='relu'),
    Dense(10, activation='softmax')
])

model.compile(optimizer='adam', loss='sparse_categorical_crossentropy', metrics=['accuracy'])
model.fit(X_train, y_train, epochs=5, batch_size=32, validation_split=0.2)

# 选择部分训练数据作为背景数据
background = X_train[np.random.choice(X_train.shape[0], 100, replace=False)]

# 初始化 DeepExplainer
explainer = shap.DeepExplainer(model, background)

# 计算 SHAP 值
shap_values = explainer.shap_values(X_test[:100])  # 计算前 100 个测试样本的 SHAP 值

# 可视化前 10 个测试样本的 SHAP 值
shap.image_plot(shap_values, -X_test[:10])

基于解释目标的方法

全局解释

用于解释整个模型的特征重要性

## 基于解释目标的方法
# 1. 全局解释

import shap
import xgboost
from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import LabelEncoder
import pandas as pd

# 加载 Ames 房价数据集
housing = fetch_openml(name="house_prices", as_frame=True)
X = housing.data
y = housing.target

# 将目标变量转换为数值类型（如果需要）
y = y.astype(float)

# 对分类变量进行标签编码
label_encoders = {}
for column in X.select_dtypes(include=['object']).columns:
    le = LabelEncoder()
    X[column] = le.fit_transform(X[column])
    label_encoders[column] = le

# 分割数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 训练 XGBoost 模型
model = xgboost.XGBRegressor(objective='reg:squarederror', random_state=42)
model.fit(X_train, y_train)
explainer = shap.TreeExplainer(model)
# 使用TreeExplainer计算全局SHAP值
shap_values = explainer.shap_values(X_test)
shap.summary_plot(shap_values, X_test, plot_type='bar')  # 绘制条形图

局部解释

用于解释单个预测的特征贡献

# 2. 局部解释

import shap
import xgboost
from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import LabelEncoder
import pandas as pd
import numpy as np

# 加载 Ames 房价数据集
housing = fetch_openml(name="house_prices", as_frame=True)
X = housing.data
y = housing.target

# 将目标变量转换为数值类型（如果需要）
y = y.astype(float)

# 对分类变量进行标签编码
label_encoders = {}
for column in X.select_dtypes(include=['object']).columns:
    le = LabelEncoder()
    X[column] = le.fit_transform(X[column])
    label_encoders[column] = le

# 分割数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 训练 XGBoost 模型
model = xgboost.XGBRegressor(objective='reg:squarederror', random_state=42)
model.fit(X_train, y_train)

# 使用 SHAP 的 TreeExplainer
explainer = shap.TreeExplainer(model)

# 使用 TreeExplainer 计算 SHAP 值
shap_values = explainer.shap_values(X_test)

# 检查 explainer.expected_value 是否为标量
if np.isscalar(explainer.expected_value):
    expected_value = explainer.expected_value
else:
    expected_value = explainer.expected_value[0]  # 如果是数组，取第一个元素

shap.force_plot(expected_value, shap_values[0, :], X_test.iloc[0, :])

基于可视化的方法

Summary Plot

用于展示所有特征的SHAP值分布

shap.summary_plot(shap_values, X_test)

Dependence Plot

用于展示单个特征的SHAP值与特征值之间的关系

# 2. Dependence Plot

shap.dependence_plot('OverallQual', shap_values, X_test)  # 以房价数据集中的OverallQual特征为例

Force Plot

用于展示单个预测的特征贡献

# 3. Force Plot

shap.force_plot(explainer.expected_value, shap_values[0, :], X_test[0, :])

基于特征交互的方法

用于分析特征之间的交互效应

## 基于特征交互的方法
# 使用TreeExplainer计算特征交互值
shap_interaction_values = explainer.shap_interaction_values(X_test)
shap.dependence_plot((0, 1), shap_interaction_values, X_test)  # 绘制特征0和特征1的交互图