使用GeoJSON数据进行SHAP值地图可视化解释ML模型

背景

随着机器学习在环境和社会经济研究中的应用越来越广泛，人们对模型的透明性和可解释性的需求也在不断增加。本文将探讨如何使用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()