Introduction to GWR and MGWR
Setup
# Load libraries
import numpy as np
import pandas as pd
pd.set_option('display.max_columns', None)
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
from pylab import rcParams
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 72
import seaborn as sns
sns.set_style("darkgrid")
sns.set_context(context="paper", font_scale=1.5, rc=None)
sns.set(font="serif")
import plotly.express as px
import plotly.graph_objects as go
import geopandas as gpd
import libpysal as ps
from libpysal import weights
from libpysal.weights import Queen
import esda
from esda.moran import Moran, Moran_Local
import splot
from splot.esda import moran_scatterplot, plot_moran, lisa_cluster, plot_local_autocorrelation
from splot.libpysal import plot_spatial_weights
from giddy.directional import Rose
import statsmodels.api as sm
import statsmodels.formula.api as smf
from stargazer.stargazer import Stargazer, LineLocation
from spreg import OLS
from spreg import MoranRes
from spreg import ML_Lag
from spreg import ML_Error
from mgwr.gwr import GWR, MGWR
from mgwr.sel_bw import Sel_BW
from mgwr.utils import shift_colormap, truncate_colormap
import warnings
warnings.filterwarnings('ignore')
import time
/opt/conda/lib/python3.9/site-packages/esda/getisord.py:636: SyntaxWarning: "is" with a literal. Did you mean "=="?
if __name__ is "__main__":
/opt/conda/lib/python3.9/site-packages/spglm/utils.py:367: SyntaxWarning: "is not" with a literal. Did you mean "!="?
if resetlist is not ():
Import data
# Load Georgia dataset
gdf = gpd.read_file(ps.examples.get_path('G_utm.shp'))
gdf.columns
fig, ax = plt.subplots(figsize=(6, 6))
gdf.plot(color = 'white', edgecolor = 'black', ax = ax)
gdf.centroid.plot(ax=ax)
ax.set_title('Map of Georgia and centroids', fontsize=12)
ax.axis("off")
#plt.savefig('myMap.png',dpi=150, bbox_inches='tight')
plt.show()
fig, ax = plt.subplots(figsize=(6, 6))
gdf.plot(column='PctBach', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=ax)
#gdf.plot(column='PctBach', cmap = 'coolwarm', linewidth=0.01, scheme = 'box_plot', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=ax)
ax.set_title('Population with BA degree (%)', fontsize=12)
ax.axis("off")
#plt.savefig('myMap.png',dpi=150, bbox_inches='tight')
plt.show()
Select variables (y, X)
y = gdf['PctBach'].values.reshape((-1,1)) # reshape is needed to have column array
y.shape
X = gdf[['PctFB', 'PctBlack', 'PctRural']].values
X.shape
Define coordinates (u,v)
u = gdf['X']
v = gdf['Y']
coords = list(zip(u,v))
Estimate GWR model
Select GWR bandwidth
%%time
gwr_selector = Sel_BW(coords, y, X)
gwr_bw = gwr_selector.search()
CPU times: user 2.38 s, sys: 313 ms, total: 2.7 s
Wall time: 3.03 s
print('GWR bandwidth =', gwr_bw)
GWR bandwidth = 117.0
Fit GWR model
gwr_results = GWR(coords, y, X, gwr_bw).fit()
gwr_results.summary()
===========================================================================
Model type Gaussian
Number of observations: 159
Number of covariates: 4
Global Regression Results
---------------------------------------------------------------------------
Residual sum of squares: 2315.466
Log-likelihood: -438.549
AIC: 885.098
AICc: 887.490
BIC: 1529.786
R2: 0.548
Adj. R2: 0.540
Variable Est. SE t(Est/SE) p-value
------------------------------- ---------- ---------- ---------- ----------
X0 14.835 1.400 10.599 0.000
X1 2.110 0.302 6.988 0.000
X2 -0.027 0.018 -1.525 0.127
X3 -0.079 0.014 -5.734 0.000
Geographically Weighted Regression (GWR) Results
---------------------------------------------------------------------------
Spatial kernel: Adaptive bisquare
Bandwidth used: 117.000
Diagnostic information
---------------------------------------------------------------------------
Residual sum of squares: 1650.860
Effective number of parameters (trace(S)): 11.805
Degree of freedom (n - trace(S)): 147.195
Sigma estimate: 3.349
Log-likelihood: -411.653
AIC: 848.915
AICc: 851.350
BIC: 888.212
R2: 0.678
Adjusted R2: 0.652
Adj. alpha (95%): 0.017
Adj. critical t value (95%): 2.414
Summary Statistics For GWR Parameter Estimates
---------------------------------------------------------------------------
Variable Mean STD Min Median Max
-------------------- ---------- ---------- ---------- ---------- ----------
X0 13.653 0.815 11.947 13.636 15.411
X1 2.200 1.078 0.568 2.564 3.415
X2 -0.014 0.027 -0.056 -0.017 0.033
X3 -0.069 0.013 -0.097 -0.065 -0.051
===========================================================================
# As reference, here is the (average) R2, AIC, and AICc
print('Mean R2 =', gwr_results.R2)
print('AIC =', gwr_results.aic)
print('AICc =', gwr_results.aicc)
Mean R2 = 0.6780742669593465
AIC = 848.9154070534352
AICc = 851.3502927844658
# Add R2 to GeoDataframe
gdf['gwr_R2'] = gwr_results.localR2
fig, ax = plt.subplots(figsize=(6, 6))
gdf.plot(column='gwr_R2', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=ax)
ax.set_title('Local R2', fontsize=12)
ax.axis("off")
#plt.savefig('myMap.png',dpi=150, bbox_inches='tight')
plt.show()
Add coefficients to data frame
gdf['gwr_intercept'] = gwr_results.params[:,0]
gdf['gwr_fb'] = gwr_results.params[:,1]
gdf['gwr_aa'] = gwr_results.params[:,2]
gdf['gwr_rural'] = gwr_results.params[:,3]
Filter/correct t-stats
# Filter t-values: standard alpha = 0.05
gwr_filtered_t = gwr_results.filter_tvals(alpha = 0.05)
pd.DataFrame(gwr_filtered_t)
0float64
6.365678624793901 - 8.475430601237301
1float64
0.0 - 8.906962997183667
0
7.575212201
2.043091858
1
7.867274106
0
2
7.721605279
0
3
7.289408655
0
4
6.99319588
8.464019731
5
7.374371638
8.580098061
6
7.520669526
8.684610812
7
6.845622711
8.866210969
8
7.374570257
2.276872322
9
7.877267616
0
# Filter t-values: corrected alpha due to multiple testing
gwr_filtered_tc = gwr_results.filter_tvals()
pd.DataFrame(gwr_filtered_tc)
0float64
6.365678624793901 - 8.475430601237301
1float64
0.0 - 8.906962997183667
0
7.575212201
0
1
7.867274106
0
2
7.721605279
0
3
7.289408655
0
4
6.99319588
8.464019731
5
7.374371638
8.580098061
6
7.520669526
8.684610812
7
6.845622711
8.866210969
8
7.374570257
0
9
7.877267616
0
Map coefficients
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18,6))
gdf.plot(column='gwr_fb', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[0])
gdf.plot(column='gwr_fb', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[1])
gdf[gwr_filtered_t[:,1] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[1])
gdf.plot(column='gwr_fb', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[2])
gdf[gwr_filtered_tc[:,1] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[2])
plt.tight_layout()
axes[0].axis("off")
axes[1].axis("off")
axes[2].axis("off")
axes[0].set_title('(a) GWR: Foreign born (BW: ' + str(gwr_bw) +'), all coeffs', fontsize=12)
axes[1].set_title('(b) GWR: Foreign born (BW: ' + str(gwr_bw) +'), significant coeffs', fontsize=12)
axes[2].set_title('(c) GWR: Foreign born (BW: ' + str(gwr_bw) +'), significant coeffs and corr. p-values', fontsize=12)
plt.show()
Test spatial stationarity
%%time
# Monte Carlo test of spatial variability: 500 iterations
gwr_p_values_stationarity = gwr_results.spatial_variability(gwr_selector, 500)
CPU times: user 6min 36s, sys: 8min 13s, total: 14min 49s
Wall time: 14min 57s
gwr_p_values_stationarity
# Note: The first p-value is for the intercept
Test local multicollinearity
LCC, VIF, CN, VDP = gwr_results.local_collinearity()
pd.DataFrame(VIF)
0float64
1.1414124272949622 - 2.2485786119603466
1float64
1.0065379436521202 - 1.0990689032141991
0
1.182403243
1.021370273
1
1.147263524
1.031686252
2
1.167970589
1.024843495
3
1.279951303
1.066151167
4
1.839177319
1.061995601
5
1.903473887
1.056656555
6
1.978994615
1.048208014
7
1.877960094
1.035297416
8
1.173657434
1.038969892
9
1.158581724
1.037341583
pd.DataFrame(VIF).describe().round(2)
0float64
1float64
count
159
159
mean
1.63
1.05
std
0.36
0.02
min
1.14
1.01
25%
1.26
1.03
50%
1.69
1.05
75%
1.92
1.06
max
2.25
1.1
As the max VIF for each variable is less than 10, there are no signs of local multicollinearity
pd.DataFrame(CN)
0float64
10.068471508458375 - 12.35240178119242
0
10.78576742
1
11.30261189
2
10.98943303
3
10.65286267
4
11.51778843
5
11.75646986
6
11.68556366
7
10.46958525
8
11.3219659
9
11.17220514
gdf['gwr_CN'] = CN
fig, ax = plt.subplots(figsize=(6, 6))
gdf.plot(column='gwr_CN', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=ax)
ax.set_title('Local multicollinearity (CN > 30)?', fontsize=12)
ax.axis("off")
#plt.savefig('myMap.png',dpi=150, bbox_inches='tight')
plt.show()
Estimate MGWR model
Standardize variables
Zy = (y - y.mean(axis=0)) / y.std(axis=0)
ZX = (X - X.mean(axis=0)) / X.std(axis=0)
Select MGWR bandwidths
mgwr_selector = Sel_BW(coords, Zy, ZX, multi=True)
mgwr_bw = mgwr_selector.search()
mgwr_bw
IPyWidgets are not supported
Fit MGWR model
mgwr_results = MGWR(coords, Zy, ZX, mgwr_selector).fit()
mgwr_results.summary()
IPyWidgets are not supported
===========================================================================
Model type Gaussian
Number of observations: 159
Number of covariates: 4
Global Regression Results
---------------------------------------------------------------------------
Residual sum of squares: 71.793
Log-likelihood: -162.399
AIC: 332.798
AICc: 335.191
BIC: -713.887
R2: 0.548
Adj. R2: 0.540
Variable Est. SE t(Est/SE) p-value
------------------------------- ---------- ---------- ---------- ----------
X0 0.000 0.054 0.000 1.000
X1 0.458 0.066 6.988 0.000
X2 -0.084 0.055 -1.525 0.127
X3 -0.374 0.065 -5.734 0.000
Multi-Scale Geographically Weighted Regression (MGWR) Results
---------------------------------------------------------------------------
Spatial kernel: Adaptive bisquare
Criterion for optimal bandwidth: AICc
Score of Change (SOC) type: Smoothing f
Termination criterion for MGWR: 1e-05
MGWR bandwidths
---------------------------------------------------------------------------
Variable Bandwidth ENP_j Adj t-val(95%) Adj alpha(95%)
X0 101.000 3.397 2.466 0.015
X1 101.000 3.512 2.479 0.014
X2 117.000 2.778 2.391 0.018
X3 157.000 1.786 2.218 0.028
Diagnostic information
---------------------------------------------------------------------------
Residual sum of squares: 50.804
Effective number of parameters (trace(S)): 11.474
Degree of freedom (n - trace(S)): 147.526
Sigma estimate: 0.587
Log-likelihood: -134.907
AIC: 294.761
AICc: 297.070
BIC: 333.041
R2 0.680
Adjusted R2 0.655
Summary Statistics For MGWR Parameter Estimates
---------------------------------------------------------------------------
Variable Mean STD Min Median Max
-------------------- ---------- ---------- ---------- ---------- ----------
X0 0.018 0.162 -0.237 0.061 0.251
X1 0.482 0.217 0.119 0.500 0.726
X2 -0.063 0.051 -0.158 -0.057 0.026
X3 -0.301 0.019 -0.345 -0.299 -0.264
===========================================================================
Show bandwidth intervals
mgwr_bw_ci = mgwr_results.get_bws_intervals(mgwr_selector)
print(mgwr_bw_ci)
[(87.0, 104.0), (87.0, 104.0), (114.0, 117.0), (114.0, 157.0)]
# Local R2 are NOT yet implemented in MGWR
#mgwr_results.localR2[0:5]
Add coefficients to data frame
#Add MGWR parameters to GeoDataframe
gdf['mgwr_intercept'] = mgwr_results.params[:,0]
gdf['mgwr_fb'] = mgwr_results.params[:,1]
gdf['mgwr_aa'] = mgwr_results.params[:,2]
gdf['mgwr_rural'] = mgwr_results.params[:,3]
Filter/correct t-stats
# Filter t-values: standard alpha = 0.05
mgwr_filtered_t = mgwr_results.filter_tvals(alpha = 0.05)
# Filter t-values: corrected alpha due to multiple testing
mgwr_filtered_tc = mgwr_results.filter_tvals()
Map coefficients
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18,6))
gdf.plot(column='mgwr_fb', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[0])
gdf.plot(column='mgwr_fb', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[1])
gdf[mgwr_filtered_t[:,1] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[1])
gdf.plot(column='mgwr_fb', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[2])
gdf[mgwr_filtered_tc[:,1] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[2])
plt.tight_layout()
axes[0].axis("off")
axes[1].axis("off")
axes[2].axis("off")
axes[0].set_title('(a) MGWR: Foreign born (BW: ' + str(mgwr_bw[1]) +'), all coeffs', fontsize=12)
axes[1].set_title('(b) MGWR: Foreign born (BW: ' + str(mgwr_bw[1]) +'), significant coeffs', fontsize=12)
axes[2].set_title('(c) MGWR: Foreign born (BW: ' + str(mgwr_bw[1]) +'), significant coeffs and corr. p-values', fontsize=12)
plt.show()
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18,6))
gdf.plot(column='mgwr_intercept', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[0])
gdf.plot(column='mgwr_intercept', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[1])
gdf[mgwr_filtered_t[:,0] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[1])
gdf.plot(column='mgwr_intercept', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[2])
gdf[mgwr_filtered_tc[:,0] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[2])
plt.tight_layout()
axes[0].axis("off")
axes[1].axis("off")
axes[2].axis("off")
axes[0].set_title('(a) MGWR: Intercept (BW: ' + str(mgwr_bw[0]) +'), all coeffs', fontsize=12)
axes[1].set_title('(b) MGWR: Intercept (BW: ' + str(mgwr_bw[0]) +'), significant coeffs', fontsize=12)
axes[2].set_title('(c) MGWR: Intercept (BW: ' + str(mgwr_bw[0]) +'), significant coeffs and corr. p-values', fontsize=12)
plt.show()
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18,6))
gdf.plot(column='mgwr_aa', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[0])
gdf.plot(column='mgwr_aa', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[1])
gdf[mgwr_filtered_t[:,2] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[1])
gdf.plot(column='mgwr_aa', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[2])
gdf[mgwr_filtered_tc[:,2] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[2])
plt.tight_layout()
axes[0].axis("off")
axes[1].axis("off")
axes[2].axis("off")
axes[0].set_title('(a) MGWR: African Ame. (BW: ' + str(mgwr_bw[2]) +'), all coeffs', fontsize=12)
axes[1].set_title('(b) MGWR: African Ame. (BW: ' + str(mgwr_bw[2]) +'), significant coeffs', fontsize=12)
axes[2].set_title('(c) MGWR: African Ame. (BW: ' + str(mgwr_bw[2]) +'), significant coeffs and corr. p-values', fontsize=12)
plt.show()
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18,6))
gdf.plot(column='mgwr_rural', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[0])
gdf.plot(column='mgwr_rural', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[1])
gdf[mgwr_filtered_t[:,3] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[1])
gdf.plot(column='mgwr_rural', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[2])
gdf[mgwr_filtered_tc[:,3] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[2])
plt.tight_layout()
axes[0].axis("off")
axes[1].axis("off")
axes[2].axis("off")
axes[0].set_title('(a) MGWR: Rural (BW: ' + str(mgwr_bw[3]) +'), all coeffs', fontsize=12)
axes[1].set_title('(b) MGWR: Rural (BW: ' + str(mgwr_bw[3]) +'), significant coeffs', fontsize=12)
axes[2].set_title('(c) MGWR: Rural (BW: ' + str(mgwr_bw[3]) +'), significant coeffs and corr. p-values', fontsize=12)
plt.show()
Test spatial stationarity
%%time
# Monte Carlo test of spatial variability: 10 iterations
mgwr_p_values_stationarity = mgwr_results.spatial_variability(mgwr_selector, 10)
IPyWidgets are not supported
IPyWidgets are not supported
IPyWidgets are not supported
IPyWidgets are not supported
IPyWidgets are not supported
IPyWidgets are not supported
IPyWidgets are not supported
IPyWidgets are not supported
IPyWidgets are not supported
IPyWidgets are not supported
CPU times: user 3min 17s, sys: 4min 14s, total: 7min 31s
Wall time: 7min 34s
mgwr_p_values_stationarity
# Note: The first p-value is for the intercept
Test local multicollinearity
mgwrCN, mgwrVDP = mgwr_results.local_collinearity()
gdf['mgwr_CN'] = mgwrCN
fig, ax = plt.subplots(figsize=(6, 6))
gdf.plot(column='mgwr_CN', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=ax)
ax.set_title('Local multicollinearity (CN > 30)?', fontsize=12)
ax.axis("off")
#plt.savefig('myMap.png',dpi=150, bbox_inches='tight')
plt.show()