2.1. Context

GW methods are primarily used to investigate spatial heterogeneity, where the form of the heterogeneity reflects the objective of the under-lying statistic or model. For example, a GW standard deviation (from GW summary statistics) investigates spatial change in data variability; a GW regression investigates spatial change in response and predictor data relationships; a GW variogram (16) investigates spatial change in spatial dependence. In all cases, a moving window weighting technique is used, where local models are calibrated at (sampled or un-sampled) locations (i.e. the window’s center). For an individual model at some calibration location, all neighboring observations are weighted according to the properties of a distance-decay kernel function, and the model is locally fitted to this weighted data. Thus, the geographical weighting solely applies to the data in all GW methods, where each local model is fitted to its own GW data (sub-) set. The size of the window over which this localized model might apply is controlled by the kernel function’s bandwidth. Small bandwidths lead to more rapid spatial variation in the results, while large bandwidths yield results increasingly close to the global model solution. The GW modeling paradigm encompasses many methods; each locally adapted from a global form.

2.2. Building the weighting matrix

Key to GW modeling is the weighting matrix, which sets the spatial dependency in the data. Here, is a n × n diagonal matrix (where n is the sample size) denoting the geographical weighting attached to each observation point, for any model calibration point i at location . Thus, a different weighting matrix is found at each calibration point. This matrix is determined according to the following three key elements: (i) the type of distance metric; (ii) the type of kernel weighting function; and (iii) the kernel weighting function’s bandwidth. For the distance metric, GWmodel permits the Minkowski family of distance metrics, where the power of the Minkowski distance p, needs to be specified. For example, if p = 2 (the default), then the usual Euclidean distance metric results; or if p = 1, then the Manhattan (or Taxi cab) distance metric results. As an example, Lu et al. (26) fit GW regressions using Euclidean and non-Euclidean distance metrics.

Four of the five kernel weighting functions available in GWmodel are defined in Table 1. Each function includes the bandwidth parameter (r or b), which controls the rate of decay. All functions are defined in terms of weighting the sample data, where j is the index of the observation point and d_ij the distance between the points indexed by i and j. For the box-car and bi-square functions, the bandwidth r can be specified beforehand (i.e. a fixed distance) or specified as the distance between the point i and its nearest neighbor, where N is specified beforehand (i.e. an adaptive distance). The bi-square function gives fractional decaying weights according to the proximity of the data to each point i, up until a fixed distance or a distance according to a specified nearest neighbor. The local search strategy for this and the box-car function is simply N neighbors within a fixed radius r or N nearest neighbors for an adaptive approach. Both functions can suffer from discontinuity, although the bi-square function can be defined with a bandwidth that uses all of the data to minimize such problems.

The Gaussian and exponential functions are continuous and use all the data. Their weights decay according to a Gaussian or exponential curve. According to the bandwidth set, data that are a long way from the point i receive virtually zero weight. The key difference between these functions is their behavior at the origin. Usually these continuous functions are defined with a fixed bandwidth b, but can be constructed to behave in an adaptive manner. The bi-square function is useful as it can provide an intermediate weighting between the box-car and the Gaussian functions. To get similar weights from the bi-square and Gaussian functions, the bandwidths rand b can be approximately related by . For all functions, if r or b is set suitably large enough, then all of the data can receive a weight of one and the corresponding global model or statistic is found.

Bandwidths can be: (a) user-specified, when there exists some strong prior belief to do so; (b) optimally (or automatically ) specified using cross-validation and related approaches, provided there exists an objective function (i.e. the method can be used as a predictor); or (c) user-specified, but guided by (b) where an automated approach is not viewed as a panacea for bandwidth selection (27). In GWmodel, automated bandwidths can be found for GW regression, GW regression with a locally compensated ridge term (to address local collinearity problems), generalized GW regression and GW PCA.

2.3. Study data-sets

GWmodel comes with five example data-sets, these are: (1) Georgia, (2) LondonHP, (3) DubVoter, (4) EWHP, and (5) USelect. For this article’s presentation of GW models, we use as case studies, the DubVoter and USelect data-sets.

3.1. GW summary statistics

