はじめに

　『Pythonで学ぶ空間データサイエンス入門』の独学ノートです。本の内容から寄り道してアレコレ考えます。
　本を読んだ上で補助的に読んでください。　

　この記事では、GWRモデルの最小二乗法について、図を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
4.3.4-7 GWRモデルの最小二乗法の可視化
- 最小二乗法の図
参考文献
おわりに

4.3.4-7 GWRモデルの最小二乗法の可視化

　地理的加重回帰モデル(GWRモデル・geographically weighted regression model)に対する最小二乗法(OLS・ordinary least square)の計算を図で確認します。
　地理的加重回帰モデルの定義や計算式については「4.3.4-7：GWRモデルの最小二乗法の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」、線形回帰モデルや最小二乗法については「【Python】3.2.1-3：線形回帰モデルの最小二乗法の可視化：1次元の場合【空間データサイエンス入門のノート】 - からっぽのしょこ」「【Python】3.2.1-3：線形回帰モデルの最小二乗法の可視化：2次元の場合【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

　利用するライブラリを読み込みます。

# ライブラリを読込
import numpy as np
import libpysal as ps
import geopandas as gpd
import matplotlib.pyplot as plt

最小二乗法の図

　まずは、GWRモデルに対するOLSの定義や性質を可視化します。

観測データの図

　例として利用する空間データを読み込んで前処理を行います。詳しくは本を参照してください。

　Georgiaデータセットを読み込みます。

# 空間データを読込
geom_df = gpd.read_file(ps.examples.get_path('G_utm.shp'))

# 利用する列を抽出
geom_df = geom_df[['PctFB', 'PctBlack', 'PctRural', 'PctBach', 'X', 'Y', 'geometry']]
geom_df

	PctFB	PctBlack	PctRural	PctBach	X	Y	geometry
