Encapsulating Spatially Varying Relationships with a Generalized Additive Model

1. Introduction

In regression, the relationship between response and predictor variables may vary with observation location. This may be due to a poor conceptual model, poor measurements with regional biases, omitted local factors, or it may be due to the presence of process spatial heterogeneity, also known as process spatial non-stationarity. In such cases, spatially varying coefficient (SVC) models provide an alternative to global or whole map regressions in which the relationships between response and predictor variables are assumed to be constant across space. SVCs relax the assumption of relationship spatial stationarity [1]. SVC models allow relationships to vary with observation location and generate spatially distributed or local coefficient estimates which can be mapped to indicate how and where statistical relationships vary. They provide a methods for quantifying process spatial non-stationarity or spatial heterogeneity - a key and increasingly common task in spatial data analysis. Importantly SVCs can provide insights by informing on the spatially varying nature of drivers of a response, or may guide further investigations.

The SVC brand leader is Geographically Weighted Regression (GWR) [2]. This uses a single moving window or kernel, whose size is optimised and determines the scale of spatial heterogeneity in the resultant outputs. GWR has been extended to a multiscale approach (MGWR) [3,4] in which multiple kernels are fitted, one for each predictor variable, thereby capturing multiple scales of spatially varying relationship. MGWR is considered by some to be the default GWR model [5]. However, there are a number of limitations associated with MGWR approaches, the main ones being that they cannot predict out of sample and that only approaches for Gaussian and Poisson responses have been developed [6].

This paper addresses these critical gaps and describes the construction of a novel SVC modelling approach that uses Generalized Additive Models (GAMs) [7,8] with Gaussian Process (GP) smooths to capture process spatial non-stationarity. GAM smooths are commonly used to accommodate varying relationships between predictor and response variables in attribute space. But if they are parameterised with observation location, then the result is an SVC as the coefficient estimates vary over geographic space. The GAM-based SVC approach proposed in this paper - the Geographical Gaussian Process GAM (GGP-GAM) - is evaluated through two simulation experiments. The first compares it with MGWR using simulated data with complex and highly localised degrees of spatial heterogeneity. Both models are fully tuned. The second investigates the ability of the GGP-GAM to handle larger datasets. Specifically, it examines whether performance is affected as the number of knots required by the GAM increases and as the number of observations and regression terms increases. Finally, an empirical case study of the UK referendum to leave the European Union in 2016 (Brexit) is used illustrate and compare MGWR and GAM-based approaches.

2. Literature Review

2.1. SVC Models

Spatially varying coefficient regression models seek to capture spatial dependencies in the relationships between target and predictor variables. Their outputs provide explicit indications of the scales of any spatial dependencies and the spatial variations in the predictor-to-response relationships. A widely used SVC model is Geographically Weighted Regression (GWR) [2]. GWR has been extended to Multiscale GWR (MGWR) [3,4] in which the spatial dependency between each predictor variable and the target variable is individually quantified. In a MGWR / GWR a series of local models are created using data subsets extracted using a moving window or kernel. For each local model the subsetted data are weighted by their distance to the location under consideration, resulting local, spatially distributed coefficient estimates. The kernel size is determined using a weighted least squares algorithm in a standard GWR and an iterative back-fitting algorithm in MGWR. Besides GWR/MGWR, other widely utilized multiscale SVC models include Bayesian Gaussian Process (GP) models that employ co-kriging via a Linear Model of Co-regionalisation (LMC) (Bayes-GP) [9,10,11,12], and Eigenvector Spatial Filtering (ESF) with Moran coefficients (ESF-MC) [13,14,15]. The ESF-MC enhances the deterministic ESF model [16] by incorporating random effects to capture stochastic spatial processes. These SVC models are directly interrelated [15], and comparative analyses indicate that no single multiscale SVC model significantly surpasses the others. Recently, other multiscale SVC models have emerged, such as a triangulation model utilizing bivariate spline estimators [17], a GP-based model adopting a frequentist approach with ML estimation [18], and a GLM-based model incorporating a reduced-rank spline [19]. However, these SVC models have not been widely adopted by the research community and the use of GWR approaches has increased rapidly [20]. There are a number of theoretical limitation with GWR approaches (including MGWR): they generate a collection of local models rather than a single one while Bayes-GP, ESF-MC, and GGP-GAM each offer a single non-stationary model formulation (e.g. [10,21]), they use the same observations multiply in different local models and MGWR approaches cannot predict out of sample. As a result it has been argued that GWR / MGWR are essentially exploratory approaches [18,22].