Although simple to calculate and map, GW summary statistics can act as a vital pre-cursor to an application of a subsequent GW model. For example, GW standard deviations will highlight areas of high variability for a given variable; areas where a subsequent application of a GW PCA or a GW regression may warrant close scrutiny. For attributes x and y at any location i where w_ij accords to a kernel function of Section 2, definitions for a GW mean, GW standard deviation and GW correlation coefficient are, respectively

(1)

(2)

and

(3)

with the GW covariance

(4)

3.2. GW PCA

In a PCA, a set of m correlated variables are transformed in to a new set of m uncorrelated variables called components. Components are linear combinations of the original variables that can allow for a better understanding of sources of variation and trends in the data. Its use as a dimension reduction technique is viable if the first few components account for most of the variation in the original data. In a GW PCA, a series of localized PCAs are computed, where the local component outputs (variances, loadings and scores) are mapped, permitting a local identification of any change in structure of the multivariate data. This local exploration can pinpoint locations where results from a PCA are inappropriate. This in turn, allows for better-informed model decisions for any analysis that may follow, such as a clustering or regression analysis when orthogonal input data are required. GW PCA can assess: (i) how (effective) data dimensionality varies spatially and (ii) how the original variables influence each spatially varying component.

More formally, if an observation location has coordinates (u_j, v_j), then a GW PCA involves regarding a vector of observed variables x_j as having a certain dependence on its location j, where and are the local mean vector and the local variance-covariance matrix, respectively. The local variance-covariance matrix is

(5)

where X is the n × m data matrix; and is a diagonal matrix of geographic weights, generated by a kernel function of Section 2. The kernel’s bandwidth can be user-specified or found optimally via cross-validation (13). To find the local principal components at location , the decomposition of the local variance-covariance matrix provides the local eigenvalues and local eigenvectors (or loading vectors) with

(6)

where is a matrix of local eigenvectors and is a diagonal matrix of local eigenvalues. A matrix of local component scores can be found using

(7)

If we divide each local eigenvalue by , then we find a localized proportion of the total variance (PTV) in the original data accounted for by each component. Thus, at each location j for a GW PCA with m variables, there are m components, m eigenvalues, m sets of component loadings (with each set m × m), and m sets of component scores (with each set n × m). We can obtain eigenvalues and their associated eigenvectors at un-observed locations, but as no data exists for these locations, we cannot obtain component scores.

3.3. Monte Carlo tests for non-stationarity

For GW summary statistics and GW PCA, Monte Carlo tests are possible that test for non-stationarity (1, 13). Tests confirm whether or not the GW summary statistic or aspects of the GW PCA are significantly different to that found by chance or artifacts of random variation in the data. Here, the sample data are successively randomized and the GW model is applied after each randomization. A basis of a significance test is then possible by comparing the true result with results from a large number of randomized distributions. The randomization hypothesis is that any pattern seen in the data occurs by chance and therefore any permutation of the data is equally likely.

As an example for GW correlation, the test proceeds as follows: (i) calculate the true GW correlation at all locations; (ii) randomly choose a permutation of the data where the coordinates are kept in the same pairs, as are the chosen attribute pairs; (iii) calculate a simulated GW correlation at all locations using the randomized data of (ii); (iv) repeat steps (ii) and (iii), say 99 times; (v) at each location i, rank the one true GW correlation with the 99 simulated GW correlations; (vi) at each location i, if the true GW correlation lies in the top or bottom 2.5% tail of the ranked distribution then the true GW correlation can be said to be significantly different (at the 95% level) to such a GW correlation found by chance. The results from this Monte Carlo test are mapped.

For a GW PCA, a similar procedure is followed where the test evaluates whether the local eigenvalues vary significantly across space. Here, the paired coordinates are successively randomized amongst the variable data-set and after each randomization, GW PCA is applied (with an optimally re-estimated bandwidth) and the standard deviation (SD) of a given local eigenvalue is calculated. The true SD of the same local eigenvalue is then included in a ranked distribution of SDs. Its position in this ranked distribution relates to whether there is significant (spatial) variation in the chosen local eigenvalue. The results from this Monte Carlo test are presented via a graph.

3.4. Examples: GW correlations

