Chapter 6 Spatial Varying Coefficient Models

6.1 Introduction

This practical extends the regression models that you constructed in Week 4 to the spatial (special?) case, and particularly considers Spatially Varying Coefficient (SVC) models.

A coefficient, \(\beta\), in a regression model describe the relationship between a driver (\(x\), predictor variable) and an outcome (\(y\) the target variable). It describes how \(y\) changes with changes in \(x\). When examined locally, that relationship may vary with location for a number of reasons:

a poor conceptual model
social geography is not physics: laws do not govern behaviours, preferences, attitudes and sentiment
bad measurements
unaccounted for local factors
OR process spatial heterogeneity (that is a process whose statistical relationships vary over space)

In the latter case, ‘whole map’ or global regression models such as the linear regressions in Week 4, may unreasonably assume relationship spatial stationarity. The basic idea behind SVC models is that the relationship between \(x\) and \(y\) may be *spatially non-stationary** (i.e. may vary in different in different locations).

A key advantage of SVCs is that they generate coefficient estimates for different locations in the study area, which can be mapped. The maps indicate how the relationship between an outcome ( \(y\) ) and factors ( \(x\)’s ) varies spatially.

This practical introduces Geographically Weighted Regression (GWR) Fotheringham, Brunsdon, and Charlton (2003) as an approach for handling spatially autocorrelated data. GWR is the SVC ‘brand leader’.

The first part of the practical undertakes a standard GWR and the second explores Multi-scale GWR (MGWR) Fotheringham, Yang, and Kang (2017). MGWR is now considered as the default GWR (Comber et al. 2023). The final optional section illustrates Generalised Additive Models with Gaussian Process splines. The latter are a recent advance (Comber, Harris, and Brunsdon 2024) and are a faster and more precise alternative SVC to the GWR and able to handle a variety of response types (i.e. different distributions of \(y\) ).

6.2 Data and Packages

This session will again use the Georgia data that was used in Week 4 and will again construct a regression model of the median income variable (MedInc), but with slightly different set of predictor variables.

You will need the following packages

library(sf)         # for spatial data
library(tidyverse)  # for everything!
library(GWmodel)    # to undertake the GWR
library(cowplot)    # for plotting multiple graphics

And you will need to load the georgia data as in Week 4 which contains census data for the counties in Georgia, USA.

download.file("https://github.com/gwverse/gw/blob/main/data/georgia.rda?raw=True",
              "./georgia.RData", mode = "wb")
load("georgia.RData")

You should examine this in the usual way - you may want to practice your EDA and spatial EDA outside of the practical - and recall that this data has a number of County level variables including:

% with college degree (PctBach)
% elderly (PctEld)
% foreign born (PctFB)
median income of the county (MedInc in dollars)

6.3 Local vs Global models

Recall that the classic OLS regression has the basic form of:

\[ y_i = \beta_{0} + \sum_{k=1}^{m}\beta_{k}x_{ik} + \epsilon_i \]

where for observations indexed by \(i=1 \dots n\), \(y_i\) is the response variable, \(x_{ik}\) is the value of the \(k^{th}\) predictor variable, \(m\) is the number of predictor variables, \(\beta_0\) is the intercept term, \(\beta_k\) is the regression coefficient for the \(k_{th}\) predictor variable and \(\epsilon_i\) is the random error term that is independently normally distributed with zero mean and common variance \(\sigma^2\). OLS is commonly used for model estimation in linear regression models.

And that this is specified using the function lm (linear model). The model we constructed in Week 5 is recreated here with the response variable \(y\) is median income (MedInc), and the 6 (i.e. \(k\)) socio-economic variables are PctBach, PctEld and PctFB are the different predictor variables, the \(k\) different \(x\)’s:

m = lm(MedInc~ PctBach + PctEld + PctFB, data = georgia) 
summary(m)

## 
## Call:
## lm(formula = MedInc ~ PctBach + PctEld + PctFB, data = georgia)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -27979  -4789   -922   2712  34805 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  51566.1     3528.2  14.616  < 2e-16 ***
## PctBach        798.0      147.9   5.397 2.49e-07 ***
## PctEld       -1789.0      231.2  -7.738 1.21e-12 ***
## PctFB        -1902.0      691.4  -2.751  0.00665 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7672 on 155 degrees of freedom
## Multiple R-squared:  0.4771, Adjusted R-squared:  0.4669 
## F-statistic: 47.13 on 3 and 155 DF,  p-value: < 2.2e-16

Recall also from Week 2 that we are interested in the following regression model properties:

Significance indicated by the p-value in the column \(Pr(>|t|)\) as a measure of how statistically surprising the relationship is;
The Coefficient Estimates (the \(\beta_{k}\) values) that describe how change in the \(x_k\) are related (linearly) to changes in \(y\);
And the Model fit given by the \(R^2\) value, that tells us how well the model predicts \(y\).

All well and good. However the implicit assumptions in linear regression that the relationships between the various \(x\) terms and \(y\) are linear, that the error terms have identical and independent normal distributions and that the relationships between the \(x_k\)’s and \(y\) are the same everywhere, are quite often a consequence of wishful thinking, rather than a rigorously derived conclusion. To take a step closer to reality, it is a good idea to look beyond these assumptions to obtain more reliable predictions, or to gain a better theoretical understanding of the process being considered.

