背景
随着机器学习在环境和社会经济研究中的应用越来越广泛，人们对模型的透明性和可解释性的需求也在不断增加。本文将探讨如何使用SHAP值在GeoJSON格式的地理数据上可视化和解释机器学习模型。本文所使用的数据来源于德拉瓦谷区域规划委员会（DVRPC）提供的气候与经济公平筛查工具v1.0，该数据集为不同区域的环境和经济因素提供了关键的洞察
数据与方法
本次分析所用的数据可从DVRPC的数据集中获取，数据格式为GeoJSON，GeoJSON是一种常见的地理数据结构编码格式，可以轻松与Python中的GeoPandas等地理空间库集成
使用该数据集预测空气污染水平，具体为PM2.5浓度（pm25f_pfs），并基于多个环境和社会经济特征作为预测变量，选择的预测变量（自变量）包括：
- ebf_pfs: 暴露于洪水的人口百分比
- hsef: 住房压力因素
- wf_pfs: 劳动力参与率
- n_clt: 气候风险因素
- n_trn: 交通通达性
- n_hsg: 住房质量
- n_pln: 规划指标
- sn_c: 社会需求指数
在本次分析中，采用XGBoost机器学习模型，这是一种因其鲁棒性和灵活性而被广泛应用于回归任务的算法，为了增强模型的可解释性，我们使用了SHAP值，这是一种常用于解释机器学习模型输出的方法，当然这里作者不会过多去考虑模型精确度重点在于SHAP值地图可视化解释ML模型
代码实现
数据读取
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['axes.unicode_minus'] = False
# 读取 GeoJSON 文件
gdf = gpd.read_file("Climate_and_Economic_Justice_Screening_Tool_v1.0.geojson")
# 打印前几行，查看字段信息
gdf.head()
读取一个GeoJSON格式的地理空间数据文件，显示其前几行数据以查看字段信息
数据集分割
from sklearn.model_selection import train_test_split
# 选择因变量和自变量
y = gdf['pm25f_pfs']
X = gdf[['ebf_pfs', 'hsef', 'wf_pfs', 'n_clt', 'n_trn', 'n_hsg', 'n_pln', 'sn_c']]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,
random_state=42)
提取自变量（X）和因变量（y），然后使用8-2比例划分数据集，确保随机种子设置以便重现性
模型构建
import xgboost as xgb
from sklearn.model_selection import GridSearchCV
# XGBoost回归模型参数
params_xgb = {
'learning_rate': 0.02, # 学习率，控制每一步的步长，用于防止过拟合。典型值范围：0.01 - 0.1
'booster': 'gbtree', # 提升方法，这里使用梯度提升树（Gradient Boosting Tree）
'objective': 'reg:squarederror', # 损失函数，这里使用平方误差，适用于回归任务
'max_leaves': 127, # 每棵树的叶子节点数量，控制模型复杂度
'verbosity': 1, # 控制 XGBoost 输出信息的详细程度，0表示无输出，1表示输出进度信息
'seed': 42, # 随机种子，用于重现模型的结果
'nthread': -1, # 并行运算的线程数量，-1表示使用所有可用的CPU核心
'colsample_bytree': 0.6, # 每棵树随机选择的特征比例，用于增加模型的泛化能力
'subsample': 0.7, # 每次迭代时随机选择的样本比例，用于增加模型的泛化能力
'eval_metric': 'rmse' # 评价指标，这里使用均方根误差（rmse）
}
# 初始化XGBoost回归模型
model_xgb = xgb.XGBRegressor(**params_xgb)
# 定义参数网格，用于网格搜索
param_grid = {
'n_estimators': [100, 200, 300, 400, 500], # 树的数量
'max_depth': [3, 4, 5, 6, 7], # 树的深度
'learning_rate': [0.01, 0.02, 0.05, 0.1], # 学习率
}
# 使用GridSearchCV进行网格搜索和k折交叉验证
grid_search = GridSearchCV(
estimator=model_xgb,
param_grid=param_grid,
scoring='neg_mean_squared_error', # 评价指标为负均方误差
cv=5, # 5折交叉验证
n_jobs=-1, # 并行计算
verbose=1 # 输出详细进度信息
)
# 训练模型
grid_search.fit(X_train, y_train)
# 输出最优参数
print("Best parameters found: ", grid_search.best_params_)
print("Best RMSE score: ", (-grid_search.best_score_) ** 0.5) # 还原RMSE
# 使用最优参数训练模型
best_model_xgboost = grid_search.best_estimator_
配置XGBoost回归器的超参数，并使用网格搜索和交叉验证来识别最佳参数。使用最佳参数在训练数据上重新训练模型
SHAP值可视化
点图（Dot Plot）
import shap
# 构建 shap解释器
explainer = shap.TreeExplainer(best_model_xgboost)
# 计算shap值
shap_values = explainer.shap_values(X)
# 特征标签
labels = X.columns
plt.figure(dpi=1200)
shap.summary_plot(shap_values, X, feature_names=labels, plot_type="dot", show=False)
plt.savefig("SHAP_1.pdf", format='pdf', bbox_inches='tight')
展示每个特征的SHAP值分布，清晰显示了每个特征对预测PM2.5水平的正向或负向影响及其幅度
- 横轴： 表示每个特征对模型输出的影响大小和方向， 负值意味着特征对模型输出的负面影响，正值则意味着正面影响
- 纵轴 ： 列出了模型中的各个特征，如wf_pfs（工作参与率）、ebf_pfs（暴露于洪水的比例）等
- 数据点的颜色 ： 代表特征的原始值大小， 红色表示特征值高，蓝色表示特征值低， 颜色变化展示了特征值的高低如何影响SHAP值
- 数据点的分布 ： 展示了特征值的不同范围如何影响模型输出， 数据点越集中，表明特征对输出的影响越大
条形图（Bar Plot）
plt.figure(figsize=(15, 5), dpi=1500)
shap.summary_plot(shap_values, X, plot_type="bar", show=False)
plt.savefig("SHAP_2.pdf", format='pdf', bbox_inches='tight')
plt.show()
根据SHAP值的平均绝对值对特征进行排序，突出显示了模型中最具影响力的预测变量， 其中wf_pfs（工作参与率）是对模型输出影响最大的特征，其次是ebf_pfs（暴露于洪水的比例）和hsef（住房压力因素）
使用GeoJSON绘制SHAP值地图
数据整理
# 将 SHAP 值转换为 DataFrame，并为列名添加 "shap_" 前缀，同时保留原始特征名
shap_df = pd.DataFrame(shap_values, columns=[f'shap_{col}' for col in X.columns])
# 将 SHAP 值和原始特征合并
gdf_with_shap = pd.concat([gdf, shap_df], axis=1)
gdf_with_shap
本次分析的一个重要方面是能够在地理地图上可视化SHAP值，这为区域内模型预测差异的洞察提供了帮助，通过将SHAP值与原始GeoJSON数据结合，创建一个包含详细空间可视化SHAP值的地理空间数据集（gdf_with_shap）
ebf_pfs特征SHAP地图
# 绘制基于 'shap_ebf_pfs' 字段的地图
gdf_with_shap.plot(column='shap_ebf_pfs',
legend=True,
cmap=shap.plots.colors.red_white_blue, # 使用 SHAP 配色
figsize=(10, 10))
# 设置标题
plt.title("SHAP Value for ebf_pfs")
plt.savefig("SHAP_3.pdf", format='pdf', dpi=1200, bbox_inches='tight')
plt.show()
这里先单独绘制一个特征的SHAP值地图，以了解特征如何在不同地理区域影响PM2.5预测水平，ebf_pfs（暴露于洪水的人口百分比）的SHAP地图显示了这一因素在特定区域增加或减少预测PM2.5水平的影响
- 颜色编码： 图中颜色从蓝色到红色渐变，表示SHAP值的范围， 红色表示正的SHAP值，意味着ebf_pfs特征在这些区域对预测值有正向的贡献； 蓝色表示负的SHAP值，意味着ebf_pfs特征对预测值有负向的贡献
- 颜色深浅： 颜色越深，表示SHAP值的绝对值越大，表明该特征对模型输出的影响越显著， 颜色较浅则表示 影响较小
- 空间分布： 可以看到，不同区域有不同的颜色，这反映了ebf_pfs特征在各个地理区域的影响差异， 例如，图中红色深的区域表示ebf_pfs的高值对模型输出的正面影响更大，而蓝色深的区域则表示相反
总结而言，这张图展示了特征ebf_pfs在不同地理区域对模型预测结果的影响，为环境决策提供了更具针对性的区域洞察。需要注意的是，这里的重点在于如何绘制SHAP地图，因此并未深入研究模型的精确度是否足够应用于实际场景
全面的SHAP地图可视化
# 提取所有以 'shap_' 开头的列名
shap_columns = [col for col in gdf_with_shap.columns if col.startswith('shap_')]
# 创建一个画布，设置子图的网格布局
num_cols = len(shap_columns)
fig, axes = plt.subplots(nrows=(num_cols + 2) // 2, ncols=2, figsize=(18, 6 * ((num_cols + 2) // 2)))
# 展平 axes 数组以便于迭代
axes = axes.ravel()
# 遍历每个 SHAP 列，绘制地图
for i, shap_col in enumerate(shap_columns):
gdf_with_shap.plot(column=shap_col,
legend=True,
cmap=shap.plots.colors.red_white_blue,
ax=axes[i]) # 在对应子图上绘制
axes[i].set_title(f"SHAP Value for {shap_col}", fontsize=12)
# 删除多余的子图（如果有的话）
for j in range(i + 1, len(axes)):
fig.delaxes(axes[j])
plt.tight_layout()
plt.savefig("SHAP_4.pdf", format='pdf', dpi=1200, bbox_inches='tight')
plt.show()
为了提供全面的视角，使用网格布局绘制了所有特征的SHAP地图，为不同的社会经济和环境因素如何影响各个区域的空气质量提供了详细的研究
往期推荐
不止 SHAP 力图：LIME 实现任意黑盒模型的单样本解释
微信号｜deep_ML
欢迎添加作者微信进入Python、ChatGPT群
进群请备注Python或AI进入相关群
无需科学上网、同步官网所有功能、使用无限制
如果你对类似于这样的文章感兴趣。
欢迎关注、点赞、转发~
个人观点，仅供参考