For a demonstration of an analysis using a GW summary statistic, we calculate GW correlations to investigate the local relationships between: (a) voter turnout (GenEl4004) and LARent and (b) LARent and Unempl. In the former case, the correlations provide a preliminary assessment of relationship non-stationarity between the dependent and an independent variable of a GW regression of Sections 4 and 5. In the latter case, the correlations provide an assessment of local collinearity between two independent variables of such a GW regression (33). In both cases, we specify a bi-square kernel. Furthermore, as the spatial arrangement of the EDs in Greater Dublin is not a tessellation of equally sized zones, it makes sense to specify an adaptive bandwidth, that we user-specify at N = 48. This entails that the bi-square kernel will change in radius, but will always include the closest 48 EDs for each local correlation. Bandwidths for GW correlations cannot be found optimally using cross-validation (see also (23)). We also conduct the corresponding Monte Carlo tests for the two GW correlation specifications. Here, we use the GWmodel functions gwss to find GW summary statistics and montecarlo.gwss to conduct the Monte Carlo tests. The output of our GW correlation analysis is visualized in Figure 1.

The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models

All authors

Binbin Lu, Paul Harris, Martin Charlton & Chris Brunsdon

https://doi.org/10.1080/10095020.2014.917453

Published online:

28 May 2014

Figure 1. GW correlations and associated Monte Carlo tests for: (a) GenEl2004 and LARent; and (b) LARent and Unempl. Global correlations are −0.68 and 0.67, respectively.

Display full size

From Figure 1(a), the relationship between turnout and LARent appears non-stationary, where this relationship is strongest in areas of central and southwest Dublin. Here, turnout tends to be low, while local authority renting tends to be high. The associated Monte Carlo test suggests many instances of unusual relationships, such as those found in the north, that are unusually weak. From Figure 1(b), consistently strong positive correlations between LARent and Unempl are found in three distinct areas of Greater Dublin; areas where local collinearity effects in the GW regression of Sections 4 and 5 are likely to be a cause for concern (see also (25)).

3.5. Examples: PCA to GW PCA

For applications of PCA and GW PCA, we investigate these eight variables: DiffAdd, LARent, SC1, Unempl, LowEduc, Age18_24, Age25_44 and Age45_64. We standardize the data and specify the PCA with the covariance matrix. The same (globally) standardized data is also used in the GW PCA calibration, which is similarly specified with (local) covariance matrices. The effect of this standardization is to make each variable have equal importance in the subsequent analysis (at least for the PCA case).² The PCA results (PTV data and loadings) are found using scale and princomp R functions, and are presented in Table 2.

From the PTV data, the first two components collectively account for 61.6% of the variation in the data. From the loadings, components one and two mainly represent older (Age45_64) and affluent (SC1) residents, respectively. However, these results may not reliably represent local social structure, and an application of GW PCA may be useful. Here, a bandwidth for GW PCA is found using cross-validation, where it is necessary to decide a priori on the number of components, k to retain, provided m ≠ k. Thus, we choose to find an optimal adaptive bandwidth using a bi-square kernel with the bw.gwpca GWmodel function, with k = 3. A bandwidth of N = 131 is returned and will be used to calibrate the GW PCA fit. Observe that we now specify all k = 8 components, but we will focus our investigation on only the first two components. This specification ensures that the PTV data is estimated correctly. The GW PCA fit is conducted using the gwpca GWmodel function.

The GW PCA outputs are visualized and interpreted, focusing on: (1) how data dimensionality varies spatially and (2) how the original variables influence the components. For the former, the spatial distribution of local PTV for the first two components can be mapped. For the latter, we look at the change in size and sign of the eight local loadings together, for a given component, at each of the 322 EDs. In this respect, we map multivariate glyphs that have spokes around a central hub in which the length of the spoke corresponds to the size of the local loading, and its color corresponds to the sign (in this case, blue signifies positive and red signifies negative). The glyphs are scaled relative to the spoke with the largest absolute loading. The variable corresponding to each local loading is always in the same place on the glyph, as follows: DiffAdd is at 0° (north); LARent is 45° (northeast); SC1 is 90° (east); Unempl is 135° (southeast), LowEduc is 180° (south), Age18_24 is 225° (southwest), Age25_44 is 270° (west) and Age45_64 is 315° (northeast).