2.2. GAM Based SVC Models

GAMs produce multiple model terms that are additively combined, and can accommodate different types of response variables [7,8,23]. In this way they intuitively model and capture complex data relationships and interactions, including non-linearities between predictor and response variables [24]. GAMs provide a framework for making predictions from complex systems, for quantifying variable relationships, and for making inferences about these relationships [7,8]. They perform as well as many machine learning models in terms of prediction accuracy and computational speed [25], and critically, in the context of varying coefficient modeling, GAMs combine this inferential power with transparency and process understanding. They have been described as providing "intrinsically understandable white-box machine learning models that provide a technically equivalent, but ethically more acceptable alternative to [machine learning] black-box models" ([26] p.2).

GAMs can incorporate smooths or splines to model nonlinear relationships, whose forms can vary depending on the problem and / or the data [8]. Smooths are composed of combinations of basis functions which can be single or multi-dimensional in terms of predictor variables. Consequently, a GAM with smooths is composed of sums of multi-dimensional basis functions that enable complex relationships to be modelled. As a result of these properties a number varying coefficient and Bayesian regression models using GAMs have been proposed [27,28,29,30]. These make few assumptions about the response variable distribution, which may be skewed, kurtotic or discrete, and their varying coefficient aspects relate to target variable properties of location, scale and shape. In these ’location’ refers to target variable central tendencies, ’scale’ to measures of spread such as standard deviation in normal distributions or parameters of a gamma distribution, and ’shape’ refers to the nature of the skew and kurtosis present in the distribution. In this way a typical GAM-based varying coefficient model approach can accommodate spatial effects, but not within in an explicitly geographic framework way with respect to location. Some recent work has proposed an SVC GLM regression model using GAMs [19] with reduced-rank thin-plate smooths, however an alternative approach for modelling predictor-to-response relationship, is to consider as a Gaussian Process (GP) over space. Rather than use a thin-plate approach to spatially varying relationships, a GP captures relationships that decay over distance, reflecting Tobler’s First Law of Geography [31]. As GPs have been found to be effective in handling spatial autocorrelation [32], the approach proposed in this paper uses decomposed or low-rank GP smooths parameterised with observation location for each covariate. It extends initial investigations of GAMs in SVC modelling using the Geographical Gaussian Process GAM (GGP-GAM) [33] by exploring GAM tuning with more complex spatial problems.

3. Methodology

The analyses described below first compared SVC models using GGP-GAM and MGWR approaches through with analyses of 100 simulated each with 625 observations using 3 predictor variables (i.e. $x_1$, $x_2$, $x_3$). The accuracy in estimating the true coefficient surfaces and in predicting the true values of y was investigated. A second tranche of analyses investigated the ability of GAM based SVC models to handle larger datasets in this case with 250,000 observations) and specifically to investigate model tuning via the knots parameter, k in order to identify possible trade-offs between model performance and computation overheads. This section describes the data, the formulation of the Geographical Gaussian Process GAM (GGP-GAM), the analyses comparing GGP-GAM with GWR, an investigation of GGP-GAM tuning with a larger dataset and an empirical case study.

3.1. Data

For the first analyses to compare the proposed GGP-GAM and MGWR SVC models, a simulated spatially varying coefficient dataset was created, with three coefficients each having different degrees of spatial heterogeneity. The simulated coefficient dataset was created by extracting Moran eigenvectors [13] generated using the spmoran R package [34]. Here the 1st, 10th and 25th surfaces were selected from the matrix of the first L eigenvectors, for a 25 by 25 grid (625 observations). Each coefficient surface, termed ${\beta }_1$, ${\beta }_2$, ${\beta }_3$ has a normal distribution and was rescaled to have a mean of 0 and standard deviation of 1. An intercept surface, ${\beta }_0$ with a constant value of 2 was also defined. These represent the true regression coefficients and are shown in Figure 1.

The values for 100 sets of simulated predictors $x_1$, $x_2$ and $x_3$ were calculated from a random normal distribution and were rescaled to have a range of [0, 1] and the values for the error term $ϵ$ from a random normal distribution rescaled to range have a of [0, 0.25]. One hundred simulated datasets were created. The simulated response y was calculated directly from these and the coefficients in Figure 1.

