地理加权回归模型——一种局部的空间回归分析方法

地理加权回归（GWR）是一种局部的空间回归分析方法，它允许模型参数在空间上变化，从而能够捕捉到空间数据的局部空间非平稳性。GWR模型的基本思想是在回归分析中引入空间权重，使得模型能够根据地理位置的邻近程度对观测值进行加权。

1. GWR模型的基本形式

$y_{i} = \beta (u_{i})\cdot X_{i}+\varepsilon _{i}$

其中， $y{_{i}}$ 是第 $i$ 个观测变量， $X_{i}$ 是第 $i$ 个观测值的自变量向量， $\beta (u_{i})$ 是空间上变化的回归系数向量， $u_{i}$ 表示第 $i$ 个观测值的地理位置， $\varepsilon _{i}$ 是第 $i$ 个观测值的误差项。

2. 空间权重矩阵

GWR模型中的关键部分是空间权重矩阵 $W$ ，它根据观测值之间的空间距离来定义权重。权重通常由空间权重函数 $W(h_{ij})$ 计算得出，其中 $h_{ij}$ 是第 $i$ 个和第 $j$ 个观测值之间的空间距离。权重函数可以是高斯函数、双平方指数函数等。

3. 模型参数估计

GWR模型的参数估计通常采用最大似然估计方法。通过最大化似然函数，可以得到回归系数的估计值。

4. 模型优缺点

4.1 优点

局部空间非平稳性：GWR能够捕捉到空间数据的局部空间非平稳性，即模型参数在空间上的变化。
灵活性：模型允许参数随位置变化，提供了一种灵活的方式来分析空间数据。
可视化：GWR模型的局部估计可以生成空间效应的可视化，帮助理解不同区域的空间关系。

4.2 缺点

模型复杂性：由于参数空间变化，模型可能会变得复杂，难以解释。
过度拟合：如果空间权重过大，可能会导致模型对局部数据过度拟合。
选择带宽：需要选择合适的带宽（bandwidth），这可能需要依赖于经验或交叉验证

GWR模型在地理学、环境科学、流行病学等领域中有着广泛的应用，尤其是在分析空间数据时，能够提供比传统全局回归模型更深入的洞察。

5. python实现

我们在这里给出一段Python代码，使用了各个省份的省会的经纬度作为空间位置的衡量，使用了bi-square函数作为核函数，使用了AIC(赤池信息准责)和BIC（贝叶斯信息准则）确定最优带宽，来帮助读者更好的理解和实现：

import os
import pandas as pd
from mgwr.gwr import GWR
from mgwr.sel_bw import Sel_BW
import numpy as np
from sklearn.utils import check_array  # 添加这一行来导入check_array
from sklearn.utils.validation import check_array

import os
os.environ['PYTHONUTF8'] = '1'
os.environ['PYTHONIOENCODING'] = 'utf-8'

# 步骤 1: 加载数据
# 使用'openpyxl'引擎来读取Excel文件
data = pd.read_excel('shu_ju_ji.xlsx', engine='openpyxl')




# 检查数据是否包含“经度”和“纬度”列
#print(data.dtypes)

# 将非数值列转换为数值类型(float或int)
# 这里需要确保所有自变量的列名都是正确的，并且存在于data中
numeric_columns = ['GDP', '人口密度', '城乡人均收入比值', '分省政府网站数量',
                   '网民人数', '互联网普及率', '互联网宽带接入端口数量', '光缆线路长度',
                   '老年人占比', '性别比', '流量', '高中及以上学历占比',
                   '本科及以上学历占比', '经度值', '纬度值']

# 确保所有自变量都是数值类型
for col in numeric_columns:
    data[col] = pd.to_numeric(data[col], errors='coerce')  # errors='coerce' 会将无法转换的值设置为NaN

# 移除包含NaN的行，如果需要的话
data = data.dropna()

# 定义自变量的列名
independent_vars = ['GDP']  # 示例中只选择'GDP'作为自变量
dependent_var = 'RSR拟合值'

# 确保所有自变量和因变量都是数值类型
for col in independent_vars + [dependent_var]:
    data[col] = pd.to_numeric(data[col], errors='coerce')  # errors='coerce' 会将无法转换的值设置为NaN

# 移除包含NaN的行
data = data.dropna(subset=independent_vars + [dependent_var])

# 提取自变量和因变量
X = data[independent_vars].values.astype('float64')  # 确保X是2D数组
y = data[dependent_var].values.astype('float64').reshape((-1, 1))  # 确保y是2D数组

# 检查X和y的大小
print(f"X shape: {X.shape}")
print(f"y shape: {y.shape}")

# 如果X或y是空的，处理数据问题
if X.size == 0 or y.size == 0:
    raise ValueError("One of the input arrays is empty. Please check the data.")

# 提取经纬度用于地理权重
coords = data[['经度值', '纬度值']].values

# 选择带宽
bw = Sel_BW(coords, y, X, fixed=False, kernel='bisquare').search()
print(f'Optimal bandwidth: {bw}')

# 拟合GWR模型

model = GWR(coords, y, X, bw=bw, fixed=False, kernel='bisquare')
results = model.fit()

# 输出地理加权回归的摘要统计信息
print(results.summary())

# 获取局部估计的参数
local_params = results.params

# 打印每个观测点的局部估计参数
print("\nLocally Estimated Parameters for Each Observation:")
for i, params in enumerate(local_params):
    print(f"Observation {i+1}: "
          f"Intercept: {params[0]:.4f}, "
          f"GDP Coefficient: {params[1]:.4f}")

# 如果您需要构建每个观测点的局部回归方程，可以这样做：
print("\nLocal Regression Equations for Each Observation:")
for i, params in enumerate(local_params):
    eqn = f"RSR拟合值 = {params[0]:.4f} + {params[1]:.4f} * GDP"
    print(f"Observation {i+1}: {eqn}")