在一篇令人振奋的 Nature Medicine 论文中，作者附上了用于生成 SHAP 图的制图代码，这确实值得庆幸。不过由于不同的环境等问题，我们照搬文章的代码不一定会得到文章中显示的图片。这确实会让我们很愚蒙，甚至一度怀疑自己。
引言
SHAP（SHapley Additive exPlanations）值是一种强大的工具，用于解释机器学习模型的预测结果。在复杂的模型中，理解每个特征的贡献变得尤为重要。然而，复现他人的研究结果往往充满挑战，尤其是在环境配置和数据处理上。
在本文中，我们将通过一篇 Nature Medicine 论文中的 SHAP 图为切入点，逐步复现作者的成果，并探讨在这一过程中可能遇到的各种问题及其解决方案。
环境准备
在开始之前，确保您的环境中安装了以下 Python 库：
pip install shap pandas numpy matplotlib sklearn
这里给出我所使用库的版本，以供大家参考：
- shap (0.46.0)
- pandas (2.1.4)
- numpy (1.26.4)
- matplotlib (3.8.0)
- sklearn(scikit-learn) (1.2.2)
请注意，您可能需要根据论文中的描述准备特定的数据集。不过这篇Nature已经给出了需要的数据集，直接下载使用即可。
GitHub - zeineb-bouzid/sharing-github-nature-medicine
论文中的 SHAP 图制图代码
作者提供的代码段如下所示（其余辅助代码可以看github的共享仓库），可能看似简单，但理解基本参数的意义至关重要。
def plot_Shapley(clf, df_test, cols, set_name, nf): #nf = number of features#Load JS visualization code to notebookshap.initjs()#Shapley explainabilityplt.figure()df = pd.DataFrame(data=df_test, columns=(cols))shap_values_list = []for calibrated_classifier in clf.calibrated_classifiers_:explainer = shap.TreeExplainer(calibrated_classifier.base_estimator)shap_values = explainer.shap_values(df)shap_values_list.append(shap_values)shap_values = np.array(shap_values_list).sum(axis=0) / len(shap_values_list)shap.summary_plot(shap_values[1], features=df, feature_names=np.array(cols), sort=True, show=False, max_display=nf, color_bar=False, plot_type='dot', plot_size=(5.00, 10.00))plt.title('Shapley values of the top '+str(nf)+'/'+str(len(cols))+' features for OMI prediction (derivation cohort)', fontweight="bold")plt.show()
我们将逐段分析这些代码，并讨论可能影响结果的环境因素。
常见错误及解决方案
在复现过程中，您可能会遇到形状不匹配或数据类型不一致等问题。我们将列出这些错误的常见原因，并提供有效的解决方案。 首先，当我们复制作者的源码并执行 Python 代码时，可能会出现如下错误，未能复现出预期的图像。仔细阅读错误信息后，我们发现当前版本中没有 base_estimator 属性。这是因为 calibrated_classifiers_ 的版本更新，导致该属性被删除。
为了解决这个问题，我们将 base_estimator 替换为新版本中的 estimator 属性。再次执行代码后，却出现了不同的错误，如下图所示。 可以看到，红框中显示的内容表明，shap.summary_plot 中期望的 shape 参数不正确。
为了定位 shape 的错误，我们在原始代码中增加了相关打印信息（见绿色框），以检查当前输出的 shape 信息。第一个打印信息用于查看数据的具体 shape，而第二个打印信息则检查 SHAP 库输出结果的 shape（即需要可视化的数据的 shape）。
执行代码后，输出结果如右侧红框所示。可以看到，shap_values 是一个三维数组，而 shap.summary_plot 期望的输入数据格式应为 [样本数, 特征数]。因此，我们需要对代码进行适当修改。
逐步复现 SHAP 图
以下是完整的代码示例，您可以按照步骤生成与论文一致的 SHAP 图。
def plot_Shapley(clf, df_test, cols, set_name, nf): # nf = number of features# Load JS visualization code to notebookshap.initjs()# Shapley explainabilityplt.figure()df = pd.DataFrame(data=df_test, columns=cols)# 打印 X_test 的形状print("X_test shape:", df.shape) # 确保它是 (806, 74)shap_values_list = []for calibrated_classifier in clf.calibrated_classifiers_:explainer = shap.TreeExplainer(calibrated_classifier.estimator)shap_values = explainer.shap_values(df)# 打印每个分类器的 SHAP 值的样本数print("shap_values", shap_values.shape)# 只保留正类的 SHAP 值shap_values_list.append(shap_values[:, :, 1]) # 选择正类的 SHAP 值# 将所有 SHAP 值进行平均shap_values = np.mean(shap_values_list, axis=0)# 确保 SHAP 值的形状和 DataFrame 的形状一致print("Final SHAP values shape:", shap_values.shape) # 应该是 (样本数, 特征数)print("Feature matrix shape:", df.shape) # 应该是 (样本数, 特征数)# 检查 SHAP 值的样本数是否与特征矩阵一致if shap_values.shape[0] != df.shape[0]:raise ValueError("SHAP values shape does not match feature matrix shape.")# 绘制 SHAP 值shap.summary_plot(shap_values,features=df,feature_names=cols,sort=True,show=False,max_display=nf,color_bar=False,plot_type='dot',plot_size=(6.00, 12.00))plt.title('Shapley values of the top ' + str(nf) + '/' + str(len(cols)) +' features for OMI prediction (derivation cohort)', fontweight="bold")plt.savefig(set_name + '.pdf', bbox_inches='tight')plt.show()
结论
通过本文的指导，您应该能够成功复现 Nature Medicine 论文中的 SHAP 图。同时，对 SHAP 值的理解将增强您在机器学习领域的可解释性能力。
参考文献
[nature medicine] Machine learning for ECG diagnosis and risk stratification of occlusion myocardial infarction.
希望这篇文章能为您在数据分析和可视化的旅程中提供帮助。如果您有任何问题或想要进一步探讨的主题，请在评论区留言。感谢您的阅读，我们下期再见！
欢迎关注公众号“菜鸟的人工智能之路”，持续更新医学数据处理相关内容。