For the second analyses, simulated data were again created from Moran eigenvectors [13] but now over a 500 by 500 grid and the 1st, 10th and 25th surfaces were extracted to create 3 true coefficients, over 250,000 observations. The coefficients have different spatial heterogeneities to those in Figure 1. Values for $x_1$, $x_2$, $x_3$, $ϵ$ and y were generated in the same way as for the smaller simulated study. The true coefficients are shown in Figure 2. Note that the range of the true coefficients is slightly larger than those in Figure 1 and their spatial variations are different.

3.2. GAM-Based SVC Models

The standard form for an OLS regression is:

$$y_i={\beta }_0+\sum _{j=1}^m{\beta }_j{x}_{ij}+{ϵ}_i$$

(1)

where for observations indexed by $i=1\cdots n$, $y_i$ is the response variable, $x_{ij}$ is the value of the j-th predictor variable, m is the number of predictor variables, ${\beta }_0$ is the intercept term, ${\beta }_j$ is the regression coefficient for the j-th predictor variable and ${ϵ}_i$ is the random error term. This can be extended to define a SVC regression model:

$$y_i={\beta }_{0_{(u_i,v_i)}}+\sum _{j=1}^m{\beta }_{j_{(u_i,v_i)}}{x}_{ij}+{ϵ}_i$$

(2)

where now $(u_i,v_i)$ are the spatial coordinates of the observations i and ${\beta }_j(u_i,v_i)$ are the coefficients estimated at those locations.

GAMs can also calibrate regression models in which the functions of the predictor variables are unknown and these take the following form:

$$y=\alpha +f_1\left(z_1\right)+f_2\left(z_2\right)+\cdots +f_j\left(z_j\right)+ϵ$$

(3)

where $\alpha$ is the intercept (${\beta }_0$ in Equations 1 and 2) and the j-th predictor $z_j$ can be a scalar or a vector.

These can be extended so that each$f_j\left(z_j\right)$ becomes a linear regression coefficient on another scalar predictor is a linear regression $x_j$:

$$y=\alpha +f_0\left(z_0\right)+x_1{f}_1\left(z_1\right)+x_2{f}_2\left(z_2\right)+\cdots +x_m{f}_m\left(z_m\right)+ϵ$$

(4)

And, if $z_0=z_1=\cdots z_m=z$ for example, and z is a vector of locations or coordinates (similar to $(u,v)$ in Equations 1 and 2) then a spatially varying coefficient model is specified:

$$y=\alpha +f_0\left(z\right)+x_1{f}_1\left(z\right)+x_2{f}_2\left(z\right)+\cdots +x_m{f}_m\left(z\right)+ϵ$$

(5)

where $f_j\left(z\right)$ represents the j-th SVC.

One approach to specify $\alpha \left(z\right)\cdots f_m\left(z\right)$ is to generate each function from a GP smooth defined as follows:

$$f_j\left(z\right)=\sum _{p=1}^P{\kappa }_p\left(z\right){\gamma }_{j,p}{\gamma }_{j,p}\sim N(0,s_j^2)$$

(6)

where ${\kappa }_p\left(z\right)$ represents a basis function, ${\gamma }_{j,p}$ is the corresponding coefficient and $s_j^2$ is the smoothing parameter. The SVC $f_j\left(z\right)$ exhibits high spatial variation if $s_j^2$ is large, and remains constant if $s_j^2=0$. In this way the smoothing parameter determines the wiggliness of the SVC, that is how quickly the function varies across its range.

The basis function ${\kappa }_p\left(z\right)$ is specified such that $f_j\left(z\right)$ has a spatially smooth pattern, in this case driven by the following covariance function:

