方法概述在空间数据分析中研究者通常需要回答三个递进的问题各解释变量对目标变量的影响是否存在空间非平稳性——由多尺度地理加权回归MGWR给出每个空间单元的局部系数。单因子及因子间的交互作用能在多大程度上解释目标变量的空间分异——由地理探测器计算 q 统计量并判断交互类型。在给定的数据集中哪种统计或机器学习模型能提供最可靠的目标变量预测——通过留一交叉验证LOOCV对七种候选模型进行精度比较自动选出最优模型。本代码将上述三个模块无缝衔接并额外计算全局Moran’s I以检验目标变量的空间自相关性全部结果输出为 Excel 表格与论文级统计图全程无需手动操作。框架本身不依赖于特定学科背景用户仅需按格式准备数据即可一键运行。2. 环境配置与数据约定依赖库import pandas as pd, numpy as np import geopandas as gpd import matplotlib.pyplot as plt import seaborn as sns from mgwr.gwr import MGWR from mgwr.sel_bw import Sel_BW from sklearn.model_selection import LeaveOneOut, cross_val_predict from sklearn.linear_model import LinearRegression, Ridge, Lasso from sklearn.svm import SVR from sklearn.ensemble import RandomForestRegressor, ExtraTreesRegressor, GradientBoostingRegressor from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error from esda.moran import Moran from libpysal.weights import Queen输入文件data.csv或.xlsx至少包含 3 列——经度lon、纬度lat、目标变量示例中为XHM以及若干特征列示例中Feat1~Feat7。shp/市.shp研究区面状矢量数据用于建立空间权重并计算 Moran’s I若不需要空间自相关分析可提供任一 shp 文件占位。自定义映射所有特征与目标变量的中文名称均在字典FEATURE_CN_MAP中配置例如FEATURE_CN_MAP { Feat1: 人均GDP, Feat2: 专利/万人, # 按实际情况修改 } TARGET_COL XHM框架会基于此映射自动生成图表中的变量标签使其具有可读性。3. 核心功能实现拆解3.1 MGWR 建模与局部系数提取MGWR 是地理加权回归GWR的扩展允许不同变量拥有各自的空间带宽从而更细致地刻画不同解释因素的空间作用尺度差异。处理流程对自变量进行Z-score 标准化消除量纲影响使用黄金分割搜索算法Sel_BW为每个变量含截距寻找最优带宽拟合 MGWR 模型提取每个样本点的预测值、残差及各变量局部系数计算模型拟合精度训练 R²、调整 R²、RMSE、MAE、AIC、AICc、BIC。关键代码片段def run_mgwr(data): x data[FEATURE_COLS].values y data[TARGET_COL].values.reshape((-1, 1)) coords list(zip(data[lon], data[lat])) scaler StandardScaler() x_scaled scaler.fit_transform(x) # 带宽搜索范围设置 n_obs len(data) n_features_with_const x_scaled.shape[1] 1 min_bw max(n_features_with_const 2, 10) max_bw n_obs - 1 selector Sel_BW(coords, y, x_scaled, multiTrue) bw selector.search(bw_minmin_bw, bw_maxmax_bw) # 模型拟合 model MGWR(coords, y, x_scaled, selector) results model.fit() pred results.predy.flatten() y_flat y.flatten() # 提取局部系数 params results.params result_data data.copy() result_data[Intercept_coef] params[:, 0] for idx, col in enumerate(FEATURE_COLS): result_data[f{col}_coef] params[:, idx 1] result_data[f{FEATURE_CN_MAP[col]}_局部系数] params[:, idx 1] # 精度指标 metrics pd.DataFrame([{ 模型: MGWR, 训练R2: r2_score(y_flat, pred), 调整R2: adjusted_r2(y_flat, pred, len(FEATURE_COLS)), RMSE: rmse(y_flat, pred), MAE: mean_absolute_error(y_flat, pred) }]) return result_data, metrics, bw_table, coef_summary注意若样本量过小小于变量数 2MGWR 无法进行有效的带宽搜索代码会提前抛出清晰错误提示。3.2 地理探测器单因子与交互探测地理探测器基于“如果自变量对因变量有重要影响那么自变量与因变量的空间分布应趋于一致”的假设通过层内方差和总方差的比值构造q 统计量$$q 1 - \frac{\sum_{h1}^{L} N_h \sigma_h^2}{N \sigma^2}$$式中$L$ 为分类数$N_h$ 与 $\sigma_h^2$ 分别为第 $h$ 层的样本数与方差。q 值范围为 $[0, 1]$值越大表示该因子对目标变量的解释力越强。连续变量离散化由于地理探测器要求输入为类型变量代码默认按分位数将连续变量分为 5 类def discretize(series, bins5): return pd.qcut(series, qbins, duplicatesdrop)交互探测通过将两因子的分类标签交叉拼接生成新的组合分类计算 $q_{12}$并与 $q_1$、$q_2$ 比较判断交互类型def classify_interaction(q1, q2, q12): if q12 min(q1, q2): return 非线性减弱 elif q12 max(q1, q2): return 单因子非线性减弱 elif np.isclose(q12, q1q2): return 独立 elif q12 q1q2: return 双因子增强 else: return 非线性增强最终输出包含 q 值排序表、对称交互矩阵、以及拆分的交互类型明细表。3.3 多模型预测精度自动比较考虑到研究样本量通常有限代码采用留一交叉验证LOOCV评估模型泛化能力。参与比较的模型涵盖线性模型OLS、Ridge、Lasso、支持向量回归RBF 核以及三种集成树模型随机森林、极端随机树、梯度提升树。def compare_prediction_models(data): x data[FEATURE_COLS].values y data[TARGET_COL].values loo LeaveOneOut() models { 线性回归: Pipeline([...]), Ridge: Pipeline([...]), # ... } for name, model in models.items(): cv_pred cross_val_predict(model, x, y, cvloo) model.fit(x, y) train_pred model.predict(x) # 记录训练R2、LOOCV_R2、RMSE、MAE accuracy pd.DataFrame(rows).sort_values(LOOCV_R2, ascendingFalse) best_model fitted_models[accuracy.iloc[0][模型]] return accuracy, best_pred_table, feature_importance最终按LOOCV_R2降序排列将最优模型的全量预测值、LOOCV 预测值及残差单独输出若最优模型为树模型则同时输出变量重要性。3.4 空间自相关检验基于 Queen 邻接权重矩阵计算目标变量的全局 Moran’s I作为空间依赖性的诊断依据def compute_moran(data, gdf): merged gdf.merge(data[[CITY_FIELD, TARGET_COL]], onCITY_FIELD) w Queen.from_dataframe(merged) w.transform r moran Moran(merged[TARGET_COL], w) return pd.DataFrame([{Moran_I: moran.I, P_value: moran.p_sim, Z_score: moran.z_sim}])4. 可视化输出与结果导出make_charts函数自动生成 11 张适用于论文展示的统计图表均为非地图类图号内容图表类型01地理探测器 q 值排序横向条形图02影响因子解释力棒棒糖图克利夫兰点图03交互探测热力图对称矩阵热力图04Top 12 交互组合 q 值分组条形图05七种模型精度对比水平点线图06最优模型真实值 vs 预测值散点图回归线07最优模型 LOOCV 残差分布直方图核密度08最优模型变量重要性若可得条形图09变量间相关系数热力图相关矩阵热力图10MGWR 局部系数分布箱线图箱线图11各城市局部系数热力图矩阵热力图所有图表以 600 dpi 分辨率保存至输出目录文件名自动处理特殊字符。同时所有中间结果表合并写入一个 Excel 工作簿并额外导出主要结果q 值、交互矩阵、精度对比、预测值等的 CSV 副本。5. 使用说明与注意事项路径与数据准备修改EXCEL_PATH、SHP_PATH和OUT_DIR三个变量至实际路径确保 CSV 中经度、纬度、目标变量及特征列名与脚本一致。特征变量灵活配置只需在FEATURE_CN_MAP字典中增删条目脚本会自动适配后续分析。样本量要求MGWR 要求样本量 ≥ 变量数含截距 2否则抛出错误地理探测器在样本量很小时 q 值可能不稳定。离散化层数discretize函数默认将连续变量分 5 类可根据研究需要调整bins参数。运行环境推荐 Python 3.8 及以上版本依赖库可通过以下命令安装pip install geopandas mgwr esda libpysal seaborn openpyxl scikit-learn自定义修改若希望调整模型超参数或增加新的比较模型可直接在compare_prediction_models函数的models字典中修改。6. 完整代码以下代码已去除特定领域背景直接复制即可运行。用户仅需配置输入路径与特征映射字典。# -*- coding: utf-8 -*- MGWR 地理探测器 精度提升框架 -------------------------------- 输出内容 1. MGWR 模型精度、带宽、预测值、残差、局部系数数据表 2. 地理探测器 q 值、交互探测矩阵与交互类型数据表 3. 非地图型论文图表q 值图、交互热力图、精度对比图、拟合图、残差图等 4. 自动比较多种预测模型给出当前数据下精度最高的模型 import os import warnings from pathlib import Path warnings.filterwarnings(ignore) import geopandas as gpd import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np import pandas as pd import seaborn as sns from esda.moran import Moran from libpysal.weights import Queen from mgwr.gwr import MGWR from mgwr.sel_bw import Sel_BW from sklearn.ensemble import ExtraTreesRegressor, GradientBoostingRegressor, RandomForestRegressor from sklearn.linear_model import Lasso, LinearRegression, Ridge from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score from sklearn.model_selection import LeaveOneOut, cross_val_predict from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler from sklearn.svm import SVR # # 1. 参数设置根据实际数据修改此部分 # EXCEL_PATH rdata.csv # 属性数据路径 SHP_PATH rshp\市.shp # 矢量面数据路径 OUT_DIR r./output # 输出文件夹 CITY_FIELD name # 城市名称字段若不存在会按经纬度匹配 LON_COL lon LAT_COL lat TARGET_COL XHM # 目标变量 # 特征变量中英文映射可按需增删 FEATURE_CN_MAP { Feat1: 人均GDP, Feat2: 专利/万人, Feat3: 对外开放度, Feat4: 产业高级化, Feat5: 科技支出占比, Feat6: 交通可达性, Feat7: 普惠金融指数, } TARGET_CN_NAME 目标变量 FEATURE_COLS list(FEATURE_CN_MAP.keys()) FEATURE_NAMES [FEATURE_CN_MAP[col] for col in FEATURE_COLS] os.makedirs(OUT_DIR, exist_okTrue) # 图形全局设置 mpl.rcParams[font.sans-serif] [SimHei, Microsoft YaHei, Arial Unicode MS] mpl.rcParams[axes.unicode_minus] False mpl.rcParams[figure.dpi] 120 mpl.rcParams[savefig.dpi] 600 mpl.rcParams[axes.linewidth] 1.1 sns.set_theme(stylewhitegrid, fontSimHei) INVALID_FILENAME_CHARS str.maketrans({char: _ for char in r\/:*?|}) def safe_filename(name): return name.translate(INVALID_FILENAME_CHARS)