SCI图表复现：优化SHAP特征贡献图展示更多模型细节
2025-01-06
背景
机器学习模型的可解释性需求：
随着机器学习模型在各个行业的广泛应用，尤其是在高风险领域（如金融和医疗），模型的可解释性变得越来越重要，简单来说，模型的预测过程需要能够清晰地向用户展示，以便增加信任度，特别是当决策会对人类产生重大影响时
为什么 SHAP 是解释性工具的前沿？
SHAP是一种先进的解释工具，它基于博弈论中的 Shapley 值，能够为机器学习模型的预测提供特征重要性分数，这些分数表示每个特征对模型预测的贡献，无论是正面影响还是负面影响，且能为模型的全局和局部行为提供清晰的解释
参考科学研究并复现经典图表
为了展示 SHAP 在实际应用中的强大功能，本篇文章将复现SCI里的shap特征贡献图表，图中展示了对两种污染物（NO 和 NO₂）进行的 SHAP 特征重要性分析：
此图表通过 SHAP 分析了模型中每个特征的重要性，特别是在全球（图a、b）和局部（图g、h、i）层面的解释，通过这个图表，读者能够清晰地理解每个特征如何影响预测结果，这里我们主要复现图表a、b，a、b图表原理一样
目标与方法
本文将基于类似的方法，使用帕金森症数据集，通过随机森林模型和 SHAP 工具，生成特征贡献度的可视化图表，展示如何使用 SHAP 解释模型预测结果，将帮助您掌握如何复现类似的科学可视化
代码实现
数据读取
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['axes.unicode_minus'] = False
df = pd.read_excel('Parkinsons Telemonitoring.xlsx')
# 划分特征和目标变量
X = df.drop(['total_UPDRS', 'motor_UPDRS'], axis=1)
y = df['total_UPDRS']
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,
random_state=42)
df.head()
加载数据，提取特征和目标变量 total_UPDRS，并将数据集按 80% 训练集和 20% 测试集的比例进行划分，以便后续模型训练和测试，数据集包含来自 42 位帕金森症患者的记录，通过 16 项声音测量特征来监测病情，特征包括抖动（Jitter）、振幅颤动（Shimmer）、噪声与谐波比（NHR）等声音特征，目标变量是帕金森病评估标准（UPDRS，Unified Parkinson’s Disease Rating Scale），分为总评分（total_UPDRS）和运动评分（motor_UPDRS），该数据集主要用于预测患者的 UPDRS 评分，以帮助评估帕金森症的严重程度，如需获取数据集进行复现操作，您可以通过添加作者微信联系获取
模型构建及评价
RF模型训练
from sklearn.ensemble import RandomForestRegressor
# 创建随机森林回归器实例，并设置参数
rf_regressor = RandomForestRegressor(
n_estimators=100, # 'n_estimators'是森林中树的数量。默认是100，可以根据需要调整。
criterion='squared_error', # 'criterion'参数指定用于拆分的质量指标。'squared_error'（默认）表示使用均方误差，另一选项是'absolute_error'。
max_depth=7, # 'max_depth'限制每棵树的最大深度。'None'表示不限制深度。
min_samples_split=2, # 'min_samples_split'指定一个节点分裂所需的最小样本数。默认是2。
min_samples_leaf=1, # 'min_samples_leaf'指定叶子节点所需的最小样本数。默认是1。
min_weight_fraction_leaf=0.0, # 'min_weight_fraction_leaf'与'min_samples_leaf'类似，但基于总样本权重。默认是0.0。
random_state=42, # 'random_state'控制随机数生成，以便结果可复现。42是一个常用的随机种子。
max_leaf_nodes=None, # 'max_leaf_nodes'限制每棵树的最大叶子节点数。'None'表示不限制。
min_impurity_decrease=0.0 # 'min_impurity_decrease'在分裂节点时要求的最小不纯度减少量。默认是0.0。
)
# 训练模型
rf_regressor.fit(X_train, y_train)代码创建并训练一个具有指定参数的随机森林回归模型，以预测 total_UPDRS
模型评价指标可视化
from sklearn import metrics
# 预测
y_pred_train =rf_regressor.predict(X_train)
y_pred_test = rf_regressor.predict(X_test)
y_pred_train_list = y_pred_train.tolist()
y_pred_test_list = y_pred_test.tolist()
# 计算训练集的指标
mse_train = metrics.mean_squared_error(y_train, y_pred_train_list)
rmse_train = np.sqrt(mse_train)
mae_train = metrics.mean_absolute_error(y_train, y_pred_train_list)
r2_train = metrics.r2_score(y_train, y_pred_train_list)
# 计算测试集的指标
mse_test = metrics.mean_squared_error(y_test, y_pred_test_list)
rmse_test = np.sqrt(mse_test)
mae_test = metrics.mean_absolute_error(y_test, y_pred_test_list)
r2_test = metrics.r2_score(y_test, y_pred_test_list)
# 将指标放入列表
metrics_labels = ['MSE', 'RMSE', 'MAE', 'R-squared']
train_metrics = [mse_train, rmse_train, mae_train, r2_train]
test_metrics = [mse_test, rmse_test, mae_test, r2_test]
# 创建柱状图
x = np.arange(len(metrics_labels)) # 横坐标位置
width = 0.35 # 柱子的宽度
fig, ax = plt.subplots()
# 训练集和测试集的柱子
bars1 = ax.bar(x - width/2, train_metrics, width, label='Train')
bars2 = ax.bar(x + width/2, test_metrics, width, label='Test')
# 添加标签和标题
ax.set_ylabel('Scores')
ax.set_title('Comparison of Train and Test Set Metrics')
ax.set_xticks(x)
ax.set_xticklabels(metrics_labels)
ax.legend()
# 在每个柱子上显示数值
def autolabel(bars):
"""在每个柱子上显示数值."""
for bar in bars:
height = bar.get_height()
ax.annotate('{}'.format(round(height, 3)),
xy=(bar.get_x() + bar.get_width() / 2, height),
xytext=(0, 3), # 3 点垂直偏移
textcoords="offset points",
ha='center', va='bottom')
autolabel(bars1)
autolabel(bars2)
fig.tight_layout()
plt.savefig("Comparison of Train and Test Set Metrics.pdf", format='pdf',bbox_inches='tight')
plt.show()
代码通过计算和比较模型在训练集和测试集上的误差和拟合优度指标（MSE、RMSE、MAE、R²），并使用柱状图可视化两者的表现，帮助评估模型的性能是否存在过拟合或欠拟合的情况，从这些指标可以看出，模型在训练集和测试集上的表现较为接近，说明模型没有严重的过拟合或欠拟合现象，虽然测试集上的误差略高于训练集，但差异并不大，表明模型具有较好的泛化能力
模型预测可视化
import seaborn as sns
# 创建一个包含训练集和测试集真实值与预测值的数据框
data_train = pd.DataFrame({
'True': y_train,
'Predicted': y_pred_train,
'Data Set': 'Train'
})
data_test = pd.DataFrame({
'True': y_test,
'Predicted': y_pred_test,
'Data Set': 'Test'
})
data = pd.concat([data_train, data_test])
# 自定义调色板
palette = {'Train': '#b4d4e1', 'Test': '#f4ba8a'}
# 创建 JointGrid 对象
plt.figure(figsize=(8, 6), dpi=1200)
g = sns.JointGrid(data=data, x="True", y="Predicted", hue="Data Set", height=10, palette=palette)
# 绘制中心的散点图
g.plot_joint(sns.scatterplot, alpha=0.5)
# 添加训练集的回归线
sns.regplot(data=data_train, x="True", y="Predicted", scatter=False, ax=g.ax_joint, color='#b4d4e1', label='Train Regression Line')
# 添加测试集的回归线
sns.regplot(data=data_test, x="True", y="Predicted", scatter=False, ax=g.ax_joint, color='#f4ba8a', label='Test Regression Line')
# 添加边缘的柱状图
g.plot_marginals(sns.histplot, kde=False, element='bars', multiple='stack', alpha=0.5)
# 添加拟合优度文本在右下角
ax = g.ax_joint
ax.text(0.95, 0.1, f'Train $R^2$ = {r2_train:.3f}', transform=ax.transAxes, fontsize=12,
verticalalignment='bottom', horizontalalignment='right', bbox=dict(boxstyle="round,pad=0.3", edgecolor="black", facecolor="white"))
ax.text(0.95, 0.05, f'Test $R^2$ = {r2_test:.3f}', transform=ax.transAxes, fontsize=12,
verticalalignment='bottom', horizontalalignment='right', bbox=dict(boxstyle="round,pad=0.3", edgecolor="black", facecolor="white"))
# 在左上角添加模型名称文本
ax.text(0.75, 0.99, 'Model = RF', transform=ax.transAxes, fontsize=12,
verticalalignment='top', horizontalalignment='left', bbox=dict(boxstyle="round,pad=0.3", edgecolor="black", facecolor="white"))
# 添加中心线
ax.plot([data['True'].min(), data['True'].max()], [data['True'].min(), data['True'].max()], c="black", alpha=0.5, linestyle='--', label='x=y')
ax.legend()
plt.savefig("TrueFalse.pdf", format='pdf', bbox_inches='tight')
plt.show()
代码通过散点图、回归线、直方图和拟合优度（R²）值的可视化方式，直观展示模型在训练集和测试集上的预测表现，对角线 x=y 表示理想状态下的预测，散点的偏离程度和回归线的拟合情况则表明了模型的实际预测能力，通过这些图表，可以很好地评估模型的准确性和泛化能力，这部分代码具体参考文章——用图表说话：如何有效呈现回归预测模型结果
shap原始特征贡献可视化
import shap
# 构建 shap解释器
explainer = shap.TreeExplainer(rf_regressor)
# 计算测试集的shap值
shap_values = explainer.shap_values(X_test)
# 特征标签
labels = X_test.columns
# 绘制SHAP值总结图（Summary Plot）
plt.figure(figsize=(15, 5))
shap.summary_plot(shap_values, X_test, plot_type="bar", show=False)
plt.title("SHAP_Feature_Importance_Raw_Output")
plt.savefig("SHAP_Feature_Importance_Raw_Output.pdf", format='pdf',bbox_inches='tight')
plt.show()
这段代码通过 SHAP 库内置的函数计算模型在测试集上每个特征的 SHAP 值，并自动生成一个条形图，总结各个特征对模型预测的重要性，条形图是由 SHAP 库的 shap.summary_plot() 函数生成的，它能够直观地展示哪些特征在模型预测中最为关键，从而提供全局层面的模型可解释性，这意味着用户无需手动绘制图表，直接利用 SHAP 的内置函数即可快速得到结果，但是它并不支持直接生成文献一样的特征贡献图，而是需要我们根据原理去进行图表绘制
数据整理
# 计算每个特征的贡献度
feature_contributions = np.abs(shap_values).mean(axis=0)
# 创建一个DataFrame，其中一列是特征名，另一列是特征贡献度
contribution_df = pd.DataFrame({
'Feature': labels,
'Contribution': feature_contributions
})
# 创建类别规则
categories = ['Basic Info', 'Jitter', 'Shimmer', 'Noise', 'Non-linear']
# 特征对应的类别
category_map = {
'age': 'Basic Info',
'sex': 'Basic Info',
'test_time': 'Basic Info',
'Jitter(%)': 'Jitter',
'Jitter(Abs)': 'Jitter',
'Jitter:RAP': 'Jitter',
'Jitter:PPQ5': 'Jitter',
'Jitter:DDP': 'Jitter',
'Shimmer': 'Shimmer',
'Shimmer(dB)': 'Shimmer',
'Shimmer:APQ3': 'Shimmer',
'Shimmer:APQ5': 'Shimmer',
'Shimmer:APQ11': 'Shimmer',
'Shimmer:DDA': 'Shimmer',
'NHR': 'Noise',
'HNR': 'Noise',
'RPDE': 'Non-linear',
'DFA': 'Non-linear',
'PPE': 'Non-linear'
}
# 将类别映射到DataFrame
contribution_df_sorted['Category'] = contribution_df_sorted['Feature'].map(category_map)
contribution_df_sorted
这里代码是计算每个特征对模型预测的平均贡献度，并根据特征类别对其进行分类，具体解释如下，特征的重要程度计算原理——shap样本值取绝对值的平均值从而得到每个特征的重要程度，然后原始数据对于每个特征是没有具体的类别划分的，这里我们根据特征的实际含义进行划分：将帕金森症数据集的特征按类型划分为五类：基本信息（如年龄、性别等）、抖动特征（Jitter）、振幅抖动特征（Shimmer）、噪声特征（Noise）以及非线性特征（Non-linear）方便后续的文献复现工作
文献复现
环形图绘制
# 按类别和贡献度对数据进行排序，确保同一类别的特征在一起，贡献度从高到低排列
contribution_df_sorted = contribution_df_sorted.sort_values(by=['Category', 'Contribution'], ascending=[True, False])
# 创建一个用于生成颜色渐变的函数
def get_color_gradient(base_color, num_shades):
# 生成从浅到深的颜色渐变
gradient = np.linspace(0.4, 1, num_shades) # 生成从较浅（0.4）到原色（1）的渐变
return [(base_color[0], base_color[1], base_color[2], shade) for shade in gradient]
# 为五个类别定义颜色
category_colors = {
'Basic Info': (0.9, 0.7, 0.2, 1), # 黄色
'Jitter': (0.6, 0.3, 0.9, 1), # 紫色
'Noise': (0.7, 0.3, 0.3, 1), # 暗红
'Non-linear': (0.2, 0.9, 0.9, 1), # 青色
'Shimmer': (0.3, 0.6, 0.9, 1), # 浅蓝
}
# 默认颜色，如果类别未定义时使用
default_color = (0.8, 0.8, 0.8, 1) # 灰色
# 获取内圈和外圈的贡献度数据
inner_contribution = contribution_df_sorted.groupby('Category')['Contribution'].sum()
outer_contribution = contribution_df_sorted.set_index('Feature')['Contribution']
# 检查是否有未定义的类别
undefined_categories = set(inner_contribution.index) - set(category_colors.keys())
if undefined_categories:
print(f"Warning: 以下类别没有定义颜色，将使用默认颜色: {undefined_categories}")
# 为每个类别在外圈创建颜色渐变
outer_colors = []
for category in inner_contribution.index:
# 选取当前类别的数据
category_df = contribution_df_sorted[contribution_df_sorted['Category'] == category]
# 获取类别的基础颜色，如果没有定义则使用默认颜色
base_color = category_colors.get(category, default_color)
# 为当前类别生成颜色渐变
gradient_colors = get_color_gradient(base_color, len(category_df))
outer_colors.extend(gradient_colors)
# 内外圈的标签准备
inner_labels = inner_contribution.index
outer_labels = outer_contribution.index
# 绘制同心饼图
fig, ax = plt.subplots(figsize=(8, 8), dpi=1200)
# 绘制内圈饼图（类别级别的饼图），显示百分比，不显示标签
ax.pie(inner_contribution, labels=['']*len(inner_contribution), autopct='%1.1f%%', radius=1,
colors=[category_colors.get(cat, default_color) for cat in inner_labels], wedgeprops=dict(width=0.3, edgecolor='w'))
# 绘制外圈饼图（特征级别的饼图），不显示标签和百分比
ax.pie(outer_contribution, labels=['']*len(outer_labels), radius=0.7, colors=outer_colors, wedgeprops=dict(width=0.3, edgecolor='w'))
# 添加白色中心圆，形成环形图
plt.gca().add_artist(plt.Circle((0, 0), 0.4, color='white'))
# 添加标题
plt.title('Feature and Category Contribution by SHAP')
plt.savefig("Feature and Category Contribution by SHAP.pdf", format='pdf',bbox_inches='tight')
# 显示图表
plt.show()
首先按特征类别和贡献度对数据进行排序，然后根据每个类别和特征的贡献度生成同心饼图，其中外圈展示各类别的总贡献度，内圈展示各个特征的具体贡献度，并通过颜色渐变区分类别内特征的重要性，最终生成一个可视化特征贡献的环形图并保存为 PDF 文件，这里作者关闭标签显示是为了一步一步演示，读者可以显示标签使得可视化更方便阅读
条形图绘制
# 按贡献度从高到低排序
contribution_df_sorted = contribution_df_sorted.sort_values(by='Contribution', ascending=False)
# 准备颜色列表
bar_colors = [category_colors.get(cat, (0.8, 0.8, 0.8, 1)) for cat in contribution_df_sorted['Category']]
# 绘制水平柱状图
fig, ax = plt.subplots(figsize=(10, 8), dpi=1200)
# 绘制条形图
ax.barh(contribution_df_sorted['Feature'], contribution_df_sorted['Contribution'], color=bar_colors)
# 添加图例
handles = [plt.Rectangle((0, 0), 1, 1, color=category_colors[cat]) for cat in category_colors]
labels = list(category_colors.keys())
ax.legend(handles, labels, loc='lower right')
# 设置标签和标题
ax.set_xlabel('Contribution')
ax.set_ylabel('Feature')
ax.set_title('Feature Contributions by Category')
# 反转y轴，以便贡献度最大的特征在顶部
ax.invert_yaxis()
plt.savefig("Feature Contributions by Category.pdf", format='pdf',bbox_inches='tight')
# 显示图表
plt.show()
生成一个水平条形图，按贡献度从高到低显示各个特征的贡献，同时用不同颜色区分特征所属的类别，并添加图例，该图表直观地展示了哪些特征对模型的影响最大
环形图和条形图组合绘图
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
# 按类别和贡献度对数据进行排序，确保同一类别的特征在一起，贡献度从高到低排列
contribution_df_sorted = contribution_df_sorted.sort_values(by=['Category', 'Contribution'], ascending=[True, False])
# 创建一个用于生成颜色渐变的函数
def get_color_gradient(base_color, num_shades):
# 生成从浅到深的颜色渐变
gradient = np.linspace(0.4, 1, num_shades) # 生成从较浅（0.4）到原色（1）的渐变
return [(base_color[0], base_color[1], base_color[2], shade) for shade in gradient]
# 为五个类别定义颜色
category_colors = {
'Basic Info': (0.9, 0.7, 0.2, 1), # 黄色
'Jitter': (0.6, 0.3, 0.9, 1), # 紫色
'Noise': (0.7, 0.3, 0.3, 1), # 暗红
'Non-linear': (0.2, 0.9, 0.9, 1), # 青色
'Shimmer': (0.3, 0.6, 0.9, 1), # 浅蓝
}
# 默认颜色，如果类别未定义时使用
default_color = (0.8, 0.8, 0.8, 1) # 灰色
# 获取内圈和外圈的贡献度数据
inner_contribution = contribution_df_sorted.groupby('Category')['Contribution'].sum()
outer_contribution = contribution_df_sorted.set_index('Feature')['Contribution']
# 检查是否有未定义的类别
undefined_categories = set(inner_contribution.index) - set(category_colors.keys())
if undefined_categories:
print(f"Warning: 以下类别没有定义颜色，将使用默认颜色: {undefined_categories}")
# 为每个类别在外圈创建颜色渐变
outer_colors = []
for category in inner_contribution.index:
# 选取当前类别的数据
category_df = contribution_df_sorted[contribution_df_sorted['Category'] == category]
# 获取类别的基础颜色，如果没有定义则使用默认颜色
base_color = category_colors.get(category, default_color)
# 为当前类别生成颜色渐变
gradient_colors = get_color_gradient(base_color, len(category_df))
outer_colors.extend(gradient_colors)
# 内外圈的标签准备
inner_labels = inner_contribution.index
outer_labels = outer_contribution.index
# 创建图形和子图
fig, ax = plt.subplots(figsize=(10, 8), dpi=1200)
# 设置背景颜色为淡灰色
ax.set_facecolor('#f0f0f0')
# 添加网格线，设置网格线样式
ax.grid(True, which='both', linestyle='--', linewidth=0.7, color='gray', alpha=0.7)
# ---- 绘制柱状图 ----
# 按贡献度从高到低排序
contribution_df_sorted = contribution_df_sorted.sort_values(by='Contribution', ascending=False)
# 准备颜色列表
bar_colors = [category_colors.get(cat, (0.8, 0.8, 0.8, 1)) for cat in contribution_df_sorted['Category']]
# 绘制条形图
ax.barh(contribution_df_sorted['Feature'], contribution_df_sorted['Contribution'], color=bar_colors)
# 添加图例
handles = [plt.Rectangle((0, 0), 1, 1, color=category_colors[cat]) for cat in category_colors]
labels = list(category_colors.keys())
ax.legend(handles, labels, loc='lower right')
# 设置标签和标题
ax.set_xlabel('Contribution')
ax.set_ylabel('Feature')
ax.set_title('Feature Contributions by Category')
# 反转y轴，以便贡献度最大的特征在顶部
ax.invert_yaxis()
# ---- 在柱状图中嵌入同心饼图 ----
# 创建嵌入的轴
inset_ax = inset_axes(ax, width=2, height=2, loc='upper right', bbox_to_anchor=(0.8, 0.35, 0.2, 0.2), bbox_transform=ax.transAxes)
# 绘制内圈饼图（类别级别的饼图），显示百分比，不显示标签
inset_ax.pie(inner_contribution, labels=['']*len(inner_contribution), autopct='%1.1f%%', radius=1,
colors=[category_colors.get(cat, default_color) for cat in inner_labels], wedgeprops=dict(width=0.3, edgecolor='w'))
# 绘制外圈饼图（特征级别的饼图），不显示标签和百分比
inset_ax.pie(outer_contribution, labels=['']*len(outer_labels), radius=0.7, colors=outer_colors, wedgeprops=dict(width=0.3, edgecolor='w'))
# 添加白色中心圆，形成环形图
inset_ax.add_artist(plt.Circle((0, 0), 0.4, color='white'))
plt.savefig("Combined_Feature_Contributions_and_Circular_Chart.pdf", format='pdf',bbox_inches='tight')
# 显示图表
plt.show()
通过组合柱状图和同心饼图的方式，直观展示各个特征在模型中的贡献度，以及各个类别对模型的整体贡献。柱状图显示特征的详细贡献，饼图则提供了类别层面的总体贡献展示，使得用户能够从全局和细节两方面理解特征的重要性，这里读者其实是不需要单独去绘制环形图和条形图的作者只是为了让读者更方便理解，其次这里的环形图和文献的环形图刚好是相反的，作者这里外圈是特征类别总贡献度，内圈是各个特征具体的贡献度映射，对于模型解读并没有区别
可视化解读：年龄（age）作为“基本信息”类别的特征贡献度最高，DFA作为“非线性”特征也具有显著贡献，外圈饼图显示“基本信息”类别占总贡献的 71.5%，而“非线性”特征占 20%，其余类别如“噪声”、“抖动”和“振幅”特征的贡献较小，内圈的特征贡献进一步细分了各类别中具体特征的贡献度，帮助直观理解特征对模型预测的重要性。需要注意的是，这个结果基于演示用的复现数据，对于实际生活中的情况并不一定具有直接的参考价值，如果读者在理解代码时遇到问题，可以联系作者并通过 ChatGPT 协助，结合 AI 更好地理解文章内容
同话题下的热门内容