$$\mathrm{Cov}(f_m\left(z\right),f_j\left(z^{\prime }\right)=s_j^{2}c\left(d_{z,z^{\prime }}\right)$$

(7)

where $d_{z,z^{\prime }}$ represents the distance between locations z and $z^{\prime }$ and the covariance function $c\left(d_{z,z^{\prime }}\right)$ decreases with increases in $d_{z,z^{\prime }}$. This is similar to kriging and MGWR. Kriging uses covariance functions via the semi-variogram and a covariance function is calibrated individually for each $f_j\left(z\right)$ in the same way as bandwidths are optimised in MGWR. Thus a key task in using GAMs to model SVCs is to estimate the smoothness parameters $s_j^2$ and thereby $f_j\left(z\right)$.

A GAM employs smooth functions of the predictor variables, assuming that the values of the explained variable y follow a some kind of exponential distribution:

$$f\left(y\right|\theta )=h\left(y\right)exp\left(\nu \left(\theta \right)T\left(x\right)-A\left(\theta \right)\right)$$

(8)

where $h(.)$, $\nu (.)$, $T(.)$ and $A(.)$ are known functions, and $\theta$ is a vector of parameters. This has a flexible form and is thus able to handle a wide range of distributions including Gaussian, Poisson, Gamma, or Binomial ones.

In summary, GAMs can estimate SVCs by modelling non-linear relationships as spatially non-stationary ones using GP splines over geographical as well as attribute space: a GP spline with spatial location is defined for each predictor variable j. Additionally an offset ${\alpha }_j$ is added to account for the GP’s zero mean, thereby creating the SVCs ${\beta }_j\left(z\right)$. In the standard 2D case, $z_i=(u_i,v_i)$, defining the Geographical Gaussian Process GAM (GGP-GAM).

3.3. Analysis I: Comparing GGP-GAM and MGWR

Gaussian Process splines can be incorporated within GAM models and here this was done using the gam function in the mgcv R package [35]. The observation spatial locations z were extracted and used to parameterise the splines with a GP smooth. The GPs constructed in this manner have a mean of zero; therefore, an additional fixed offset term is included for each covariate alongside the spatially smoothed terms.

There are some important tuning considerations in fitting GAMs with mgcv including the specification of spline smoothing parameters and the number of knots. These parameters influence the nature of the varying coefficients and the characterization of response-to-predictor relationships across space. The smoothing parameter regulates the degree of data smoothing through a correlation function, and the GAM function in the mgcv package optimizes this using an estimation method. In this case, REML was specified for optimization instead of a cross-validation estimator. The number of knots, k, determines the basis dimensions, the maximum number of base functions, used in the smooth and thus the degree of sensitivity in fitting the model to the data. For example, consider the models predicting y specified with different numbers of knots in Figure 3. Together with the smoothing parameter, determines the degree of relationship non-linearity in the smooth. These are heuristically optimised by most GAM implementations in order to balance over-fitting with capturing relationship complexity. In the approach suggested by Wood [35], k can be manually and iteratively determined: an initial value is set and the resultant splines are examined to determine whether the smoothing optimization has converged and whether the basis dimensions defined by k are adequate. If they not, they are increased. For the analyses comparing GGP-GAM and MGWR, $k=155$ was specified after investigation.

The GWmodel R package [36,37] was used to construct the MGWR models and were specified with an adaptive bisquare kernel. Note that fixed kernels were also investigated and the comparisons with the GGP-GAM results reported below were found to be similar. MGWR can accommodate and quantify varying degrees of spatial heterogeneity in the predictor-to-response relationships by determining a kernel bandwidth for each covariate. This requires a series of optimisation problems to be solved and the implementation of MGWR uses a back-fitting algorithm to do this, as does all implementations of MGWR [37,38,39]. Here, the convergence of the back-fitting procedure was evaluated through cross validation of the residual sum of squares and the threshold for convergence was set at ${10}^{-5}$ or 2000 convergence iterations.

In the first analyses, the GGP-GAM and MGWR models were fitted to the 100 simulations over the 625 surface points. The SVC estimates were extracted for each model, along with predictions for y (i.e. $\widehat{y}$) for each simulation. These were used to calculate measures of coefficient fit by comparing the model coefficient estimates with the true ones in Figure 1, including R² for ${\beta }_1$, ${\beta }_2$ and ${\beta }_2$ (${\beta }_0$ is stationary), and RMSE and MAE for ${\beta }_0$, ${\beta }_1$, ${\beta }_2$ and ${\beta }_3$. Model fit diagnostics of AIC, R², RMSE and MAE were also found. Finally, the model residuals ($y-\widehat{y}$) were also extracted to compare the spatial autocorrelation in the residuals arising from the models. To ensure consistency, for AIC, MAE, RMSE and R², were coded and applied to each GAM and MGWR model, rather than using the values generated by the outputs of the mgcv and the GWmodel packages.

3.4. Analysis II: GGP-GAM Tuning with a Larger Dataset

A second set of analyses investigated the impact of increased observation number on tuning via the knots parameter k in order to identify any trade-offs between performance and computation. To do this seven GGP-GAMs were specified each with an increasing number knots. The effects of the increasing k were evaluated through computation time, coefficient accuracy and predictive performance. Again these were implemented in the mgcv package but using the bam function with parallel processing. Each model had the same input data but with the k parameter varying.

3.5. Empirical Example: Brexit Vote

A final analysis to illustrate the approach constructed a GGP-GAM SVC model of the 2016 UK referendum on leaving the EU (the Brexit vote) in England, Wales and Scotland (the Brexit vote) [40]. The spatial distribution of the Leave vote share over 380 Local Authority Districts is shown in Figure 4. This suggests that the overall outcome of a 51.9% majority in favour of ’Leave’ conceals notable regional patterns. A spatially varying coefficient model identifies where the drivers of the outcome vary geographically. For this example, the independent variables from the 2011 UK Census were: Christian the proportion stating their religion as Christian, Degree the proportion with a bachelors degree, No Car the proportion who do not own a car or van, and Younger the proportion of the electorate between the ages of 20 and 44 years old. Each of the dependent variables could influence the Leave vote proportion in a spatially heterogeneous way, in a spatially homogeneous way or not at all. The variables were not rescaled.

4. Results

4.1. Comparing GPP-GAM and MGWR

The comparisons of the estimated coefficients with the true coefficients arising from the 100 GGP-GAM and MGWR models are shown in Figure 5. The boxplots have the outliers removed to emphasise the differences between the median values in the fit measures. For each fit measure (RMSE, MAE, R²), the GPP-GAM provides more accurate estimates of the true coefficients than MGWR, with the accuracy difference increasing with the degree of spatial heterogeneity (i.e. from ${\beta }_1$ to ${\beta }_3$).

The distributions of model fit indicated by AIC show that the GGP-GAMs consistently have lower AIC values compared to the MGWR models. Investigations of the large differences in AIC values revealed them to be driven by differences in the residuals, rather than the number of model parameters. These are summarised in Table 1.

The spatial variation in the residuals was also compared. A Moran’s I statistic was generated for the residuals arising from each GGP-GAM and MGWR model. Their distributions are summarised in Table 2. The MGWR residuals have higher positive values indicating stronger positive spatial autocorrelation (i.e. clustering), while the GGP-GAM residuals have more negative values indicating stronger negative spatial autocorrelation (i.e. dispersion).

The accuracy of the predicted measures of y (i.e. $\widehat{y}$) arising from the two models are compared in Figure 6, using RMSE, MAE and R². The GGP-GAMs generate more accurate predictions, reflecting their better of estimation of the true coefficients.

Model calibration can be investigated. The GGP-GAM spline smoothing parameters are part of a covariance function that penalizes model complexity, while the covariate-specific bandwidths in MGWR describe the scale of each predictor-to-response relationship, thereby reflecting the scale of process spatial heterogeneity. These are summarised in Table 3. The homogeneity of these, as indicated by the inter-quartile ranges, reflect the convergence of the bandwidth optimisation in the MGWR models and the optimisation of the spline smoothing parameters in the GGP-GAM models.

4.2. A Single GGP-GAM in Detail

It is instructive to examine the nature of a GGP-GAM and its coefficients in detail. The 51st model was randomly selected and its properties examined. Diagnostics summaries of the model SVCs are shown in Table 4. These are generated by the gam.check function in the mgcv package, summarise the results of the model optimisation procedures and allow the basis dimensions defined by k to be evaluated. Wood [24] notes that a “low p-value and a k-index of $<1$) may indicate that k is too low, especially if effective degrees of freedom (EDF) are close to k”. Similar results were found for other models that were investigated. It is evident that for these models, a k value of 155 is adequate: the EDFs are lower than k and all the p-values are high.