Regression models (and most other models) are aspatial and because they consider all of the data in the study area in one universal model, they are referred to as Global models. For example, the coefficient estimates of the model Median Income in Georgia assume that the contributions made by the different variables in predicting MedInc are the same across the study area. In reality this assumption of spatial invariance is violated in many instances when geographic phenomena are considered.
For example, the association between PctBach and MedInc, might be different in different places. This is a core concept for geographical data analysis, that of process spatial heterogeneity. That is, how and where processes, relationships, and other phenomena vary spatially, as introduced in Week 1 and reinforced in the lecture this week.

The standard OLS regression model (above) can be extended to a Spatially Varying Coefficient (SVC) model. The form is the same except that now location is included and a value for each regression coefficient is assigned to a location \((u,v)\) in space (e.g. \(u\) is an Easting and \(v\) a Northing), so that we denote it \(\beta_k(u,v)\):

\[ y = \beta_{0_(u_i, v_i)} + \sum_{m}^{k=1}\beta_{k}x_{ik_{(u_i, v_i)}} + \epsilon_{(u_i, v_i)} \]

In this session you will explore GWR approaches for solving this equation (and thereby for handling the impacts of any spatial dependencies in the error terms) and for quantifying process spatial heterogeneity using Geographically Weighted Regression

SVC approaches use relative location (i.e. distance between observations) and because degrees have different dimensions depending on where you are on the earth¹, the data need to be converted to a Geographic projection. Here the the georgia data are transformed to a the EPSG 102003 USA Contiguous Albers Equal Area Conic projection.

georgia = st_transform(georgia, "ESRI:102003")

6.4 Geographically Weighted Regression

6.4.1 Description

In overview, geographically weighted (GW) approaches use a moving window or kernel that passes through the study area. At each location being considered, data under the kernel are subsetted, weighted by their distance to the kernel centre and used to make a local calculation of some kind. In GWR this is a regression. In this way GW approaches construct a series of models at discrete locations rather than a single one for whole study area.

There are a number of options for the shape of the kernel - see page 6 in Gollini et al. (2015) - and the default is to use a bi-square kernel. This generates higher weights at locations very near to the kernel centre relative to those towards the edge.

The size or bandwidth of the kernel is critical as it determines the degree of smoothing and the scale over which the GWR operates. This can and should be optimally determined rather than just picking a bandwidth. An optimal bandwidth for GWR can be found by minimising a model fit diagnostic. Options include a leave-one-out cross-validation (CV) score which optimises model prediction accuracy and the Akaike Information Criterion (AIC) which optimises model parsimony - that is it trades off prediction accuracy with model complexity (number of variables used). Note also that the kernel can be defined as a fixed value (i.e. a distance), or in an adaptive way (i.e. with a specific number of nearby data points). Gollini et al. (2015) provide a description of this.

Thus there are 2 basic steps to any GW approach:

Determine the optimal kernel bandwidth (i.e. how many data points to include in each local model or the kernel / window size);
Fit the GW model.

6.4.2 Undertaking a GWR

These steps are undertaken in the code block below with an adaptive bandwidth, noting that georgia needs to be converted to an sp format data (the GWmodel package has not been updated to take sf format objects).

The bandwidth - the kernel size - is critical: it how many data points / observations are included in each local model. However, determining the optimal bandwidth requires different potential bandwidths to be evaluated, which in turn requires some measure of model fit to be applied. There are 2 options for this in a GWR (as indicated above):

leave-one-out Cross Validation (CV). This essentially removes each observation one at a time for a potential bandwidth, and provides a measure of how well the local model with that bandwidth predicts the value of each removed data point;
Akaike’s Information Criterion (AIC or corrected AIC, AICc). This provides a measure of model parsimony, a kind of balance between model accuracy and model complexity. A parsimonious model is one that has a explanatory power but is as simple as possible.

In determining the optimal GWR bandwidth, AIC is preferred to CV because it reflects model parsimony and its use tends to avoid over-fitting GWR models (thus AIC bandwidths tend to be larger than bandwidths found using CV).

The GWmodel package (Lu et al. 2014) has tools for determining the optimal GWR bandwidth, but requires the data to be transformed to sp format from sf:

# convert to sp
georgia.sp = as(georgia, "Spatial")
# determine the kernel bandwidth
bw <- bw.gwr(MedInc~ PctBach + PctEld + PctFB,
             approach = "AIC",
             adaptive = T,
             data=georgia.sp)

You should have a look at the bandwidth: it is an adaptive bandwidth (kernel size) indicating that its size varies but the nearest 48 observations will be used to compute each local weighted regression. Here it, has a value of 48 indicating that 30% of the county data will be used in each local calculation (there are 159 records in the data):

bw

## [1] 48

Try running the below to see the fixed distance bandwidth

# determine the kernel bandwidth
bw2 <- bw.gwr(MedInc~ PctBach + PctEld + PctFB, 
             approach = "AIC",
             adaptive = F,
             data=georgia.sp) 
bw2

This indicates a fixed bandwidth of ~130km. For comparison the distances between each of the counties is summarised below, suggesting that 130km is about 25% of the maximum distance:

summary(as.vector(st_distance(georgia)))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0   85620  156531  166238  237457  517837

You should now create the GWR model with adaptive bandwidth stored in bw, and have a look at the outputs:

