本次分析所用的數(shù)據(jù)可從DVRPC的數(shù)據(jù)集中獲取，數(shù)據(jù)格式為GeoJSON，GeoJSON是一種常見的地理數(shù)據(jù)結(jié)構(gòu)編碼格式，可以輕松與Python中的GeoPandas等地理空間庫集成
使用該數(shù)據(jù)集預(yù)測空氣污染水平，具體為PM2.5濃度（pm25f_pfs），并基于多個環(huán)境和社會經(jīng)濟(jì)特征作為預(yù)測變量，選擇的預(yù)測變量（自變量）包括：
在本次分析中，采用XGBoost機(jī)器學(xué)習(xí)模型，這是一種因其魯棒性和靈活性而被廣泛應(yīng)用于回歸任務(wù)的算法，為了增強(qiáng)模型的可解釋性，我們使用了SHAP值，這是一種常用于解釋機(jī)器學(xué)習(xí)模型輸出的方法，當(dāng)然這里作者不會過多去考慮模型精確度重點在于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格式的地理空間數(shù)據(jù)文件，顯示其前幾行數(shù)據(jù)以查看字段信息
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比例劃分?jǐn)?shù)據(jù)集，確保隨機(jī)種子設(shè)置以便重現(xiàn)性
import xgboost as xgb
from sklearn.model_selection import GridSearchCV
# XGBoost回歸模型參數(shù)
params_xgb = {
'learning_rate': 0.02, # 學(xué)習(xí)率，控制每一步的步長，用于防止過擬合。典型值范圍：0.01 - 0.1
'booster': 'gbtree', # 提升方法，這里使用梯度提升樹（Gradient Boosting Tree）
'objective': 'reg:squarederror', # 損失函數(shù)，這里使用平方誤差，適用于回歸任務(wù)
'max_leaves': 127, # 每棵樹的葉子節(jié)點數(shù)量，控制模型復(fù)雜度
'verbosity': 1, # 控制 XGBoost 輸出信息的詳細(xì)程度，0表示無輸出，1表示輸出進(jìn)度信息
'seed': 42, # 隨機(jī)種子，用于重現(xiàn)模型的結(jié)果
'nthread': -1, # 并行運(yùn)算的線程數(shù)量，-1表示使用所有可用的CPU核心
'colsample_bytree': 0.6, # 每棵樹隨機(jī)選擇的特征比例，用于增加模型的泛化能力
'subsample': 0.7, # 每次迭代時隨機(jī)選擇的樣本比例，用于增加模型的泛化能力
'eval_metric': 'rmse' # 評價指標(biāo)，這里使用均方根誤差（rmse）
}
# 初始化XGBoost回歸模型
model_xgb = xgb.XGBRegressor(**params_xgb)
# 定義參數(shù)網(wǎng)格，用于網(wǎng)格搜索
param_grid = {
'n_estimators': [100, 200, 300, 400, 500], # 樹的數(shù)量
'max_depth': [3, 4, 5, 6, 7], # 樹的深度
'learning_rate': [0.01, 0.02, 0.05, 0.1], # 學(xué)習(xí)率
}
# 使用GridSearchCV進(jìn)行網(wǎng)格搜索和k折交叉驗證
grid_search = GridSearchCV(
estimator=model_xgb,
param_grid=param_grid,
scoring='neg_mean_squared_error', # 評價指標(biāo)為負(fù)均方誤差
cv=5, # 5折交叉驗證
n_jobs=-1, # 并行計算
verbose=1 # 輸出詳細(xì)進(jìn)度信息
)
# 訓(xùn)練模型
grid_search.fit(X_train, y_train)
# 輸出最優(yōu)參數(shù)
print("Best parameters found: ", grid_search.best_params_)
print("Best RMSE score: ", (-grid_search.best_score_) ** 0.5) # 還原RMSE
# 使用最優(yōu)參數(shù)訓(xùn)練模型
best_model_xgboost = grid_search.best_estimator_配置XGBoost回歸器的超參數(shù)，并使用網(wǎng)格搜索和交叉驗證來識別最佳參數(shù)。使用最佳參數(shù)在訓(xùn)練數(shù)據(jù)上重新訓(xùn)練模型
import shap
# 構(gòu)建 shap解釋器
explainer = shap.TreeExplainer(best_model_xgboost)
# 計算shap值
shap_values = explainer.shap_values(X)
# 特征標(biāo)簽
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值分布，清晰顯示了每個特征對預(yù)測PM2.5水平的正向或負(fù)向影響及其幅度
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()
根據(jù)SHAP值的平均絕對值對特征進(jìn)行排序，突出顯示了模型中最具影響力的預(yù)測變量，其中wf_pfs（工作參與率）是對模型輸出影響最大的特征，其次是ebf_pfs（暴露于洪水的比例）和hsef（住房壓力因素）
數(shù)據(jù)整理
# 將 SHAP 值轉(zhuǎn)換為 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值，這為區(qū)域內(nèi)模型預(yù)測差異的洞察提供了幫助，通過將SHAP值與原始GeoJSON數(shù)據(jù)結(jié)合，創(chuàng)建一個包含詳細(xì)空間可視化SHAP值的地理空間數(shù)據(jù)集（gdf_with_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))
# 設(shè)置標(biāo)題
plt.title("SHAP Value for ebf_pfs")
plt.savefig("SHAP_3.pdf", format='pdf', dpi=1200, bbox_inches='tight')
plt.show()
這里先單獨繪制一個特征的SHAP值地圖，以了解特征如何在不同地理區(qū)域影響PM2.5預(yù)測水平，ebf_pfs（暴露于洪水的人口百分比）的SHAP地圖顯示了這一因素在特定區(qū)域增加或減少預(yù)測PM2.5水平的影響
總結(jié)而言，這張圖展示了特征ebf_pfs在不同地理區(qū)域?qū)δＰ皖A(yù)測結(jié)果的影響，為環(huán)境決策提供了更具針對性的區(qū)域洞察。需要注意的是，這里的重點在于如何繪制SHAP地圖，因此并未深入研究模型的精確度是否足夠應(yīng)用于實際場景
# 提取所有以 'shap_' 開頭的列名
shap_columns = [col for col in gdf_with_shap.columns if col.startswith('shap_')]
# 創(chuàng)建一個畫布，設(shè)置子圖的網(wǎng)格布局
num_cols = len(shap_columns)
fig, axes = plt.subplots(nrows=(num_cols + 2) // 2, ncols=2, figsize=(18, 6 * ((num_cols + 2) // 2)))
# 展平 axes 數(shù)組以便于迭代
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]) # 在對應(yīng)子圖上繪制
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()
為了提供全面的視角，使用網(wǎng)格布局繪制了所有特征的SHAP地圖，為不同的社會經(jīng)濟(jì)和環(huán)境因素如何影響各個區(qū)域的空氣質(zhì)量提供了詳細(xì)的研究
文章轉(zhuǎn)自微信公眾號@Python機(jī)器學(xué)習(xí)AI