Figure 2(a) presents the local PTV map, where there is clear spatial variation in the PTV data. Higher percentages are located in the south and in central Dublin, whilst lower percentages are located in the north. The PTV data is also generally higher in the local case, than in the global case (at 61.6%). Figure 2(b) presents a multivariate glyph map for the loadings on the first component, where a spatial preponderance of glyphs of one color or another, or larger spokes on the same variables provide a general indication of the structures being represented at each of the 322 EDs. This map is not intended to be scrutinized in detail, but clearly indicates strong local trends in social structure across Greater Dublin. To provide support for our chosen GW PCA specification, the associated Monte Carlo test evaluates whether the local eigenvalues for the first component vary significantly across space. The results are given in Figure 2(c), where the p-value for the true SD of the eigenvalues is calculated at 0.02. Thus, an application of GW PCA is considered worthy as the null hypothesis of local eigenvalue stationarity is firmly rejected at the 95% level, for the dominant first component.

The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models

All authors

Binbin Lu, Paul Harris, Martin Charlton & Chris Brunsdon

https://doi.org/10.1080/10095020.2014.917453

Published online:

28 May 2014

Figure 2. GW PCA results: (a) PTV data for the first two components; (b) multivariate glyphs of the loadings for all eight components; and (c) Monte Carlo test for first component only.

Display full size

4.1. Basic and mixed GW regression

The basic form of the GW regression model is

(8)

where y_i is the dependent variable at location i; x_ikis the value of the kth independent variable at location i; m is the number of independent variables; β_i0 is the intercept parameter at location i; β_ikis the local regression coefficient for the kth independent variable at location i; and ɛ_iis the random error at location i. A key assumption in this basic (and related forms of) GW regression is that the local coefficients vary at the same scale and rate across space (depending on the particular kernel weighting function that is specified). However, some coefficients (and relationships) may be expected to have different degrees of variation over the study region. In particular, some coefficients (and relationships) are viewed as constant (or stationary) in nature, whilst others are not. For these situations, a mixed GW regression can be specified (1, 2). This semi-parametric model treats some coefficients as global (and stationary), whilst the rest are treated as local (and non-stationary), but with the same rate of spatial variation. This model’s general form can be written as

(9)

where are the k_a global coefficients; are the k_b local coefficients; are the independent variables associated with global coefficients; and are the independent variables associated with local coefficients. In a vector-matrix notation, Equation (9) can be rewritten as

(10)

where y is the vector of dependent variables; X_a is the matrix of globally fixed variables; a is the vector of k_a global coefficients; X_b is the matrix of locally varying variables; and b is the matrix of local coefficients. To calibrate this model in GWmodel, we follow that of Brunsdon et al. (2), where a back-fitting procedure is adopted (34). If we define the hat matrix for the global regression part of the model, as S_a; and that for the GW regression part, as S_b; then Equation (9) can be rewritten as

(11)

where the two components, and can be expressed as

(12)

(13)

and thus the calibration procedure can be briefly described in the following six steps:

An alternative mixed GW regression fitting procedure is given in Fotheringham et al. (1), that uses a method from Speckman (35). This alternative is not as computationally intensive as that presented, and is under consideration for future mixed GW regressions in GWmodel.

4.2. Monte Carlo tests for regression coefficient non-stationarity

For a mixed GW regression, difficulties arise when deciding whether a relationship should be fixed globally or allowed to vary locally. Here, Fotheringham et al. (36) adopt a stepwise procedure, where all possible combinations of global and locally varying relationships are tested, and an optimal mixed model is chosen according to a minimized AIC value. This approach is comprehensive, but computationally expensive, and is utilized in the GW regression 4.0 executable software (37). Alternatively, a Monte Carlo approach can be used to test for significant (spatial) variation in each regression coefficient (or relationship) from the basic GW regression model, where the null hypothesis is that the relationship between dependent and independent variable is constant (1, 3, 4). The procedure is analogous to that presented for the local eigenvalues of a GW PCA in Section 3.3, where for the basic GW regression the true variability in each local regression coefficient is compared to that found from a series of randomized data-sets. If the true variance of the coefficient does not lie in the top 5% tail of the ranked results, then the null hypothesis cannot be rejected at the 95% level; and the corresponding relationship should be globally fixed when specifying the mixed GW regression model. Observe, that if all relationships are viewed as non-stationary, then the basic GW regression should be preferred. Conversely, if all relationships are viewed as stationary, then the standard global regression should be preferred. Advances on the mixed GW regression model, where the relationships can be allowed to vary at different rates across space can be found in the study of Yang et al. (38).