# fit the GWR model
m.gwr <- gwr.basic(MedInc~ PctBach + PctEld + PctFB, 
                adaptive = T,
                data = georgia.sp,
                bw = bw)  
m.gwr

##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2024-10-17 16:38:43.867396 
##    Call:
##    gwr.basic(formula = MedInc ~ PctBach + PctEld + PctFB, data = georgia.sp, 
##     bw = bw, adaptive = T)
## 
##    Dependent (y) variable:  MedInc
##    Independent variables:  PctBach PctEld PctFB
##    Number of data points: 159
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##    Min     1Q Median     3Q    Max 
## -27979  -4789   -922   2712  34805 
## 
##    Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
##    (Intercept)  51566.1     3528.2  14.616  < 2e-16 ***
##    PctBach        798.0      147.9   5.397 2.49e-07 ***
##    PctEld       -1789.0      231.2  -7.738 1.21e-12 ***
##    PctFB        -1902.0      691.4  -2.751  0.00665 ** 
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 7672 on 155 degrees of freedom
##    Multiple R-squared: 0.4771
##    Adjusted R-squared: 0.4669 
##    F-statistic: 47.13 on 3 and 155 DF,  p-value: < 2.2e-16 
##    ***Extra Diagnostic information
##    Residual sum of squares: 9123479259
##    Sigma(hat): 7623.079
##    AIC:  3301.791
##    AICc:  3302.183
##    BIC:  3183.48
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: bisquare 
##    Adaptive bandwidth: 48 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: Euclidean distance metric is used.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                  Min.  1st Qu.   Median  3rd Qu.     Max.
##    Intercept 36563.45 51275.92 57838.03 65500.51 87342.20
##    PctBach    -117.33   248.02   493.77   893.53  1937.83
##    PctEld    -3491.73 -2512.30 -2049.18 -1605.08  -793.53
##    PctFB     -8282.57 -3704.10 -1873.14  -408.23  2581.19
##    ************************Diagnostic information*************************
##    Number of data points: 159 
##    Effective number of parameters (2trace(S) - trace(S'S)): 39.54012 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 119.4599 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 3228.433 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 3178.961 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 3145.552 
##    Residual sum of squares: 3695444020 
##    R-square value:  0.7881856 
##    Adjusted R-square value:  0.7174852 
## 
##    ***********************************************************************
##    Program stops at: 2024-10-17 16:38:43.880479

The printout gives quite a lot of information: the standard OLS regression model, a summary of the GWR coefficients and a whole range of fit measures (\(R^2\), AIC, BIC).

6.4.3 Evaluating the results

The code above prints out the whole model summary. You should compare the summary of the local (i.e. spatially varying) GWR coefficient estimates with the global (i.e. fixed) outputs of the ordinary regression coefficient estimates. Note the first line of code below. The output of the GWR model has lots of information including a slot called SDF or spatial data frame and this is an sp format object.

# GWR outputs
names(m.gwr)

## [1] "GW.arguments"  "GW.diagnostic" "lm"            "SDF"          
## [5] "timings"       "this.call"     "Ftests"

Have a look at the SDF (a spatial object in sp format) and notice the use of @data to access the data table in the code below (it contains lots of information). The code below also extracts the coefficients from the OLS model:

head(m.gwr$SDF@data)
# the GWR betas (in the first 4 columns)
summary(m.gwr$SDF@data[, 1:4])
# the OLS coefficients
coef(m)

# get a summary of the GWR coefficients
coefs_gwr = apply(m.gwr$SDF@data[, 1:4], 2, summary)
# link to bandwidths and create a table
tab.gwr = data.frame(t(coefs_gwr), coef(m))
# rename 2nd and 5th columns
names(tab.gwr)[c(2,5)] = c("Q1", "Q3")
# have a look!
tab.gwr

You could write this out to a csv file for inclusion in a report or similar as in the GWR coefficient example above.

write.csv(tab.gwr, file = "gwr_ceofs.csv", row.names = F)

Table 6.1: A comparison of the GWR and OLS regression coefficient estimates.
	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.	Global
Intercept	36563.4	51275.9	57838.0	57684.0	65500.5	87342.2	51566.1
PctBach	-117.3	248.0	493.8	597.3	893.5	1937.8	798.0
PctEld	-3491.7	-2512.3	-2049.2	-2095.0	-1605.1	-793.5	-1789.0
PctFB	-8282.6	-3704.1	-1873.1	-2260.0	-408.2	2581.2	-1902.0

So the Global coefficient estimates are similar to the Mean and Median values of the GWR output. But you can see the variation in the association between median income (MedInc) and the different coefficient estimates.

So for example, for PctPov (poverty) :

the global coefficient estimate is NA.
the local coefficient estimate varies from 3.02072^{4} to 4.31069^{4} between the 1st and 3rd quartiles.

6.4.4 Mapping GWR outputs

The variation in the GWR coefficient estimates can be mapped to show how different factors are varying associated with Median Income in different parts of Georgia. These data are held in the SDF of the GWmodel output which is converted to sf format in the code below and examined:

gwr_sf = st_as_sf(m.gwr$SDF)
gwr_sf

There are lots of information in the GWR output, but here we are interested in are the coefficient estimates and the t-values whose names are the same as the inputs (and the t-values with _TV).

The code below creates maps where large variation in the relationship between the predictor variable and median income was found as in Figure 6.1:

ggplot(gwr_sf) +
  geom_sf(aes(fill = PctEld)) +
  scale_fill_viridis_c() +
  theme_minimal()

Figure 6.1: Maps of the GWR coefficient estimates.

It is also possible to create maps showing coefficients whose values flip from positive to negative in different parts of the state, as in Figure 6.2.

# PctFB
p1 = ggplot(gwr_sf) +
  geom_sf(aes(fill = PctFB)) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_minimal()
# PctBach
p2 = ggplot(gwr_sf) +
  geom_sf(aes(fill = PctBach)) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_minimal()
plot_grid(p1, p2, ncol = 2)

Figure 6.2: Maps of GWR coefficient estimates, whose values flip from positive to negative.

So what we see here are the local variations in the degree to which changes in different variables are associated with changes in Median Income. The maps above are not maps of the predictor variables (PctBach, PctPov and PctBlack). Rather, they show how the relationships between the inputs and Median Income vary spatially - the spatially vary strength of the relationships,

Now a critical aspect is that coefficient significance through the p-value is now local rather than a global statistic. Here we need to go into the t-values from which p-values are derived². Significant coefficient estimates are those whose t-values are less than -1.96 or greater than +1.96. It possible to identify and map these locations from the GWR outputs which have the t-values (_TV) for each input variable.