0	0.64	20.76	75.6	8.2	941396.6	3521764	POLYGON ((931869.062 3545540.5, 934111.625 354...
1	1.58	26.86	100.0	6.4	895553.0	3471916	POLYGON ((867016.312 3482416, 884309.375 34826...
2	0.27	15.42	61.7	6.6	930946.4	3502787	POLYGON ((914656.875 3512190, 924718.375 35125...
3	0.11	51.67	100.0	9.4	745398.6	3474765	POLYGON ((744258.625 3480598.5, 765025.062 348...
4	1.43	42.39	42.7	13.3	849431.3	3665553	POLYGON ((832974.188 3677273.5, 834048.688 367...
...	...	...	...	...	...	...	...
154	2.55	4.06	70.0	12.0	686891.4	3855274	POLYGON ((684405.938 3873304.5, 684661.562 387...
155	0.09	31.76	100.0	7.6	838551.5	3538547	POLYGON ((845655.625 3557787.75, 846281.25 355...
156	0.43	45.94	59.6	10.4	891228.5	3749769	POLYGON ((902279.938 3768771.75, 906628.562 37...
157	0.29	41.99	100.0	8.8	858796.9	3637891	POLYGON ((867125.062 3651962.75, 867417.75 365...
158	0.59	30.71	71.1	6.3	801018.1	3487328	POLYGON ((798619.812 3525322, 799558.25 352473...

159 rows × 7 columns

　libpysal ライブラリに含まれるジョージア州の空間データを利用します。
　(利用する列を取り出す必要はないのですが、後ほど列を追加した際に確認しにくくなるので、不要な列を取り除いておきます。)

　説明変数と被説明変数を設定します。

# 説明変数を取得
tmp_X = geom_df[['PctFB', 'PctBlack', 'PctRural']].values.astype(np.float64)

# データ数・パラメータ数を取得
N, K = tmp_X.shape
print(N, K)

# 定数項用の列を結合
X = np.hstack(
    [np.ones(shape=(N, 1)), tmp_X]
)
print(X[:5].round(1))
print(X.shape)

# 被説明変数を取得
y = geom_df['PctBach'].values.astype(np.float64)

# パラメータ数を取得
K = X.shape[1] - 1

159 3
[[  1.    0.6  20.8  75.6]
 [  1.    1.6  26.9 100. ]
 [  1.    0.3  15.4  61.7]
 [  1.    0.1  51.7 100. ]
 [  1.    1.4  42.4  42.7]]
(159, 4)

　Pandasデータフレームから各次元 $k$ の説明変数 $x_{1k}, \cdots, x_{Nk}$ として利用する列を取り出し、 $0$ 番目の次元(定数項用)の値 $x_{j0} = 1$ と列方向に結合して、各地点 $j$ の説明変数 $\mathbf{x}_j = (x_{j0}, x_{j1}, \cdots, x_{jK})$ とします。また、 $N$ 個の地点を行方向に結合して説明変数 $\mathbf{X} = (x_{10}, \cdots, x_{NK})$ とします。
　同様に、 $N$ 個の被説明変数 $\mathbf{y} = (y_1, \cdots, y_N)$ として用いる列を取り出して、NumPy配列に変換します。
　データフレームの行数がデータ数 $N$ 、取り出した列数がパラメータ数 $K$ に対応します。

　説明変数と被説明変数のコロプレス図を作成します。

作図コード(クリックで展開)

# 列名を指定
col_label_lt = ['PctBach', 'PctFB', 'PctBlack', 'PctRural']

# 行数・列数を設定
R, C = 2, 2

# 観測データを作図
fig, axes = plt.subplots(nrows=R, ncols=C, 
                         figsize=(11, 10), dpi=100, facecolor='White', 
                         constrained_layout=True)
fig.suptitle('observation data', fontsize=20)

# 次元番号を初期化
k = 0

# 変数を描画
for r in range(R):
    for c in range(C):
        
        # 列名を取得
        col_label = col_label_lt[k]

        # ラベル用の文字列を作成
        if k == 0:
            var_label = '$y_j$'
        else:
            var_label = f'$x_{{j{k}}}$'

        # 変数を描画
        ax = axes[r, c]
        geom_df.plot(ax=ax, column=col_label, 
                     edgecolor='white', linewidth=0.5, 
                     legend=True, legend_kwds={'label': var_label}) # y_j, x_jk
        ax.set_xlabel('$u_j$')
        ax.set_ylabel('$v_j$')
        ax.set_title(col_label)
        ax.grid()
        ax.set_aspect('equal', adjustable='box')

        # 次元番号をカウント
        k += 1

plt.show()

　 $N$ 個の地点について、被説明変数 $\mathbf{y} = (y_1, \cdots, y_j, \cdots, y_N)$ と各次元の説明変数 $(x_{1k}, \cdots, x_{jk}, \cdots, x_{Nk})$ の値をそれぞれコロプレス図で示します。

　説明変数と被説明変数を標準化します。

# 変数を標準化
tmp_X = (tmp_X - tmp_X.mean(axis=0)) / tmp_X.std(axis=0)
y = (y - y.mean()) / y.std()

# 標準化した変数を格納
geom_df[['PctFB', 'PctBlack', 'PctRural']] = tmp_X
geom_df['PctBach'] = y

# 定数項用の列を結合
X = np.hstack(
    [np.ones(shape=(N, 1)), tmp_X]
)
print(X[:5].round(1))
print(X.mean(axis=0))
print(X.std(axis=0))
print(y[:5].round(1))
print(y.mean())
print(y.std())

[[ 1.  -0.4 -0.4  0.2]
 [ 1.   0.4 -0.   1.1]
 [ 1.  -0.7 -0.7 -0.3]
 [ 1.  -0.8  1.4  1.1]
 [ 1.   0.2  0.9 -1. ]]
[ 1.00000000e+00 -1.34064667e-16 -2.23441112e-17 -8.04388003e-16]
[0. 1. 1. 1.]
[-0.5 -0.8 -0.8 -0.3  0.4]
2.2344111187424534e-16
1.0

　この例では作図(の配色の共通化)用に、 $\mathbf{X}$ の次元(列)ごとと $\mathbf{y}$ を標準化しておきます。ただし、0番目の次元(定数項用の変数)については、標準偏差が0のため、0除算になり計算できません。

　標準化した説明変数と被説明変数をグラフで確認します。

　 $\mathbf{X}$ の各次元(グラフ)と $\mathbf{y}$ ごとに標本平均が0、標本標準偏差が1になるように標準化しました。各変数の値は変わりますが、大小関係などは変わりません。

　位置データを作成します。

# 位置情報を取得
coord_arr = np.stack(
    [geom_df['X'].values.astype(np.float64),
     geom_df['Y'].values.astype(np.float64)], 
    axis=1
)
print(coord_arr[:5])
print(coord_arr.shape)

[[ 941396.6 3521764. ]
 [ 895553.  3471916. ]
 [ 930946.4 3502787. ]
 [ 745398.6 3474765. ]
 [ 849431.3 3665553. ]]
(159, 2)

　Pandasデータフレームから経度 X 列と緯度 Y 列(に対応する値?これは何ですか?)を取り出してNumPy配列に格納します。
　各地点 $j$ の経度を $u_i$ 、経度を $v_i$ として、位置を $(u_i, v_i)$ で表します。

　ジョージア州における各地点のグラフを作成します。

# 位置データを作図
fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', 
                       constrained_layout=True)
geom_df.plot(ax=ax, color='white', edgecolor='black', linewidth=0.5) # 行政区界
ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], 
           color='black', s=5) # 各地点
ax.set_xlabel('$u_j$')
ax.set_ylabel('$v_j$')
ax.set_title(f'$N = {N}$', loc='left')
fig.suptitle('State of Georgia ', fontsize=20)
ax.grid()
ax.set_aspect('equal', adjustable='box')
plt.show()

　区域ごとの点を各地域の位置として用います。

　基準地点と各地点の距離の関係をアニメーションで確認します。

作図コード(クリックで展開)

# 基準地点を指定
i = 0

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', 
                       constrained_layout=True)
fig.suptitle('State of Georgia ', fontsize=20)

# 作図処理を定義
def update(j):
    
    # 前フレームのグラフを初期化
    ax.cla()

    # 位置を抽出
    u_i, v_i = coord_arr[i]
    u_j, v_j = coord_arr[j]

    # 距離を計算
    dist = np.sqrt(np.sum((coord_arr[j] - coord_arr[i])**2))

    # 位置データを描画
    geom_df.plot(ax=ax, color='white', edgecolor='black', linewidth=0.5) # 行政区界
    ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], 
               color='black', s=5) # 各地点
    ax.scatter(x=u_j, y=v_j, s=50, color='black', 
               label=f'$(u_j, v_j) = ({u_j:.1f}, {v_j:.1f})$') # 各地点
    ax.scatter(x=u_i, y=v_i, s=50, color='red', 
               label=f'$(u_i, v_i) = ({u_i:.1f}, {v_i:.1f})$') # 基準地点
    ax.set_xlabel('$u_j$')
    ax.set_ylabel('$v_j$')
    ax.set_title(f'$N = {N}, i = {i+1}, j = {j+1}$', loc='left')
    ax.grid()
    ax.legend(loc='upper left')
    ax.set_aspect('equal', adjustable='box')

# 動画を作成
ani = FuncAnimation(fig=fig, func=update, frames=N, interval=500)

# 動画を書出
ani.save(
    filename='georgia_dist.mp4', 
    progress_callback = lambda i, n: print(f'frame: {i} / {n}')
)

　基準地点 $i$ を赤色の点、各地点 $j$ を黒色の点、2点間の距離 $d_{ij}$ を青色の線分で示します。
　GeoSeries.centroid.plot() は重心の点を描画するので、X, Y 列の点とはズレが生じます。

重み行列の図

　重み行列については「【Python】4.3.3：GWRモデルの重み行列の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」、重み行列の作成については「【Python】4.3.2-3：GWRモデルの重み行列の実装【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。
　アニメーションの作図コードについては「ols.py」を参照してください。

　バンド幅を指定して、重み行列を作成します。

# バンド幅を指定:(バンド幅を固定する場合)
#b = 270000.0

# バンド幅内の地点数を指定:(地点数を固定する場合)
b_idx = 116

# 受け皿を初期化
b_lt = []
W_lt = []

# 基準地点ごとに処理
for i in range(N):
    
    # 距離を計算
    dist_vec = np.sqrt(np.sum((coord_arr - coord_arr[i])**2, axis=1))

    # バンド幅を設定:(地点数を固定する場合)
    b = np.sort(dist_vec)[b_idx]
    
    # 対角要素を計算:(bi-square型)
    w_diag = (1.0 - (dist_vec / b)**2)**2
    w_diag[dist_vec >= b] = 0.0
    
    # 重み行列を作成
    W_i = np.diag(w_diag)

    # 結果を格納
    b_lt.append(b)
    W_lt.append(W_i)
print(W_lt[0].round(2))
print(W_lt[0].shape)
print(len(W_lt))

[[1.   0.   0.   ... 0.   0.   0.  ]
 [0.   0.88 0.   ... 0.   0.   0.  ]
 [0.   0.   0.99 ... 0.   0.   0.  ]
 ...
 [0.   0.   0.   ... 0.07 0.   0.  ]
 [0.   0.   0.   ... 0.   0.53 0.  ]
 [0.   0.   0.   ... 0.   0.   0.52]]
(159, 159)
159

　基準地点 $i$ ごとに、各地点 $j$ とのユークリッド距離 $d_{ij} = \sqrt{(u_j - u_i)^2 + (v_j - v_i)^2}$ を計算し、重み行列 $\mathbf{W}_i = (w_{i11}, \cdots, w_{iNN})$ を作成して、リストに格納していきます。

　全ての基準地点で共通のバンド幅 $b$ を指定します。この場合、基準地点ごとにバンド幅の範囲内の地点数が変わります。
　または、バンド幅の範囲内に含まれる地点数 b_idx を指定して、基準地点ごとにバンド幅 $b$ を設定します。基準地点と各地点の距離を昇順に並べ替えて b_idx 番目の距離をバンド幅 $b$ として用います。これは mgwr ライブラリにおける処理と同じ操作です。
　バンド幅を指定(固定)する場合は、コメントアウトしてください。

　重み行列の対角要素のコロプレス図を作成します。

作図コード(クリックで展開)

# 基準地点を指定
i = 0

# 重み行列を取得
W_i = W_lt[i]

# 重みを格納
geom_df['weight'] = np.diag(W_i)

# バンド幅を取得
b = b_lt[i]

# 基準地点の位置を取得
u_i, v_i = coord_arr[i]

# バンド幅の範囲を計算
t_vec = np.linspace(start=0.0, stop=2.0*np.pi, num=361) # ラジアン
u_vec = u_i + b * np.cos(t_vec)
v_vec = v_i + b * np.sin(t_vec)

# グラフサイズを設定
u = 50000 # 切り下げ・切り上げの単位を指定
u_min = np.floor(geom_df.bounds.minx.min() /u)*u
u_max = np.ceil(geom_df.bounds.maxx.max()  /u)*u
v_min = np.floor(geom_df.bounds.miny.min() /u)*u
v_max = np.ceil(geom_df.bounds.maxy.max()  /u)*u

# 重みのコロプレス図を作成
fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', 
                       constrained_layout=True)
geom_df.plot(ax=ax, column='weight', 
             cmap='viridis', vmin=0.0, vmax=1.0, 
             edgecolor='white', linewidth=0.5, 
             legend=True, legend_kwds={'label': '$w_{ijj}$'}) # 対角要素
ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], 
           color='black', s=5) # 各地点
ax.scatter(x=u_i, y=v_i, color='red', s=50, 
           label='$(u_i, v_i)$') # 基準地点
ax.plot(u_vec, v_vec, 
        color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.set_xlabel('$u_j$')
ax.set_ylabel('$v_j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('weight matrix: (bi-square kernel)', fontsize=20)
fig.legend()
ax.grid()
ax.set_xlim(xmin=u_min, xmax=u_max)
ax.set_ylim(ymin=v_min, ymax=v_max)
ax.set_aspect('equal', adjustable='box')
plt.show()

　基準地点 $i$ の重み行列の対角要素 $\mathrm{diag}(\mathbf{W}_i) = (w_{i11}, \cdots, w_{ijj}, \cdots, w_{iNN})$ の値をコロプレス図で示します。
　基準地点 $i$ の重みが $w_{iii} = 1$ で、離れた地点ほど値が0に近付きます。この例では、bi-square型の重み関数により重み付けしているので、バンド幅の範囲外は0になります。

　重み付けした説明変数と被説明変数をグラフで確認します。

　上段は元の変数 $y_j, x_{jk}$ 、下段は重み付けした変数 $w_{ijj} y_j, w_{ijj} x_{jk}$ を表します。分かりやすいように、値が0のとき白色になるように調整しています。

　基準地点ごとに各地点の重みが変わるので、各地点のモデルへの影響の度合いも変わります。

係数パラメータの推定

　係数パラメータの推定値を計算します。

# 係数パラメータの推定値を初期化
Beta_hat = np.zeros(shape=(N, K+1))

# 基準地点ごとに処理
for i in range(N):

    # 重み行列を取得
    W_i = W_lt[i]
    
    # 係数パラメータの推定値を計算
    Beta_hat[i] = np.linalg.inv(X.T @ W_i @ X) @ X.T @ W_i @ y
print(Beta_hat[:5].round(2))
print(Beta_hat.shape)

[[-0.23  0.23  0.06 -0.43]
 [-0.28  0.17  0.1  -0.41]
 [-0.25  0.2   0.07 -0.43]
 [-0.23  0.15  0.05 -0.36]
 [ 0.19  0.72 -0.17 -0.24]]
(159, 4)

　全ての基準地点をまとめた係数パラメータの推定値 $\hat{\boldsymbol{\beta}} = (\hat{\beta}_{10}, \cdots, \hat{\beta}_{NK})$ として、基準地点 $i$ ごとに、係数パラメータの推定値 $\hat{\boldsymbol{\beta}} = (\hat{\beta}_{i0}, \hat{\beta}_{i1}, \cdots, \hat{\beta}_{iK})$ を $\hat{\boldsymbol{\beta}}_i = (\mathbf{X}^{\top} \mathbf{W}_i \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{W}_i \mathbf{y}$ で計算して、2次元配列に格納していきます。

　係数パラメータの推定値のコロプレス図を作成します。

作図コード(クリックで展開)

# 列名を設定
col_var_label_lt   = ['intercept', 'PctFB', 'PctBlack', 'PctRural']
col_param_label_lt = ['param_itercept', 'param_PctFB', 'param_PctBlack', 'param_PctRural']

# 係数パラメータの推定値を格納
geom_df[col_param_label_lt] = Beta_hat

# 行数・列数を設定
R, C = 2, 2

# 係数パラメータを作図
fig, axes = plt.subplots(nrows=R, ncols=C, 
                         figsize=(11, 10), dpi=100, facecolor='White', 
                         constrained_layout=True)
fig.suptitle('estimated coefficient parameters', fontsize=20)

# 次元番号を初期化
k = 0

# 係数を描画
for r in range(R):
    for c in range(C):

        # 列名を取得
        col_var_label   = col_var_label_lt[k]
        col_param_label = col_param_label_lt[k]

        # 配色の範囲を設定
        beta_max = max(Beta_hat[:, k].max(), abs(Beta_hat[:, k].min()))

        # ラベル用の文字列を作成
        param_label = f'$\\beta_{{i{k}}}$'
        
        # 係数パラメータを描画
        ax = axes[r, c]
        geom_df.plot(ax=ax, column=col_param_label, 
                     #cmap='seismic', vmin=-beta_max, vmax=beta_max, 
                     edgecolor='white', linewidth=0.5, 
                     legend=True, legend_kwds={'label': param_label}) # β_ik
        ax.set_xlabel('$u_j$')
        ax.set_ylabel('$v_j$')
        ax.set_title(col_var_label)
        ax.grid()
        ax.set_aspect('equal', adjustable='box')

        # 次元番号をカウント
        k += 1

plt.show()

　各地点(区域)が基準地点 $i$ に対応し、 $K+1$ 個のグラフはそれぞれ各次元 $k$ の係数パラメータ $\hat{\beta}_{1k}, \cdots, \hat{\beta}_{Nk}$ の値をコロプレス図で示します。

　被説明変数の推定値(理論値)を計算します。

# 被説明変数の推定を計算
Y_hat = X @ Beta_hat.T
print(Y_hat.round(2))
print(Y_hat.shape)

[[-0.43 -0.46 -0.44 ... -0.08 -0.13 -0.44]
 [-0.62 -0.68 -0.65 ...  0.13  0.11 -0.62]
 [-0.3  -0.33 -0.31 ... -0.09 -0.16 -0.33]
 ...
 [-0.13 -0.11 -0.12 ... -0.24 -0.3  -0.13]
 [-0.81 -0.77 -0.8  ... -0.75 -0.72 -0.71]
 [-0.34 -0.35 -0.34 ... -0.15 -0.2  -0.34]]
(159, 159)

　全ての基準地点をまとめた被説明変数の推定値 $\hat{\mathbf{Y}} = (\hat{y}_{11}, \cdots, \hat{y}_{NN})$ を $\hat{\mathbf{Y}} = \mathbf{X} \hat{\boldsymbol{\beta}}$ で計算します。

　被説明変数の推定値のコロプレス図を作成します。

作図コード(クリックで展開)

# 被説明変数の推定値を格納
geom_df['explained'] = np.diag(Y_hat)

# 理論値を作図
fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', 
                       constrained_layout=True)
geom_df.plot(ax=ax, column='explained', 
             edgecolor='white', linewidth=0.5, 
             legend=True, legend_kwds={'label': '$\\hat{y}_{ii}$'}) # 対角要素
ax.set_xlabel('$u_j$')
ax.set_ylabel('$v_j$')
ax.set_title(f'$N = {N}$', loc='left')
fig.suptitle('explained variables', fontsize=20)
ax.grid()
ax.set_aspect('equal', adjustable='box')
plt.show()