4.3. Example: mixed GW regression model specification

We now demonstrate modeling building for mixed GW regression using the DubVoter data. First we calibrate a basic GW regression. We then conduct the Monte Carlo test on this model’s outputs, to gauge for significant variation (or non-stationarity) in each coefficient, including the intercept term. Finally, we fit a mixed GW regression according to the results of the Monte Carlo test. Our regressions investigate the local/global relationships between the response: GenEl2004 and these eight predictors: DiffAdd, LARent, SC1, Unempl, LowEduc, Age18_24, Age25_44 and Age45_64. For both GW regressions, we specify a bi-square kernel with an adaptive bandwidth.

The optimal bandwidth for the basic GW regression is found at N = 109 in accordance to an automatic AICc approach via the GWmodel function bw.gwr. This bandwidth is then used to calibrate the basic GW regression via the GWmodel function gwr.basic. We then conduct the Monte Carlo test (via the GWmodel function montecarlo.gwr), where the results are presented in Table 3.

The results suggest (say, at the 95% level) that the Intercept term together with the DiffAdd, LARent, LowEduc, Age25_44, and Age45_64 variables, should all be fixed as global in the mixed model. Accordingly, the mixed model is calibrated using the GWmodel function gwr.mixed with the same adaptive bandwidth as that found for the basic model. Observe that the geographically varying coefficients in the mixed model are less variable than the corresponding coefficients from the basic model, although the same bandwidth is used. The estimated coefficients from the functions gwr.basic and gwr.mixed can be summarized in a manner that imitates the layout of the GW regression 3.0 executable software (39), as shown in Tables 4–6.

As an example visualization, Figure 3(a–b) present the coefficient surfaces corresponding to Unempl found from the basic and mixed GW regressions, respectively. The spatial variation in this coefficient is clearer greater when using the basic GW regression. Differences in the coefficient surfaces primarily occur in the northwest and southwest areas of Dublin.

The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models

All authors

Binbin Lu, Paul Harris, Martin Charlton & Chris Brunsdon

https://doi.org/10.1080/10095020.2014.917453

Published online:

28 May 2014

Figure 3. (a) Basic; (b) mixed; and (c) heteroskedastic GW regression coefficients estimates for Unempl.

Display full size

5.1. Multiple hypothesis tests with GW regression

For GW regression, pseudo t-values can be used to test, in a purely informal sense, for evidence of local coefficient estimates that are significantly different from zero (e.g. 22). For each coefficient estimate, at location i, the pseudo t-value can be calculated using:

(14)

where is the standard error of . For details on the standard error calculations, see (1). However, as GW regression yields a separate model at each location i, (where each model is calibrated with the same observations but with different weighting schemes), we need to conduct a large number of simultaneous t-tests. This operation is likely to result in high order multiple inference problems, where the likelihood of increased type I errors needs to be controlled. In this respect, various standard approaches are available that adjust each test’s p-value; these include: (a) Benjamini-Hochberg (40), (b) Benjamini-Yekutieli (41), and (c) Bonferroni (42) approaches. All three approaches can be used with GW regression, but a fourth, the Fotheringham-Byrne (43) approach is specifically designed for this purpose. This Bonferroni-style adjustment can be briefly described as follows. Let the probability of rejecting one or more true null hypothesis (i.e. the family-wise error rate, FWER) be denoted by (with m the number of tests). Then the FWER for testing hypotheses about GW regression coefficients can be controlled at or less, by selecting

(15)

where α is the probability of a type I error in the ith test; p_e is the effective number of parameters in the GW regression; np is the number of parameters in each individual local regression (i.e. the same as that found in the global regression); and n is the sample size.

