Section 12.2 Geographically Weighted Regression
The basic idea behind Geographically Weighted Regression (GWR) is to explore how the relationship between a dependent variable (\(Y\)) and one or more independent variables (\(X\)s) might vary geographically. Instead of assuming that a single model can be fitted to the entire study region, it looks for geographical differences in the relationship. This is achieved by fitting a regression equation to every feature in the dataset. We construct these separate equations by incorporating the dependent and explanatory variables of the features falling within the neighbourhood of each target feature.
GWR can be used for a variety of applications, including the following:
-
Is the relationship between homicide rate and resource deprivation consistent across the study area?
-
What are the key variables that explain high homicide rates?
-
Where are the counties in which we see high homicide rates? What characteristics seem to be associated with these? Where is each characteristic most important?
-
Are the factors influencing higher homicide rates consistent across the study area?
GWR provides three types of regression models: Continuous, Binary, and Count. These types of regression are known in statistical literature as Gaussian, Logistic, and Poisson, respectively. The Model Type for your analysis should be chosen based on the level of measurement of the Dependent Variable \(Y\text{.}\)
So to illustrate, let us first build our linear model once again, regressing homicide rate in the 1970s for each county against resource deprivation (RD70), population structure (PS70), median age (MA70), divorce rate (DV70), and unemployment rate (UE70).
# make linear model
model <- lm(HR70 ~ RD70 + PS70 + MA70 + DV70 + UE70, data = ncovr)
summary(model)
## ## Call: ## lm(formula = HR70 ~ RD70 + PS70 + MA70 + DV70 + UE70, data = ncovr) ## ## Residuals: ## Min 1Q Median 3Q Max ## -21.12 -3.14 -1.02 1.90 69.75 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 11.1376 0.7602 14.65 <2e-16 *** ## RD70 4.1541 0.1147 36.23 <2e-16 *** ## PS70 1.1449 0.1165 9.82 <2e-16 *** ## MA70 -0.2080 0.0233 -8.94 <2e-16 *** ## DV70 1.4066 0.1064 13.22 <2e-16 *** ## UE70 -0.4268 0.0497 -8.59 <2e-16 *** ## --- ## Signif. codes: ## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 6.03 on 3079 degrees of freedom ## Multiple R-squared: 0.328, Adjusted R-squared: 0.327 ## F-statistic: 301 on 5 and 3079 DF, p-value: <2e-16
And now, to get an idea of how our model might perform differently in different areas, we can map the residuals across each county:
# map residuals
ncovr$resids <- residuals(model)
qtm(ncovr, fill = "resids")
If we see that there is a geographic pattern in the residuals, it is possible that an unobserved variable may be influencing the dependent variable, and there might be some sort of spatial variation in our model’s performance. This is what GWR allows us to explore, whether the relationship between our variables is stable over space.
Subsection 12.2.1 Calculate kernel bandwidth
The big assumption on which the GWR model is based is the kernel bandwitdh. As mentioned earlier, in order to calculate local coefficients for our study regions, GWR takes into account the data points within our specified bandwith, weighting them appropriately. The bandwidth we choose will determine our results. If we choose something too large, we might mask variation and get results similar to the global results. If we choose something too small, we might get spikes, and lots of small scale variation, but also possibly large standard errors and unreliable estimates. Therefore, the process of choosing an appropriate bandwidth is paramount to GWR.
It is possible that we choose bandwidth manually, based on some prior knowledge, or a theoretically based argument. However, if we do not have such ideas to start, we can use data-driven processes to inform our selection. The most frequently used approach is to use cross-validation. Other approaches include Akaike Information Criteria, and maximum likelihood.
Here we use
gwr.sel() which uses a cross-validation method for selecting an optimal method for GWR. The function finds a bandwidth for a given geographically weighted regression by optimizing a selected function. For cross-validation, this scores the root mean square prediction error for the geographically weighted regressions, choosing the bandwidth minimizing this quantity.
With the
adapt parameter you can choose whether you want the bandwidth to adapt (i.e., find the proportion between 0 and 1 of observations to include in weighting scheme (k-nearest neighbours)), or whether you want to find a global bandwidth. A global bandwidth might make sense where the areas of interest are of equal size, and equally spaced; however, if we have variation, for example like in the case of the counties of the United States, it makes sense to adjust the bandwidth for each observation. Here we set to adapt.
Another important consideration is the function for the geographical weighting. We mentioned this above. In the
gwr.sel() function you can set this with the gweight= parameter. By default, this is set to a Gaussian function. We will keep it that way here.
With an sf object (which our data are in), one of the parameters required is a coordinate point for each one of our observations. In this case, what we can do is grab the centroid of each one of our polygons. To do this, we can use the
st_coordinates() function:
ncovr <- ncovr %>%
mutate(cent_coords = st_coordinates(st_centroid(.)))
Now that we have our centroids, we can calculate our optimal bandwith, and save this into an object called
gwr_bandwidth.
gwr_bandwidth <- gwr.sel(HR70 ~ RD70 + PS70 + MA70 + DV70 + UE70,
data = ncovr,
coords = ncovr$cent_coords,
adapt = T)
## Adaptive q: 0.382 CV score: 109039 ## Adaptive q: 0.618 CV score: 110549 ## Adaptive q: 0.2361 CV score: 107560 ## Adaptive q: 0.1459 CV score: 105925 ## Adaptive q: 0.09017 CV score: 104176 ## Adaptive q: 0.05573 CV score: 102283 ## Adaptive q: 0.03444 CV score: 100199 ## Adaptive q: 0.02129 CV score: 98707 ## Adaptive q: 0.01316 CV score: 98388 ## Adaptive q: 0.01161 CV score: 98447 ## Adaptive q: 0.01475 CV score: 98373 ## Adaptive q: 0.0145 CV score: 98390 ## Adaptive q: 0.01725 CV score: 98566 ## Adaptive q: 0.01571 CV score: 98411 ## Adaptive q: 0.01512 CV score: 98356 ## Adaptive q: 0.01534 CV score: 98364 ## Adaptive q: 0.01508 CV score: 98358 ## Adaptive q: 0.01517 CV score: 98354 ## Adaptive q: 0.01524 CV score: 98351 ## Adaptive q: 0.01528 CV score: 98356 ## Adaptive q: 0.01524 CV score: 98351
So we now have our bandwidth. It is:
gwr_bandwidth
## [1] 0.01524
As our data are in WGS84, this value is in degrees.
Subsection 12.2.2 Building the geographically weighted model
Now we can build our model. The
gwr() function implements the basic geographically weighted regression approach to exploring spatial non-stationarity for given bandwidth and chosen weighting scheme. We specify the formula and our data source, as with the linear and spatial regressions, and the coordinates of the centroids, as we did in the bandwidth calculation above; but now we also specify our bandwith, which we created earlier, but we do this with the adapt= parameter. This is because, by default, the function takes a global bandwidth with the bandwidth= parameter, and we have created an adaptive bandwidth, so we shall use this adapt= parameter to include our GWRbandwidth object. We also set the hatmatrix = parameter to TRUE, to return the hatmatrix as a component of the result, and the se.fit= parameter to TRUE to return local coefficient standard errors - as we set hatmatrix to TRUE, two effective degrees of freedom sigmas will be used to generate alternative coefficient standard errors.
gwr_model = gwr(
HR70 ~ RD70 + PS70 + MA70 + DV70 + UE70,
data = ncovr,
coords = ncovr$cent_coords,
adapt = gwr_bandwidth,
hatmatrix = TRUE,
se.fit = TRUE
)
gwr_model
## Call: ## gwr(formula = HR70 ~ RD70 + PS70 + MA70 + DV70 + UE70, data = ncovr, ## coords = ncovr$cent_coords, adapt = gwr_bandwidth, hatmatrix = TRUE, ## se.fit = TRUE) ## Kernel function: gwr.Gauss ## Adaptive quantile: 0.01524 (about 47 of 3085 data points) ## Summary of GWR coefficient estimates at data points: ## Min. 1st Qu. Median 3rd Qu. Max. ## X.Intercept. -4.2897 4.4361 7.4395 12.3212 30.1335 ## RD70 -0.3625 2.0168 2.8788 3.8206 6.6405 ## PS70 -4.0630 0.2914 1.0987 1.9019 3.4998 ## MA70 -0.6831 -0.2072 -0.1028 -0.0186 0.2971 ## DV70 -0.7513 0.5368 0.9652 1.4351 3.1027 ## UE70 -1.4895 -0.4286 -0.1787 -0.0220 1.8486 ## Global ## X.Intercept. 11.14 ## RD70 4.15 ## PS70 1.14 ## MA70 -0.21 ## DV70 1.41 ## UE70 -0.43 ## Number of data points: 3085 ## Effective number of parameters (residual: 2traceS - traceS'S): 244 ## Effective degrees of freedom (residual: 2traceS - traceS'S): 2841 ## Sigma (residual: 2traceS - traceS'S): 5.405 ## Effective number of parameters (model: traceS): 171.2 ## Effective degrees of freedom (model: traceS): 2914 ## Sigma (model: traceS): 5.337 ## Sigma (ML): 5.186 ## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 19276 ## AIC (GWR p. 96, eq. 4.22): 19082 ## Residual sum of squares: 82986 ## Quasi-global R2: 0.5022
This output tells us the 5-number summary for the coefficients (\(Y\)) for each predictor variable \(X\text{.}\) We can see for all predictor variables, some counties have negative values, while other counties have positive values. Evidently there is some variation locally in the relationships between these independent variables and our dependent variable.
For this to be informative, however, we should map this variation - this is where we can gain insight into how the relationships may change over our study area. In our
gwr_model object, we have a Spatial Data Frame. This is a SpatialPointsDataFrame or SpatialPolygonsDataFrame object which contains fit.points, weights, GWR coefficient estimates, R square, and coefficient standard errors in its "data" slot. Let’s extract this into another object called gwr_results:
gwr_results <- as.data.frame(gwr_model$SDF)
To map this, we can join these results back to our original data frame.
gwr_results <- cbind(ncovr, gwr_results)
Now we can begin to map our results. First, let’s see how well the model performs across the different counties of the United States. To do this, we can map the
localR2 value which contains the local \(R^2\) value for each observation.
qtm(gwr_results, fill = "localR2" )
## old-style crs object detected; please recreate object with a recent sf::st_crs() ## old-style crs object detected; please recreate object with a recent sf::st_crs() ## old-style crs object detected; please recreate object with a recent sf::st_crs() ## old-style crs object detected; please recreate object with a recent sf::st_crs() ## old-style crs object detected; please recreate object with a recent sf::st_crs()
We can see some really neat spatial patterns in how well our model performs in various regions of the United States. Our model performs really well in the North East for example, but not so well in the South.
We can also map the coefficients for specific variables. Let’s say we’re interested in how the relationship between resource deprivation (RD70) and homicide rate.
NOTE: In ourgwr_resultsdata frame we have a variable calledRD70; this contains the observed resource deprivation value for the county. We then haveRD70.1which is the coefficient value for each county. In the results contained in the SDF element of our model (gwr_model$SDF) it was namedRD70but when we joined the two objects, the.1was appended to allow us to differentiate between the two.
qtm(gwr_results, fill = "RD70.1" )
We can see in this map again some regional patterns in the coefficient of resource deprivation. The relationship is stronger in the counties in the Northeast, while in other areas it is weak or even 0. And in some counties in Texas, we can see that the direction of the relationship flips. In these areas the coefficients are negative. This means that resource deprivation is associated with lower homicide rates in these counties.
This may be interesting and may raise questions; however, it is important to consider not only the coefficients, but also the standard errors associated with them, so we can get some ideas behind the reliability of these estimates. So let us have a look by mapping the standard errors for these coefficients, which are contained in the
RD70_se variable.
qtm(gwr_results, fill = "RD70_se" )
Higher standard errors mean that our estimates may be less reliable. So it is important to keep in mind both maps: the one with the coefficients, and the one showing the standard errors, when drawing conclusions from our GWR results.
Subsection 12.2.3 Using GWR
Overall, GWR is a good exploratory tool which can help illustrate how a particular model might apply differently in different regions of your study area, or how relationships between variables may vary spatially. This technique can be extended, for example to count models, or logit models, but (following the lead of Luc Anselin) we stop here, and pause to highlight some issues. Specifically, our results are very sensitive to our choice of bandwidth. The results we achieve with this approach might not show up when we use other bandwiths. Additionally, GWR does not really explain anything. What we do is demonstrate that the relationship between our variables is not stable over space. But we cannot with this technique explain why it is not stable over space.
Nevertheless, it is a great technique to identify spatial variation, which, unlike spatial regimes, does not rely on our a priory segmentation of our study region, but instead it allows for the variation to emerge from the data themselves.