The code below recreates one of the maps above maps but highlights areas where PctBach is a significant predictor of MedInc (Figure 6.3. It shows that coefficients are highly localized.

# determine which are significant
tval = gwr_sf |>dplyr::select(all_of("PctBach_TV")) |> st_drop_geometry()
signif = tval < -1.96 | tval > 1.96
# map the counties
ggplot(gwr_sf) +
  geom_sf(aes(fill = PctBach), col = NA) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_minimal()  +
  # now add the tvalues layer
  geom_sf(data = gwr_sf[signif,], fill = NA, col = "black")

Figure 6.3: A map of the spatial variation of the PctBach coefficient estimate, with significant areas highlighted.

6.4.5 GWR Summary

GWR investigates if and how relationships between response and predictor variables vary geographically. It is underpinned by the idea that whole map (constant coefficient) regressions such as those estimated by ordinary least squares may make unreasonable assumptions about the stationarity of the regression coefficients under investigation. Such approaches may wrongly assume that regression relationships (between predictor and target variables) are the same throughout the study area.

GWR outputs provides measures of process heterogeneity – geographical variation in data relationships – through the generation of mappable and varying regression coefficients, and associated statistical inference. It has been extensively applied in a wide variety of scientific and socio-scientific disciplines.

There are a number of things to note about GWR:

It reflects a desire to shift away from global, whole map (Openshaw 1996) and one-size-fits-all statistics, towards approaches that capture local process heterogeneity. This was articulated by Goodchild’s (2004) who proposed a 2nd law of geography: the principle of spatial heterogeneity or non-stationarity, noting the lack of a concept of an average place on the Earth’s surface comparable, for example, to the concept of an average human (Goodchild 2004).
There have been a number of criticisms of GWR (see Section 7.6 in Gollini et al. (2015)).
There have been a number of advances to standard GWR notably Multi-Scale which is now considered to be the default GWR. These have all been summarised in a GWR Route map paper (Comber et al. 2023).

Consideration of developments and limitations will be useful for your assignment write up.

Task 1

Create a function for creating a ggplot of coefficients held in an sf polygon object, whose values flip, highlighting significant observations. Note that functions are defined in the following way:

my_function_name = function(parameter1, paramter2) { # note the  curly brackets 
  # code to do stuff
  # values to return
}

A very simple example is below. This simply takes variable name from an sf object, creates a vector of values and returns the median and the mean:

median_mean_from_sf = function(x, var_name) {
  # do stuff
  x = x |> dplyr::select(all_of(var_name)) |> st_drop_geometry() |> unlist()
  mean_x = mean(x, na.rm = T)
  median_x = median(x, na.rm = T)
  # values to return
  return(c(Mean = mean_x, Median = median_x))
}
# test it
median_mean_from_sf(x = georgia, var_name = "MedInc")

##     Mean   Median 
## 37147.09 34307.00

You should be able to use the code that was used to generate the plot above for the do stuff bit:

# determine which are significant
tval = gwr_sf |> dplyr::select(all_of("PctBach_TV")) |> st_drop_geometry()
signif = tval < -1.96 | tval > 1.96
# map the counties
ggplot(gwr_sf) +
  geom_sf(aes(fill = PctBach), col = NA) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_minimal()  +
  # now add the tvalues layer
  geom_sf(data = gwr_sf[signif,], fill = NA, col = "black")

NB You are advised to attempt the tasks after the timetabled practical session. Worked answers are provided in the last section of the worksheet.

6.5 Multi-scale GWR

6.5.1 Description

In a standard GWR only a single bandwidth is determined and applied to each predictor variable. In reality some relationships in spatially heterogeneous processes may operate over larger scales than others>f this is the case the best-on-average scale of relationship non-stationarity determined by a standard GWR may under- or over-estimate the scale of individual predictor-to-response relationships.

To address this, Multi-scale GWR (MGWR) can be used (Yang 2014; Fotheringham, Yang, and Kang 2017; Oshan et al. 2019). It determines a bandwidth for each predictor variable individually, thereby allowing individual response-to-predictor variable relationships to vary. The individual bandwidths provide insight into the scales of relationships between response and predictor variables, allowing for local to global bandwidths. Recent work has suggested that MGWR should be the default GWR (Comber et al. 2023), with a standard GWR only being used in specific circumstances.

6.5.2 Undertaking a MGWR

A MGWR implementation in GWmodel is undertaken in 2 stages: first the bandwidths are found and then they are used to create the final model. Note both stages use the same function, gwr.multiscale but with different arguments passed to it.

The code below undertakes illustrates the search for the optimal bandwidths, and then having found them passes into a second MGWR operation. These will generate a lot of text to your screen!

gw.ms <- gwr.multiscale(MedInc~ PctBach + PctEld + PctFB,
                        data = georgia.sp,
                        adaptive = T, max.iterations = 1000,
                        kernel = "bisquare",
                        bws0=c(100,100,100),
                        verbose = F, predictor.centered=rep(F, 3))

We can examine the bandwidths for the individual predictor variables:

bws = gw.ms$GW.arguments$bws
bws

## [1]  80  53 157  18

These provide indication of the spatial variation (localness) of the scale of each predictor-to-response relationship. Recall that there are 159 counties in Georgia and in that context the adaptive bandwidths printed out above indicate these scales:

# proportions of data points used in each local regression
bws /nrow(georgia.sp)

## [1] 0.5031447 0.3333333 0.9874214 0.1132075

And pass these into a second model.

gw.ms <- gwr.multiscale(MedInc~ PctBach + PctEld + PctFB, 
                        data = georgia.sp,
                        adaptive = T, max.iterations = 1000,
                        kernel = "bisquare",
                        bws0=bws,
                        bw.seled = rep(T, length(bws)),
                        verbose = F, predictor.centered=rep(F, 3))

6.5.3 Evaluating the results

The coefficient estimate spatial variation can be examined (hint: have a look at what is returned by names(gw.ms):

# get a summary of the MGWR coefficients
coefs_msgwr = apply(gw.ms$SDF@data[, 1:4], 2, summary)
# link to bandwidths and create a table
tab.mgwr = data.frame(Bandwidth = bws, 
                 t(round(coefs_msgwr,1)))
# rename 3rd and 6th columns
names(tab.mgwr)[c(3,6)] = c("Q1", "Q3")
# have a look!
tab.mgwr

Table 6.2: The MSGWR variable specific bandwidths and coefficient estimates.
	Bandwidth	Min.	X1st.Qu.	Median	Mean	X3rd.Qu.	Max.
Intercept	80	52701.7	53320.5	54451.5	55197.2	57567.4	58228.9
PctBach	53	563.0	718.9	1011.9	1092.1	1536.1	1601.1
PctEld	157	-2161.5	-2158.7	-2152.8	-2149.8	-2140.2	-2132.9
PctFB	18	-14197.2	-6338.5	-4054.0	-4525.3	-2737.6	482.0

You should examine these, particularly the degree to which the coefficient estimates vary in relation to their bandwidths. Notice that smaller bandwidths (a smaller moving window) results in greater variation and this reduces as the bandwidths tend towards the global (recall that there are \(n = 159\) observations in the georgia data). So some of the bandwidths indicate highly local relationships with the target variable, while other are highly global, with bandwidths approaching \(n = all\).

Notice the coefficients and the degree to which vary, alongside the bandwidths - high bandwidths of course result in less variation in the coefficient estimates.

These tell us something about the local-ness of the interactions of different covariates with the target variable MedInc: PctFB, PctPov and PctBlack have small, highly local bandwidths, using few data points in each local regression. Whilst PctRural, PctBach and PctEld have a global relationship with MedInc and nearly all of their values are used in each local regression model.

In the standard GWR PctFB flipped and it does so here potentially because of the more localised bandwidths (recall that the standard GWR model had a bandwidth of 48).

As before, you could write this out to a csv file for inclusion in a report or similar as in the GWR coefficient example above. And the coefficients can be mapped. To a degree there is no point in mapping the coefficients whose bandwidths tend towards the global (i.e. are close to \(n\)), but the others are potentially interesting.

Again these are held in the SDF of the GWmodel output which is converted to sf format in the code below and examined:

mgwr_sf = st_as_sf(gw.ms$SDF)
mgwr_sf

6.5.4 Mapping MGWR outputs

The code below creates the maps for the local covariates in the same way as before

p1 = 
  ggplot(mgwr_sf) +
  geom_sf(aes(fill = PctBach)) +
  scale_fill_gradient2(low = "red", high = "blue") +  
  theme_minimal()  
# coefficients whose values flip 
p2 = 
  ggplot(mgwr_sf) +
  geom_sf(aes(fill = PctFB)) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_minimal()
plot_grid(p1, p2, ncol = 2)

Figure 6.4: Maps of the MGWR coefficient estimates.

You may have noticed that t-values are also included in the outputs and gain significant areas can be identified and highlighted. This suggests a different spatial extent to that indicated by a standard GWR.

# determine which are significant
tval = mgwr_sf$PctFB_TV
signif = tval < -1.96 | tval > 1.96
ggplot(mgwr_sf) +
  geom_sf(aes(fill = PctFB), col = NA) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_minimal()  +
  # now add the tvalues layer
  geom_sf(data = mgwr_sf[signif,], fill = NA, col = "black")

Figure 6.5: A map of the spatial variation of the PctFB coefficient estimate arising from a MGWR, with significant areas highlighted.

6.5.5 Comparing GWR and MGWR

It is instructive to compare the GWR and MGR outputs, both in terms of the numeric distributions and the spatial ones. The inter-quartile ranges of the numeric distributions give an idea about the different spatial variations and we can examine these with respect to the global coefficient estimates from the OLS model and the MGWR bandwidths:

# OLS
coef(m)

## (Intercept)     PctBach      PctEld       PctFB 
##  51566.0950    798.0094  -1788.9654  -1901.9645

# GWR
apply(m.gwr$SDF@data[,1:4], 2, summary)

##         Intercept   PctBach     PctEld      PctFB
## Min.     36563.45 -117.3257 -3491.7348 -8282.5749
## 1st Qu.  51275.92  248.0238 -2512.3026 -3704.1031
## Median   57838.03  493.7701 -2049.1778 -1873.1398
## Mean     57683.97  597.2625 -2094.9583 -2259.9854
## 3rd Qu.  65500.51  893.5312 -1605.0809  -408.2344
## Max.     87342.20 1937.8310  -793.5303  2581.1860

# MGWR
apply(gw.ms$SDF@data[,1:4], 2, summary)

##         Intercept   PctBach    PctEld       PctFB
## Min.     52701.73  562.9541 -2161.470 -14197.2401
## 1st Qu.  53320.46  718.9348 -2158.735  -6338.5340
## Median   54451.49 1011.8921 -2152.782  -4053.9998
## Mean     55197.23 1092.0693 -2149.825  -4525.3304
## 3rd Qu.  57567.38 1536.1112 -2140.215  -2737.5840
## Max.     58228.88 1601.1373 -2132.857    482.0177

The maps similarly may show large differences in both the value of the coefficient estimates (see Figure 6.6) in different places but also the significant areas (see Figure 6.7).

# gwr
p3 = 
  ggplot(gwr_sf) +
  geom_sf(aes(fill = PctBach), col = NA) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_minimal() + labs(title = "GWR")
# mgwr
p4 = 
  ggplot(mgwr_sf) +
  geom_sf(aes(fill = PctBach), col = NA) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_minimal() + labs(title = "MGWR")
plot_grid(p3, p4, ncol = 2)

Figure 6.6: Maps of the coefficients arising from a GWR and MGWR.

# gwr
tval = gwr_sf$PctFB_TV
signif = tval < -1.96 | tval > 1.96
p5 = 
  ggplot(gwr_sf) +
  geom_sf(aes(fill = PctFB), col = NA) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_map()  + labs(title = "GWR") +
  # now add the tvalues layer
  geom_sf(data = gwr_sf[signif,], fill = NA, col = "black")

# mgwr
tval = mgwr_sf$PctFB_TV
signif = tval < -1.96 | tval > 1.96
p6 = 
  ggplot(mgwr_sf) +
  geom_sf(aes(fill = PctFB), col = NA) +
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_map()  + labs(title = "MGWR") +
  # now add the tvalues layer
  geom_sf(data = mgwr_sf[signif,], fill = NA, col = "black")
plot_grid(p5, p6, ncol = 2)

Figure 6.7: Maps of the signiticant coefficients arising from a GWR and MGWR.

6.5.6 MGWR Summary

A single bandwidth is used in standard GWR under the assumption that the response-to-predictor relationships operate over the same scales for all of the variables contained in the model. This may be unrealistic because some relationships can operate at larger scales and others at smaller ones. A standard GWR will nullify these differences and find a best-on-average scale of relationship non-stationarity (geographical variation), as described through the bandwidth.

MSGWR was proposed to allow for a mix of relationships between predictor and response variables: from highly local to global. In MGWR, the bandwidth for each relationship is determined separately, allowing the scale of individual response-to-predictor relationships to vary.

For these reasons, the MGWR has been recommended as the default GWR. MGWR plus an OLS regression provides solid basis for making a decision to use a GWR and for determining the variation predictor variable to response variable relationships: the MGWR bandwidths indicate their different scales. A standard GWR should only be used in circumstances where all the predictors are found to have the same bandwidth.

6.6 Optional: Generalised Additive Models

6.6.1 Description

Generalised Additive Models (GAMs) (Hastie and Tibshirani 1990) generate multiple model terms (like mini-models) which are added together and to fit non-linear and data with varying relationships and complex interactions (Wood 2017). Consider the relationship between \(x\) and \(y\) in the figure below: it varies in a discontinuous way and the GAM model is able to fit to that.

Figure 6.8: An example of a non-linear relationship between x and y.

To model non-linear relationships this GAMs use something called splines (also called smooths) which can be of different forms depending on the problem and data (Hastie and Tibshirani 1990). A spline can be represented as a linear combination of functions called basis functions. Each of these is assigned a coefficient and these are linearly combined to generate the predictions of \(y\), \(\hat{y}\). If the splines are parameterised by observation location then they can be used to generate spatially varying coefficients (Comber, Harris, and Brunsdon 2024), and if the spline have a Gaussian Process (GP) form then they capture distance observation autocorrelation and decay in the manner suggested by Tobler’s 1st Law of Geography.

Additionally GAMs based SVCs with GP splines address a number of long standing criticisms of GWR based approaches, heat are described in full in Comber et al. (2023) and summarised in Comber, Harris, and Brunsdon (2024):

GWR and MGWR suffer from a lack of generic framework and a clear theoretical model
GWR and MGW generate a collection of local models rather than a single model
MGWR can only handle Gaussian and Poisson responses
MGWR has no options for handling or penalizing collinearity, down-weighting outliers, handling heteroskedastic and autocorrelated error terms.
there are concerns about GWR’s reliability as a spatial predictor

6.6.2 Undertaking a GGP-GAM

The code below extracts the coordinates from the projected georgia data and converts these to kilometres. It also creates a dummy variable for the intercept.

library(mgcv)       # for GAMs
library(stgam)      # for GAM based SVCs
library(tidyverse)  # for everything else
coords = st_coordinates(st_centroid(georgia))

## Warning: st_centroid assumes attributes are constant over geometries

# head(coords)
df.gam = 
  georgia |> 
  mutate(Intercept = 1, 
         X = coords[,"X"]/1000, 
         Y = coords[,"Y"]/1000) |>
  st_drop_geometry() |>
  as_tibble()

First, just for illustration the gam function form the mgcv package is used to generate a standard parametric model. This generates the same OLS model as the that created by lm above and stored in m:

gam_lm = gam(MedInc~ PctBach + PctEld + PctFB, data = df.gam)
summary(gam_lm)

Then a GGP-GAM can be fitted, specifying a GP spline, s for each covariate. Notice the stricture of the spline in the code below:

gam.m = gam(MedInc~ 0 + 
              Intercept + s(X,Y,bs='gp', by=Intercept) + 
              PctBach + s(X,Y,bs='gp', by=PctBach) +
              PctEld + s(X,Y,bs='gp', by=PctEld) +
              PctFB + s(X,Y,bs='gp', by=PctFB),
              data = df.gam)

It is possible to check that the GAM has fitted correctly (if you want to unpick this see the help for the gam.check() function):

gam.check(gam.m)

6.6.3 Evaluating the results

Next we can inspect the model:

summary(gam.m)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## MedInc ~ 0 + Intercept + s(X, Y, bs = "gp", by = Intercept) + 
##     PctBach + s(X, Y, bs = "gp", by = PctBach) + PctEld + s(X, 
##     Y, bs = "gp", by = PctEld) + PctFB + s(X, Y, bs = "gp", by = PctFB)
## 
## Parametric coefficients:
##           Estimate Std. Error t value Pr(>|t|)    
## Intercept  56001.1     3787.8  14.785   <2e-16 ***
## PctBach      453.9     3609.9   0.126    0.900    
## PctEld      -273.6     2853.8  -0.096    0.924    
## PctFB      17806.1   107586.4   0.166    0.869    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                     edf Ref.df     F p-value   
## s(X,Y):Intercept 14.807 17.653 1.185 0.29120   
## s(X,Y):PctBach    8.635 10.657 1.433 0.16957   
## s(X,Y):PctEld     7.653  9.268 0.945 0.49627   
## s(X,Y):PctFB     22.984 25.572 2.369 0.00145 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 132/135
## R-sq.(adj) =  0.805   Deviance explained = 99.1%
## GCV = 3.345e+07  Scale est. = 2.1547e+07  n = 159

This returns summaries of the GGP-GAM fixed (parametric) coefficient estimates which can be interpreted as the linear terms in the model, in a similar way as the results of the standard linear regression model. Here only the Intercept has a significant global effect, but the other covariates do not.

Summaries of the spline smooth terms are also reported. The edf (effective degrees of freedom) summarises the complexity of the spline smooths, with an edf value of 1 indicating a straight line, 2 a quadratic curve etc. Higher edf values indicate increasing non-linearity in the relationship between the covariate and the response. Here, these are for each covariate over a 2-dimensional space defined by \((X,Y)\). The p-values relate to splines / smooths defined over this geographic space and their significance can be interpreted as indicating whether they vary locally over space. In contrast to the parametric coefficients, only the GP spline for PctFB is locally significant (i.e. its relationship with the target variable \(y\) varies locally), whereas the other GP splines are not.

This is a very different model to those suggested by GWR and MGWR: The GAM model suggests a large and significant global Intercept and a significant local effect associated with PctFB.

6.6.4 GGP-GAM SVCs

It also is possible to generate the GGP-GAM SVCs and how the relationships between \(y\) and the \(x\)s vary. The calculate_vcs function from the stgam package below does this for all of the covariates.

# function
terms = c("Intercept", "PctBach", "PctEld", "PctFB")
gam_svc = calculate_vcs(df.gam, gam.m, terms)
# check the results
gam_svc |>
  select(starts_with("b_")) |> 
  sapply(summary) |> 
  t() |> round(1)

##                 Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
## b_Intercept  36238.8 51592.3 56266.4 56001.1 62292.5 70493.2
## b_PctBach     -206.2   164.9   514.7   569.7   942.0  1455.1
## b_PctEld     -2827.7 -2294.6 -1809.4 -1872.6 -1613.1  -463.6
## b_PctFB     -12970.0 -5558.8 -1805.4 -2833.9  -135.3  5261.5

The SVC results and areas of significance can be mapped. The code below assigns the georgia geometry to the results of the SVC and maps the SVC of the significant covariates. Note the assignment of the geometry to the SVC object:

The code below creates the maps for the local covariates:

# assign geometry 
st_geometry(gam_svc) <- st_geometry(georgia)
p7 = ggplot(gam_svc) +
  geom_sf(aes(fill = b_PctFB), col = NA) + 
  scale_fill_gradient2(low = "red", high = "blue") +
  theme_map()
p7

Figure 6.9: Map of the GGP-GAM SVCs for PctFB.

You could compare the GWR, MGWR and GGP-GAM outputs. They indicate subtly different spatial processes.

# GAM
plot_grid(p5, p6, p7 + labs(title = "GAM"), ncol = 3)

6.6.5 GGP-GAM summary

This GAM with GPs generates a very different output the to MGWR, with many more coefficients (try running coef(gam.m)). The GAM SVC is quick and has a clear theoretical framework. However they lack an explicit quantification of the scale of the relationship between target and response variables, as is provided by the GWR and MGWR bandwidths. In contrast, they do however provide a measure of the spatial complexity of the predictor to response relationship and generate more reliable models (compare the \(R^2\) fit metrics of the 3 models).

Some useful descriptions of GAMs can be found at:

And this is a great video if you have a spare hour:

https://www.youtube.com/watch?v=q4_t8jXcQgc

6.7 Answers to Tasks

Task 1

Create a function for creating a ggplot of coefficients held in an sf polygon object, whose values flip, highlighting significant observations. The function below takes three inputs:

x an sf object with the GWR coefficients
var_name the name of a variable in the GWR output
vr_name_TV the name of the t-values for the variable in the GWR output

my_map_signif_coefs_func = function(x, var_name, var_name_TV) {
  # determine which are significant
  tval = x |> dplyr::select(all_of(var_name_TV)) |> st_drop_geometry()
  signif = tval < -1.96 | tval > 1.96
  # map the counties
  p_out = 
    ggplot(x) +
    geom_sf(aes(fill = PctBach), col = NA) +
    scale_fill_gradient2(low = "red", high = "blue") +
    theme_minimal() +
    # now add the tvalues layer
    geom_sf(data = x[signif,], fill = NA, col = "black")
  # return the output
  return(p_out)
}

Test it!

my_map_signif_coefs_func(x = gwr_sf, var_name = "PctFB", var_name_TV = "PctFB_TV")
my_map_signif_coefs_func(x = mgwr_sf, var_name = "PctFB", var_name_TV = "PctBach_TV")

6.8 References

References

Comber, Alexis, Christopher Brunsdon, Martin Charlton, Guanpeng Dong, Richard Harris, Binbin Lu, Yihe Lü, et al. 2023. “A Route Map for Successful Applications of Geographically Weighted Regression.” Geographical Analysis 55 (1): 155–78.

Comber, Alexis, Paul Harris, and Chris Brunsdon. 2024. “Multiscale Spatially Varying Coefficient Modelling Using a Geographical Gaussian Process GAM.” International Journal of Geographical Information Science 38 (1): 27–47.

Fotheringham, A Stewart, Chris Brunsdon, and Martin Charlton. 2003. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. John Wiley & Sons.

Fotheringham, A Stewart, Wenbai Yang, and Wei Kang. 2017. “Multiscale Geographically Weighted Regression (MGWR).” Annals of the American Association of Geographers 107 (6): 1247–65.

Gollini, Isabella, Binbin Lu, Martin Charlton, Christopher Brunsdon, Paul Harris, et al. 2015. “GWmodel: An r Package for Exploring Spatial Heterogeneity Using Geographically Weighted Models.” Journal of Statistical Software 63 (i17).

Goodchild, Michael F. 2004. “The Validity and Usefulness of Laws in Geographic Information Science and Geography.” Annals of the Association of American Geographers 94 (2): 300–303.

Hastie, T, and R Tibshirani. 1990. “Generalized Additive Models. Chapman Hall & CRC.” Monographs on Statistics & Applied Probability. Chapman and Hall/CRC 1.

Lu, Binbin, Paul Harris, Martin Charlton, and Chris Brunsdon. 2014. “The GWmodel r Package: Further Topics for Exploring Spatial Heterogeneity Using Geographically Weighted Models.” Geo-Spatial Information Science 17 (2): 85–101.

Openshaw, Stan. 1996. “Developing GIS-Relevant Zone-Based Spatial Analysis Methods.” Spatial Analysis: Modelling in a GIS Environment, 55–73.

Oshan, Taylor M, Ziqi Li, Wei Kang, Levi J Wolf, and A Stewart Fotheringham. 2019. “Mgwr: A Python Implementation of Multiscale Geographically Weighted Regression for Investigating Process Spatial Heterogeneity and Scale.” ISPRS International Journal of Geo-Information 8 (6): 269.

Wood, Simon N. 2017. Generalized Additive Models: An Introduction with r. CRC press.

Yang, Wenbai. 2014. “An Extension of Geographically Weighted Regression with Flexible Bandwidths.” PhD thesis, University of St Andrews.

see https://www.thoughtco.com/degree-of-latitude-and-longitude-distance-4070616 ↩︎
There is a nice description here at the start of this webpage - ignore all the MiniTab stuff! https://blog.minitab.com/blog/statistics-and-quality-data-analysis/what-are-t-values-and-p-values-in-statistics ↩︎