To demonstrate this topic, the results from the basic GW regression of Section 4 are investigated. Here all four adjustment approaches are provided in the GWmodel function gwr.t.adjust, which simply needs to be specified with the results from the GW regression run. Accordingly, we can map the competing test results, where as an example, Figure 4(a) displays the surface of the original (un-adjusted) p-values for the Unempl coefficient; and Figures 4(b–c) and 5(a–b) present the surfaces for the corresponding adjusted p-values (all significant results are colored red). Observe that the standard (Benjamini-Hochberg, Benjamini-Yekutieli, Bonferroni) approaches adjust the p-values to zero or one (i.e. in-significant and significant), whilst the Fotheringham-Byrne approach provides outputs from zero to one. From the five maps, it can be observed that: (1) the un-adjusted p-values are similar to that found from Benjamini-Hochberg and Benjamini-Yekutieli adjustments; and (2) the Bonferroni and Fotheringham-Byrne adjustments provide similar results.

The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models

All authors

Binbin Lu, Paul Harris, Martin Charlton & Chris Brunsdon

https://doi.org/10.1080/10095020.2014.917453

Published online:

28 May 2014

Figure 4. The p-values associated with Unempl from the basic GW regression of Section 4: (a) un-adjusted; (b) adjusted by Benjamini-Hochberg; and (c) adjusted by Benjamini-Yekutieli.

Display full size

The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models

All authors

Binbin Lu, Paul Harris, Martin Charlton & Chris Brunsdon

https://doi.org/10.1080/10095020.2014.917453

Published online:

28 May 2014

Figure 5. The p-values associated with Unempl from the basic GW regression of Section 4: (a) adjusted by Bonferroni; and (b) adjusted by Fotheringham-Byrne.

Display full size

5.2. Local collinearity diagnostics for a basic GW regression

The problem of collinearity amongst the predictor variables of a regression model has long been acknowledged and can lead to a loss of precision and power in the coefficient estimates (44). This issue is heightened in GW regression since: (A) its effects can be more pronounced with the smaller samples that are used to calibrate each local regression; and (B) if the data is spatially heterogeneous in terms of its correlation structure, some localities may exhibit collinearity when others do not. In both cases, (local) collinearity may cause serious problems in GW regression, when none are found in the corresponding global regression. Diagnostics to investigate local collinearity in a GW regression model, include finding: (i) GW correlations amongst pairs of independent variables; (ii) local variance inflation factors (VIFs) for each independent variable; (iii) local variance decomposition proportions (VDPs); and (iv) local (design matrix) condition numbers (CN). Accordingly, the following rules of thumb can be taken to indicate likely local collinearity problems in the GW regression model: (a) absolute GW correlation values greater than 0.8 for a given independent variable pair; (b) VIFs greater than 10 for a given independent variable; (c) VDPs greater than 0.5; and (d) CN values greater than 30. Observe that all diagnostics are found at the same spatial scale as each local regression of the GW regression model and can thus be mapped. These diagnostics are considered an integral part of an analytical toolkit that should always be employed in any GW regression analysis. Figure 6(a) depicts the levels of complexity associated with such investigations (excluding the use of VDPs). Details on the use and merit of these diagnostics, possible model solutions, and critical discussions on this issue with GW regression, can be found in (6, 25, 33, 45–49). For this study, our objective is to simply demonstrate some of the collinearity diagnostics used. The application of a (possible) model solution with GWmodel, for example, via a locally compensated ridge GW regression, is demonstrated in Gollini et al. (25) using the GWmodel function gwr.lcr.

The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models

All authors

Binbin Lu, Paul Harris, Martin Charlton & Chris Brunsdon

https://doi.org/10.1080/10095020.2014.917453

Published online:

28 May 2014

Figure 6. (a) Levels of complexity for different localized collinearity diagnostics; (b) GW correlations; (c) local VIFs; and (d) local condition numbers.

Display full size

For a given GW regression specification, GW correlations, local VIFs, local VDPs and the local CNs can be found using the GWmodel function gwr.collin.diagno. The same local CNs can also be found using the function gwr.lcr. Example maps presenting GW correlations, local VIFs and local CNs are given in Figure 6(b–d), reflecting diagnostics for same the basic GW regression of Section 4. Scales of each map are chosen to highlight our critical thresholds. Clearly, significant collinearity is present in our study GW regression model, where DiffAdd appears to be a major cause with respect to its relationship to Age25_44 in central areas of Dublin. As the local CNs are large everywhere, the simple removal of one variable from the analysis may go some way in alleviating this problem; before proceeding to a more locally focused analysis with some locally compensated model. This of course brings into question the validity of the results of Kavanagh et al. (28), where basic GW regression was applied to this data.