It is also possible to examine the fixed (parametric) coefficient estimates. These can be considered as the linear model terms, similar to the outputs of a standard OLS regression and are shown in Table 5. The effect of the large k is evident: the Intercept has a significant global effect and the global relationships between the response and the covariates are reduced to zero.

Table 6 summarises the smooth terms, the geographic GP splines (for each SVC). The full set of coefficients are not printed because there many coefficients for each spline, one for each basis function. The edf (effective degrees of freedom) indicates the spline complexity, with higher edf values indicating increasing non-linearity in the predictor-to-response relationship. For example, an edf of 1 suggests a linear relationship, an edf of 2 a quadratic one, etc. In this context, these values apply to each covariate across a 2-dimensional space defined by $(X,Y)$, and the p-values indicate the significance of any spatial variation in the coefficient estimates - i.e. whether they vary significantly over space. In this case the SVCs for ${\beta }_1$, ${\beta }_2$ and ${\beta }_3$ are locally significant but the Intercept (${\beta }_0$) is not.

It also is possible to examine whether and how the relationships between y and the covariates vary spatially - i.e. the SVCs. Table 7 summarises these for the single GGP-GAM model. These reflect the relatively low spatial variation in the Intercept and its globally significant effect, and the relatively high spatial variation in the coefficient estimates for $x_1$, $x_2$ and $x_3$.

