此圖表通過 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、R2),并使用柱狀圖可視化兩者的表現,幫助評估模型的性能是否存在過擬合或欠擬合的情況,從這些指標可以看出,模型在訓練集和測試集上的表現較為接近,說明模型沒有嚴重的過擬合或欠擬合現象,雖然測試集上的誤差略高于訓練集,但差異并不大,表明模型具有較好的泛化能力
模型預測可視化
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()
代碼通過散點圖、回歸線、直方圖和擬合優度(R2)值的可視化方式,直觀展示模型在訓練集和測試集上的預測表現,對角線 x=y 表示理想狀態下的預測,散點的偏離程度和回歸線的擬合情況則表明了模型的實際預測能力,通過這些圖表,可以很好地評估模型的準確性和泛化能力,這部分代碼具體參考文章——用圖表說話:如何有效呈現回歸預測模型結果
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 更好地理解文章內容