????在進行數據分析和機器學習時經常用到shap，本文對shap相關的操作進行演示。波士頓數據集鏈接在這里。

SHAP Analysis Guide

Set up

導入必要包

import pandas as pd
import numpy as np
import lightgbm as lgb
import matplotlib
import matplotlib.pyplot as plt
import shap
import seaborn as sns
import warningsplt.style.use("style.mplstyle")

Load Data

數據集原始論文
關于特征B

import pandas as pd# 讀取 txt 文件
file_path = "shap/boston.txt"  # 替換為你的文件路徑
with open(file_path, "r") as file:lines = file.readlines()# 初始化空列表存儲數據
data = []# 按兩行一組處理數據
for i in range(0, len(lines), 2):# 第一行：前 11 個數據點row1 = list(map(float, lines[i].strip().split()))# 第二行：后 3 個數據點row2 = list(map(float, lines[i + 1].strip().split()))# 將兩行合并為一組完整數據（共 14 個數據點）data.append(row1 + row2)# 定義列名
columns = ["crime rate","% residential zone","% industrial zone","Charles River","NOX concentration","number of rooms","% built before 1940","remoteness","connectedness","tax rate","pupil-teacher ratio","B","% working class",'target'
]# 轉換為 DataFrame
df = pd.DataFrame(data, columns=columns)# 查看數據
print(df.head())

可視化部分數據

刪除特征B

# 刪除 "B" 列
df = df.drop("B", axis=1)

練習修改數據

# 修改數據（例如將 "% working class" 列乘以 2）
df["% working class"] = df["% working class"] * 2
df['target'] = df['target'] * 1000
X = df.drop("target", axis=1).copy()
# 查看處理后的數據
print(df.head())

打印房間預測值

保存經過處理后的數據集

df.to_csv("shap/data.csv", index=False)

打印據相關的數值信息

修改字體(可選)

# 設置字體為系統支持的字體
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = ['Arial', 'DejaVu Sans', 'Liberation Sans', 'Helvetica']

數據分布可視化

fig, axs = plt.subplots(ncols=5, nrows=3, figsize=(14, 9))
gs = axs[1, 2].get_gridspec()
for ax in axs[:, -1]:ax.remove()
axbig = fig.add_subplot(gs[:, -1])
axs = axs.flatten()for col, i in zip(X.columns, [0, 1, 2, 3, 5, 6, 7, 8, 10, 11, 12, 13]):axs[i].hist(X[col])axs[i].set_xlabel(col, size=16)axs[i].grid()axs[i].set_ylim(0, 490)d = df.copy()
d["house price"] = df["target"]
sns.boxplot(y=d["house price"], color="#d45087")
axbig.set_yticklabels(axbig.get_yticks(), rotation=45)
axbig.get_yaxis().set_major_formatter(matplotlib.ticker.FuncFormatter(lambda x, p: format(int(x), ","))
)
axbig.set_xlabel("house price ($)")
axbig.set_ylabel("")
axbig.grid(axis="y")fig.subplots_adjust(hspace=0.5)
fig.text(-0.01, 0.5, "count of instances", va="center", rotation="vertical", size=20)
fig.tight_layout()
fig.show()

Train Model

訓練一個LGBM模型

m = lgb.LGBMRegressor()
m.fit(X, y)

預測并評估模型

Compute SHAP Values

explainer = shap.Explainer(m)
shap_values = explainer(X)

打印數據維度，可以看出shap值的維度與原始數據的維度相同

在 SHAP 的回歸場景下，base value 等于模型對訓練集預測結果的期望

繪制shap分析圖

下圖是第一種圖，蜂巢圖也是summary圖，信息最多。紅色表示數值大，藍色表示數值小，點的密集程度表示數據的分布。特征的重要性從上到下排序，橫軸左側表示負相關，右側表示正相關。

計算分位點

i_med = np.argsort(y_pred)[len(y_pred)//2]
i_max = np.argmax(y_pred)
i_80 = np.argsort(y_pred)[int(len(y_pred)*0.8)]
i_60 = np.argsort(y_pred)[int(len(y_pred)*0.6)]
i_40 = np.argsort(y_pred)[int(len(y_pred)*0.4)]
i_20 = np.argsort(y_pred)[int(len(y_pred)*0.2)]
i_min = np.argmin(y_pred)

下面是第二種圖，瀑布圖，可以展示不同特征對最終預測結果具體的影響。

下圖展示了力圖，紅色越長表示正向作用力越強，模型預測越偏大，反之藍色會讓模型預測偏小。


下面是條狀圖，橫軸表示shap值絕對值的平均值，展示了每一種特征對最終結果的影響程度。

值得注意的是，shap相關的api也是可以進行更改的，用戶可以按照自己的喜好更改api進行繪圖。如下所示，可以人為的繪制出條狀圖，并可以將蜂巢圖的輸出轉化成絕對值的形式。

plt.subplot(2, 1, 1)
# plt.gcf()
shap.plots.bar(shap_values.abs.max(0), max_display=99, show=False)
plt.subplot(2, 1, 2)
shap.plots.beeswarm(shap_values.abs, color="shap_red", max_display=99, show=False, plot_size=None
)
ax = plt.gca()
masv = {}
for feature in ax.get_yticklabels():name = feature.get_text()col_ind = X.columns.get_loc(name)mean_abs_sv = np.mean(np.abs(shap_values.values[:, col_ind]))masv[name] = mean_abs_sv
ax.scatter(masv.values(),[i for i in range(len(X.columns))],zorder=99,label="Mean Absolute SHAP Value",c="k",marker="|",linewidths=3,s=100,
)
ax.legend(frameon=True)
plt.tight_layout()
plt.show()

下圖展示的是依賴圖

n = 5
fig, ax = plt.subplots(1, n, figsize=(15, 3))for i, (k, v) in enumerate(sorted(masv.items(), key=lambda x: x[1], reverse=True)):if i < n:shap.plots.scatter(shap_values[:, k], ax=ax[i], show=False, alpha=0.6)ax[i].grid(axis="y")if i != 0:ax[i].set_ylabel("")ax[i].spines["left"].set_visible(False)ax[i].set_ylim(ax[0].get_ylim())ax[i].set_yticklabels(["" for _ in range(len(ax[0].get_yticks()))])else:ax[i].set_ylabel("SHAP value")
plt.show()

以前5個特征為例，可以看出每個特征和shap值之間的關系，圖一圖二可以看出比較明顯的線性關系。圖四表示距離較近時也有一定的線性關系，但是隨著距離增加以后就線性無關了。

下面展示最后一種散點圖，會自動給出對當前特征(工薪階級)相互作用強的特征作為顏色依據。可以明顯的看出，工薪階級百分比與自身 SHAP 值呈顯著負相關，并且工薪階級占比越大，NOX的排放濃度也越大，房價也越低；反之，工薪階級占比少，NOX濃度高，房價也高。

fig, ax = plt.subplots()
shap.plots.scatter(shap_values[:, "% working class"], color=shap_values, ax=ax)
plt.show()