5.3. Heteroskedastic GW regression

Basic GW regression assumes that the error term is normally distributed with zero mean and constant (stationary) variance over the study region (). An extension of GW regression is possible, which allows a non-stationary error variance . This (possibly more realistic) model was first proposed in Fotheringham et al. (1) and has been further extended to a predictive form in Harris et al. (18). Details of such heteroskedastic GW regressions can be found in (1, 18), where an iterative modeling technique is used, requiring the model to converge to some pre-specified level of tolerance. An alternative (parametric) heteroskedastic GW regression can be found the work of Paez et al. (50–55).

For GWmodel, the function gwr.hetero allows the non-parametric version to be specified. Currently, the model is only given in a very rudimentary form, where the kernel function used to control the coefficient estimates is also used to control the local error variances; which are themselves approximated by the squared residuals (e²). Outputs from gwr.hetero are the regression coefficients only, which can be compared to the corresponding coefficients found using a basic GW regression (gwr.basic). As an example, Figure 3(c) displays the coefficient surface for Unempl from the heteroskedastic model. Clearly, there is little difference from the coefficient surface of the basic model (Figure 3(a)). Thus, modeling with a stationary error variance appears reasonable. The use of a non-stationary error variance can be useful, however, for improved measures of uncertainty (18) and for outlier detection (22).

6.1. GW discriminant analysis

The theoretical context for DA is briefly described. Suppose a population, of which each object belongs to k possible categories; and a training set X, where each row vector x_i indicates an observation belonging to category l ; then for an observation vector x, the discrimination rule is to assign x to the lth category with the maximum probability that x belongs to this category, say

(16)

Here, p_jf_j(x) is proportional to the posterior probability of an observation arising from population j once the value of x is known. Now a linear DA (LDA) assumes that the distribution for each category is multivariate normal with a discrimination rule of

(17)

where ∑ is the covariance matrix, q is the number of independent variables in x, and is the mean for population j. This function can be simplified by taking logs and changing signs

(18)

LDA assumes that ∑ is identical for each category, whilst an alternative, quadratic DA (QDA) assumes that ∑ is different for each category j, where its discrimination rule replaces ∑ in Equation (18) with .

GW DA is a direct local adaption of DA, with the chosen discrimination rule varying across space. Here the stationary mean and covariance estimates of DA are replaced with respective GW mean estimates (Equation (1)) and GW covariance estimates (Equation (4)) in the discrimination rules for both LDA and QDA. Thus, a local LDA or QDA can be found at any location (u,v) using GW DA. Bandwidth selection follows a cross-validation approach, where an optimum bandwidth is identified by minimizing this score:

(19)

where is the proportion of incorrect assignments when the observation i is removed from the sample data.

6.2. Example: GW discriminant analysis

To demonstrate a GW DA, we use the USelect data described in Section 2. Here, we calibrate a GW DA using the function gwda in GWmodel, together with a standard DA (LDA) using the function lda from MASS (51). The GW DA is conducted with an adaptive bandwidth (bi-square kernel) with its optimum found using the GWmodel function bw.gwda. The resultant confusion matrices are presented in Table 7. Here, the DA classification accuracy is 72.5%, whilst GWDA provides a slightly improved classification accuracy of 74.0%. An interesting feature is that the global model predicts only one county in the ‘Borderline’ category. Results of the actual presidential election are mapped at the county level in Figure 7(a), where Bush was a clear winner in most counties, while the election was more competitive in areas like Wisconsin and Maine. The classification results using DA and GW DA are mapped in Figure 7(b–c), where the spatial pattern in the GW DA classifications appears marginally closer to the true results than that found with the DA.

The GWmodel R package: further topics for exploring spatial heterogeneity using geographically weighted models

All authors

Binbin Lu, Paul Harris, Martin Charlton & Chris Brunsdon

https://doi.org/10.1080/10095020.2014.917453

Published online:

28 May 2014

Figure 7. (a) Results of 2004 US presidential election. Classification results using: (b) DA; and (c) GW DA.

Display full size