A final investigation maps the estimated GGP-GAM and MGWR SVCs. These are shown for the 51st GGP-GAM and MGWR models in Figure 7 alongside the true coefficients. The GGP-GAM and MGWR coefficients are estimated from the same input data, and the shading range was intentionally set to the same as Figure 1. The estimated coefficient surfaces for ${\beta }_0$, ${\beta }_1$, ${\beta }_2$ and ${\beta }_3$ indicate better performance the GAM based approach over MGWR. The grey areas in the MGWR coefficient estimates indicate predicted values outside of range of -2.2 to +2.2.

4.3. GGP-GAM Tuning with a Larger Dataset

A second set of analyses investigated the impact of increased observation number and tuning via the knots parameter for 250,000 observations over a 500 by 500 grid (Figure 2. A single set of simulated predictors $x_1$, $x_2$ and $x_3$ were again calculated from random normal distributions rescaled to have a range of [0, 1] with a random normally distributed error term $ϵ$ rescaled to [0, 0.25]. The response variable y was calculated directly from these and the coefficients in Figure 2. Seven GGP-GAM models were constructed but specified with different vales of k: 100, 250, 500, 750, 1000, 1500 and 2000. The evaluations included computation time, coefficient accuracy and predictive performance.

Model characteristics and model fits are shown in Figure 6. Model deviance (unpenalized), the effective model residual degrees of freedom and the estimated variance parameter all decrease with increased k as might be expected. The model fits also increase with as the degree of tuning increases with k. There is a distinct elbow to many of these trends suggesting suggests that a trade-off is possible between more complex models and increased computation times with higher values of k.

Figure 8. The GGP-GAM model metrics and fits with increasing number of knots k.

Figure 8. The GGP-GAM model metrics and fits with increasing number of knots k.

Table 8 summarises the residuals arising from the different GGP-GAMs. They have similar central tendencies and the inter-quartile ranges indicates that although the models with higher values of k have a lower variation, the differences are not dramatic. Again this reinforces the possibility of a model complexity trade-off with accuracy.

The fixed parametric coefficient estimates of the seven different models were investigated to examine the degree to which the effects observed in Table 5 for the small single GGP-GAM model are also found (i.e. a significant intercept and insignificant covariates whose coefficients were zero). These are shown for models constructed from the larger dataset in Table 9. Whilst the intercept ( ${\beta }_0$ ) has a similar value for each model and is always significant, the trend in the other covariates is more variable, especially at lower values of k. However, they are not significant for the models when $k\ge 1000$.

Table 10 summarises the smooth terms for the seven models, the geographic GP splines. Recall that the effective degrees of freedom (edf) summarises the complexity of the spline smooths and that p-values indicate whether the coefficients are locally significant and whether their respective covariate relationship with the target variable y varies locally. In the model constructed for the single smaller dataset summarised in Table 6, the spatially varying intercept was not significant locally, but all the other covariates were. The same pattern is found with the larger data, but notice how the effective degrees of freedom increases with k.

The GGP-GAM spline smoothing parameters from the different models were examined. Recall that the homogeneity of these indicate the convergence of the models as determined by k which define the spline basis dimensions. These are summarised in Table 11 for the different models. Two trends are evident in the smoothing parameter values across all values of k: the heterogeneity of the Intercept across the different values of k, reflecting its global significance and and local in-significance, in Table 9 and Table 10 respectively, and the homogeneity within each covariate, with values of the same order of magnitude across values of k.

In summary, series of GGP-GAMs were constructed using a larger simulated dataset, with different values of k. As k increases the number of spline bases increase and the models take longer to fit, indicating the trade-off between model complexity and computation time. A distinct elbow was found with increasing k in model deviance, residual degrees of freedom, and estimated variance as well as in the model fit measures (AIC, R², RMSE, MAE) suggesting possible trade-offs at around $k=500$ to $k=750$ (Figure 6). The residuals of this model were found to be similar to those of the models with high values of k (Table 8), where greater tuning might be expected to result in stronger model performance. Except for the Intercept, the global parametric coefficients flatten out or are insignificant at all values of k (Table 9) and all the covariates are locally significant at all values of k (Table 10), supporting trade-offs at values of k between 500 and 750 in this case.

4.4. Empirical Example: Brexit Vote

A final empirical analysis compared GGP-GAM and MGWR SVC models of the 2016 UK referendum on leaving the EU in England, Wales and Scotland (Figure 4). The model summaries are shown in Table 12 and in this case both models perform well in terms of fit metrics (R², AIC, MAE, RMSE) with the MGWR model marginally out-performing the GGP-GAM model, and the GGP-GAM model running more efficiently. However perhaps of greater interest are the minor differences in coefficient estimates as summarised in Table 13 and Table 14 and mapped in Figure 9. The tables include diagnostic metrics of the significance of the GGP-GAM local coefficients and the MGWR bandwidths whose size indicate the degree of localness (their theoretical maximum is 1198 km).

Overall, when the central tendencies and inter-quartile ranges are examined, both sets of coefficient estimates have similar values and ranges. Considering each in turn, a number of observations can be made:

• Intercept (${\mathit{\beta }}_{\mathbf{0}}$): This is not locally significant in the GGP-GAM model but is in its parametric form (not shown). The MGWR model indicates that it has a highly localised (i.e. spatially varying) relationship with a relatively small bandwidth. Both sets of coefficient estimates are positive and have similar values and ranges.
• Christian: Both sets of coefficient estimates indicate a negative association with the Leave vote in Scotland and parts of North Wales and a positive one in England. It is locally significant in the GGP-GAM model and exhibits moderate local variation in the MGWR model, with a bandwidth of 159 km.
• Degree: This is locally significant in the GGP-GAM model and is negatively associated with the Leave vote throughout the study area in both models. It indicates moderate local variation in the MGWR bandwidth (204 km).
• No Car: This is locally significant in the GGP-GAM model. It is mostly negatively associated with the Leave vote share in both models and indicates similar areas of positive association with the Leave vote share in the North. The MGWR bandwidth indicates moderate local variation (172 km).
• Younger: This is not locally significant in the GGP-GAM model In the MGWR, its bandwidth (1196 km) indicates that it is globally (evenly) associated with the Leave vote share.

In summary the two models have similarly high fit and accuracy metrics but they suggest subtly different process spatial heterogeneities and non-stationarities in the relationships between the Brexit Leave vote share as the target variable and the different predictor variables. These are more apparent in the extremes of the study area, potentially reflecting the mechanics of the moving window approach in GWR-based models.

5. Discussion

There is increasing interest in SVC models because of their ability to accommodate, capture and describe process spatial non-stationarity via regression coefficient estimates that are allowed to vary over space. Data are increasingly spatial (i.e. the location of observation is included in some form) and consequently more researchers are working with spatial data. SVC models are attractive because they provide an indication of how and where statistical relationships vary over space and the results can be mapped, providing insights about the spatial heterogeneity of the process being examined.

Some initial work has proposed SVC models using GAMs with Gaussian Process smooths that include observation location - a Geographical GP GAM (GGP-GAM) - and a comparative study with MGWR was undertaken [33]. This used a different, less spatially complex, set of simulated data to those presented here and accepted the default settings for GAMs in the mgcv package [35]. In this study, GGP-GAM and MGWR were applied to a more spatially complex set of simulated data (i.e. with much greater spatial heterogeneity) and the GGP-GAMs were tuned by specifying a sufficiently large number of knots such that the properties of the smooths converged. This was because the number of knots defines the basis dimension in the GP splines and the high value ensured sufficient degrees of freedom. in each of the smooths. MGWR tuning was undertaking through bandwidth optimisation.

For the first part of this study the GGP-GAMs performed better than MGWR across a range of fit measures, although a higher degree of variation was found in the GGP-GAM smoothing parameters than in the MGWR bandwidths (Table 3). The GGP-GAM generated more accurate coefficient estimates than MGWR, a finding that was confirmed when the MGWR analyses were re-run with fixed distance bandwidths rather than the adaptive ones. On investigation, the difference in performance was found to be due to the relatively higher residuals arising from the MGWR models, although objectively both models performed well. This resulted in the MGWR model struggling to estimate the known, true coefficients (Figure 7) and suggests that MGWR models (and potentially other geographically weighted approaches) may struggle with more complex, highly localised spatial heterogeneities, due to their rudimentary kernel / moving window based formulation. Further investigation is required to confirm this assertion.

The GGP-GAM was applied to a larger simulated dataset and a series of models were constructed with varying values of k, defining the number of knots. The results indicate excellent measures of fit at all values of k, and a degree stability in the GGP spline smooths was found for $k\ge 750$, suggesting that acceptable degrees of trade-off are possible between model complexity, computation times and GGP-GAM stability. Identifying the value of k at which model performance stabilises has operational benefits: the value of k determines the upper limit on the effective degrees of freedom linked to spline smoothness and, consequently, the spline bases. Increasing k also prolongs computation time. However, determining the optimal k algorithmically for a given case study is challenging. A common approach is to explore various k values and increasing and decreasing them as necessary. Here it was shown that the gains at values of $k\ge 750$ resulted in marginal performance gains at the cost of increases in computation time (a few minutes vs. a few hours).

The final analysis constructed and compared GGP-GAM and MGWR SVC models of the UK’s Brexit referendum vote to Leave the European Union in 2016 with a small number of socio-economic variables. Both approaches generated similarly accurate models in terms of the ability to model the Leave vote share. The spatially varying coefficient estimates they generated indicated minor differences locally, with perhaps the greatest differences evident in the extremes of the study area. Consider, for example, the coefficients for Christian and No Car in the Shetland Islands in Figure 9: these are positive in the GGP-GAM model and negative in the MGWR model, reflecting the smoothing effect of the moving window approach and potentially the adaptive distance kernel. There were also difference in the predictor variables that were found to have a local effect on the target variable. In sum, the MGWR and GGP-GAM approaches generated similar understandings of the process in this empirical case study but with differences in the geographical extremes. This may be due to the way that GWR based approaches capture spatial non-stationarity: they move a weighted window across observations and estimate coefficients at each location rather than quantifying local variation in ’data’ space [22]. Further work will seek to unpick these inferential differences in the face of each model’s predictive power.

Other of areas of work are suggested by the findings of this study. Future work will explore how the GGP-GAM approach may be further calibrated through out of sample predictions, over randomly selected spatial data subsets. MGWR approaches are unable to do this because of their iterative back-fitting approach to model calibration. As a consequence they cannot be applied to ’new’ data. It will also continue to extend the stgam R package [41]. Initial work on spatially varying coefficient modelling using GGP-GAMs simply coded the analysis in R. But as the analysis progressed, the stgam package was developed to contain contains tools to create varying coefficient models based on GAM GP smooths. The functionality of this toolkit will continue to be expanded. One area of development includes explicit metrics of the scale of spatial dependencies in the coefficient estimates. In MGWR this is provided by the kernel bandwidth and in GGP-GAMs a similar indication of scale could be to determine the variogram range for each spatially varying coefficient. Additionally, the opportunities arising from the core GAM functionality of the proposed approach [7,8,35] will be explored. GGP-GAMs are more flexible than GWR based approaches. They are inherently multiscale through their specification, they are able to handle different responses as well as having options for handling collinearity, outliers, heteroskedastic and autocorrelated error terms (options that are not readily available for MGWR, aside from recent work on MGWR with autocorrelated errors [42] and Poisson responses [6]. Finally, extensions into space-time will be explored to investigate how to model space-time dependencies [43], potentially using model averaging [44] to determine optimal space-time model form. The ability to construct coefficient models in which the coefficients are allows to vary over space and / or time is key in, for example, analysis of resilience to climate change where the aim is to determine the varying effects drivers of changes in resilience and to predict tipping points.

6. Conclusion

This study describes spatially varying coefficient modelling using GAMs with Gaussian Process splines parameterised with observation location, referred to as the Geographical Gaussian Process GAM (GGP-GAM). In a series of analyses, the GGP-GAM was applied to simulated data with known spatial heterogeneities and found to performed better than MGWR, the SVC brand leader. This raised questions about the ability of kernel based approaches like GWR to handle highly localised process. The GGP-GAM was then applied to a larger datasets to investigate model tuning. Processes for determining acceptable trade-offs between model complexity, computational efficiency and stability were identified, with respect to the number of knots used in the GP smooths. Finally, the GGP-GAM and MGWR were applied to en empirical case study of the 2016 UK Brexit referendum vote. Both models performed well in terms of fit metrics with the MGWR model marginally out-performing the GGP-GAM model. The coefficient estimates generated by the two approaches for each covariate were broadly similar in terms of their sign (positive / negative), their mapped spatial distributions, and the degree of local heterogeneity they suggested, However, differences were found in the geographical extremes likely due to the mechanics and smoothing effect of the moving window approach in GWR based models. Further will continue to the develop the GGP-GAM approach, including the R package that has been created out of this work [41], handling responses with non-Gaussian distributions, and the extension into the temporal domain to generate spatially and temporally varying coefficient models.