Geographically Weighted Regression Analysis for Spatial Economics Data: A Bayesian Recourse

Abstract

The geographically weighted regression (GWR) is a well-known statistical approach to explore spatial non-stationarity of the regression relationship in spatial data analysis. In this paper, we discuss a Bayesian recourse of GWR. Bayesian variable selection based on spike-and-slab prior, bandwidth selection based on range prior, and model assessment using a modified deviance information criterion and a modified logarithm of pseudo-marginal likelihood are fully discussed in this paper. Usage of the graph distance in modeling areal data is also introduced. Extensive simulation studies are carried out to examine the empirical performance of the proposed methods with both small and large number of location scenarios, and comparison with the classical frequentist GWR is made. The performance of variable selection and estimation of the proposed methodology under different circumstances are satisfactory. We further apply the proposed methodology in analysis of a province-level macroeconomic data of thirty selected provinces in China. The estimation and variable selection results reveal insights about China’s economy that are convincing and agree with previous studies and facts.

Keywords

MCMC model assessment spatial econometrics variable selection

Introduction

For geographically sparse data with inherent spatial variability, estimating coefficients of a regression model for a particular location based on only observations from this location is not feasible due to the small number of observations. The geographically weighted regression (GWR; Brunsdon, Fotheringham, and Charlton 1996) is an important tool to explore spatial non-stationarity of the regression relationship in spatial data analysis. It has been applied to a variety of fields, including geology, environmental science, epidemiology, and econometrics. Fotheringham, Brunsdon, and Charlton (2002) has summarized the basic theory, statistical inference, and bandwidth selection for GWR, and proposed natural extensions of GWR under the generalized linear model framework. The basic idea of GWR is to make use of information from nearby locations. The weighting idea is a natural strategy to use in light of Tobler’s first law that “near things are more related than distant things” (Tobler 1970). In estimating the parameter for one specific location, subjects in the data are weighed according to their distance from this location, with greater weight for closer subjects. Páez, Uchida, and Miyamoto (2002a, 2002b) proposed estimation and inference for GWR under a maximum-likelihood-based framework. Mei, Wang, and Zhang (2006) proposed a mixed geographically weighted regression model which included both spatially-varying and non-spatially-varying coefficients, and gave a testing procedure of important explanatory variables. N. Wang, Mei, and Yan (2008) proposed a local linear-based GWR for the spatially varying coefficient models which can significantly improve GWR. More recently, da Silva and Fotheringham (2016) discussed the multiple testing issue, and proposed a solution which outperforms other solutions such as the Bonferroni procedure under the GWR framework. The aforementioned works discussed the GWR in the frequentist fashion. From the Bayesian perspective, LeSage (2004) proposed a Bayesian GWR, which gives a prior distribution on the parameter vector depending on historical knowledge. The proposed model, however, used cross-validation for bandwidth selection, which relies on a user-specified grid of bandwidth, and is computationally intensive like other cross-validation based methods.

In this paper, we propose Bayesian techniques for the GWR using the likelihood-based approach in Páez, Uchida, and Miyamoto (2002a, 2002b). In addition to Bayesian estimation and inference, a spike and slab (Ishwaran and Rao 2005) prior is applied for variable selection for Bayesian GWR. Furthermore, bandwidth selection and weighting scheme selection are discussed based on prior selection and Bayesian model selection criteria. An introduction to the implementation of GWR based on nimble (de Valpine et al. 2017), a relatively new and powerful R package for Bayesian inference, is presented as an open source repository on GitHub. Our simulation studies showed the promising empirical performance of the proposed methods in both non-spatially varying and spatially varying cases. In addition, our proposed Bayesian approach reveals interesting features of the province-level macroeconomic data in China.

The rest of this paper is organized as follows. In Section “Geographically Weighted Regression,” the GWR and its weighting schemes are discussed. Section “Bayesian Recourse for GWR” gives a detailed discussion of Bayesian inference, variable selection, bandwidth selection, and model assessment for the GWR. Extensive simulation studies are conducted in Section “Simulation Studies” to investigate the empirical performance of the proposed methods. In Section “Real Data Analysis,” we implement our model using province-level macroeconomic data in China from year 2012 to year 2016. Finally, Section “Discussion” contains a brief summary of this paper.

Geographically Weighted Regression

From Brunsdon, Fotheringham, and Charlton (1998), the GWR model can be written as:

y (s) = β_{1} (s) x_{1} (s) + \dots + β_{p} (s) x_{p} (s) + ϵ (s)

where $y (s)$ is the response variable at location s, $β_{i} (s)$ , $i = 1, \dots, p$ are the coefficients of independent variables at location s, and $ϵ (s)$ is the random effect at location s, assumed to follow $N (0, σ^{2})$ . In addition, we also assume $c o v (ϵ (ℓ), ϵ (m)) = 0$ for any $ℓ \neq m$ . Given a weighting function, the weights of each observation can be calculated with the distance between that observation and s. Estimation of coefficients at location s can be formulated in a way similar to the weighted least squares:

\hat{β} (s) = {(X^{⊤} W (s) X)}^{- 1} X^{⊤} W (s) Y

(2)

where X is the $n \times p$ matrix of covariates, Y is the vector n × 1 of responses, and $W (s) = diag (w_{1} (s), \dots, w_{n} (s))$ is a diagonal matrix of the weights.

Spatial Weighting Functions and Distances

We first introduce several spatial weighting functions that can be used in GWR. Notice that the weighting scheme of ordinary least squares can be defined in the following form:

w_{j k} = 1, \forall j, k

(3)

where j represents the location of the observations, and k represents the location for which parameters are estimated. In a global model where observations from all locations are used to estimate one vector of coefficients, each observation is assigned a weight of unity.

A first step to consider locality is to include observations that are only within a certain distance d from the target location, i.e.

w_{j k} = \{\begin{matrix} 1, d_{j k} \leq d \\ 0, o t h e r w i s e \end{matrix}

(4)

where $d_{j k}$ is the distance between locations j and k. This weighting scheme is one of the simplest to calculate. It is, however, a discontinuous function of distance, which can sometimes lead to undesired jumps in the estimated parameter surface. In order to get a continuous weighting function, the exponential function and the Gaussian function can also be used. The exponential weighting scheme can be written as:

w_{j k} = exp (- d_{j k} / b)

(5)

where b is the bandwidth that can be chosen appropriately to control the decay with respect to distance. The Gaussian weighting scheme can be written as:

w_{j k} = exp (- (d_{j k} / b))

(6)

Both (5) and (6) are decreasing functions of d_jk, which, intuitively, indicates that an observation very far away from the location of interest contributes little in the estimation of parameters at this location. In order to provide a continuous, near-Gaussian weighting function up to distance b from the estimation point, and then zero weights for any data point beyond b, Brunsdon, Fotheringham, and Charlton (1996, 1998); Fotheringham, Charlton, and Brunsdon (1998) proposed the bi-square function:

w_{j k} = \{\begin{matrix} {(1 - {(\frac{d_{j k}}{b})}^{2})}^{2} i f d_{j k} \leq b \\ 0 o t h e r w i s e \end{matrix}

(7)

In this function, b serves two purposes simultaneously. It controls the speed of weight decay as the distance $d_{j k}$ increases, and also, as a threshold, as points beyond b are not be considered in estimation for the location of interest. By tuning the threshold b, one can control the number of neighbors that are used to estimate the parameters for the location of interest. The weighting schemes mentioned above are the most popular schemes used in the GWR.

We then briefly discuss different choices of distance function. The Euclidean distance, defined as:

d_{j k} = \sqrt{{({latitude}_{j} - {latitude}_{k})}^{2} + {({longitude}_{j} - {longitude}_{k})}^{2}},

is one of the most popular choices when the precise (latitude, longitude) location of each observation is available. However, in some public health and epidemiology studies, or some socioeconomics studies, data are collected and summarized on a higher level than single observations, such as wards (Brunsdon, Fotheringham, and Charlton 1996) or counties (Xue, Schifano, and Hu 2019), which produces areal data instead of point-reference data. All observations within the same area are assigned the same (latitude, longitude). For example, Hu and Huffer (2020) attributed each county’s observations to its centroid. Note that the Euclidean distance is easily affected by the areas of the administrative divisions, which additionally complicates the process of parameter tuning as there is no golden benchmark measure of distance.

An alternative distance measure when we have areal data is the graph distance (Müller et al. 1987; Bhattacharyya and Bickel 2014). The administrative devisions are regarded as vertices of a graph, denoted as $v_{1}, \dots, v_{n}$ . The graph also includes a set of edges, $E (G) = \{e_{1}, \dots, e_{m}\}$ , where each edge connects a pair of vertices. The graph distance is defined as:

d_{v_{i} v_{j}} = \{\begin{matrix} |V (e)| i f e i s t h e s h o r t e s t p a t h c o n n e c t i n g v_{i} a n d v_{j} \\ \infty i f v_{i} a n d v_{j} a r e n o t c o n n e c t e d \end{matrix}

(8)

where $|V (e)|$ denotes the number of edges in e.

While it remains a subjective problem in choosing appropriate bandwidths and thresholds for the geographical distance-based methods, i.e. one has to decide “how close is close enough,” a natural definition of closeness would derive from the graph distance. Counties sharing a common boundary, i.e. having graph distance 1, are close, while having graph distance greater than 1 indicates “not close,” and observations in these far neighboring counties need to be weighed down. A graph distance-based weighting function would be

w_{j k} = \{\begin{matrix} 1 i f d_{j k}^{G} \leq b \\ f (d_{j k}^{G} b) o t h e r w i s e \end{matrix},

(9)

where d denotes the graph distance, and f is a certain weighting function with bandwidth b. In this work, we choose $f ()$ to be the negative exponential function, i.e.

w_{j k} = \{\begin{matrix} 1 i f d_{j k}^{G} \leq 1 \\ e x p (- d_{j k}^{G} / b) o t h e r w i s e \end{matrix}

(10)

Note that while the bi-square weighting scheme (7) also allows for bandwidth tuning, it excludes observations beyond the value of the bandwidth, or, threshold. For areal data in economics such as the Georgia housing data to be presented later, this could mean the exclusion of a fair amount of observations. The parameter estimates produced based on a small number of observations would nevertheless be highly unstable, which should be avoided. The weighting scheme in (10), by restricting the immediate neighbors to be assigned the same weight and all others positive weight, compared with (7), adds more ensurance of parameter estimation stability.

Bayesian Recourse for GWR

In this section, we propose the posterior estimation, variable selection, and bandwidth selection for the Bayesian GWR model. The proposed methods are implemented with the powerful R package nimble. The code and documentation can be found at GitHub.

Bayesian Estimation for GWR

According to Boscardin and Gelman (1993) and Páez, Uchida, and Miyamoto (2002a, 2002b), we can get the estimation of GWR using Bayesian computation. The likelihood function of this model can be written as:

Y | β (s), X, W (s), σ^{2} (s) \sim MVN (X β (s), σ^{2} (s) W^{- 1} (s))

(11)

where MVN indicates the multivariate normal distribution. In order to have a conjugate posterior distribution, we can set the priors of $β (s)$ and $σ^{2} (s)$ as:

β (s) | Σ_{β} \sim N (0, Σ_{β})

(12)

where $Σ_{β}$ is a diagonal matrix, and

σ^{2} (s) \sim IGamma (α_{1}, α_{2}), j = 1, \dots, p,

(13)

where $α_{1}$ , $α_{2}$ are the hyper-parameters for distributions of $σ^{2} (s)$ . One set of non-informative choices of hyper-parameters is to set $Σ_{β} = 100 I_{p}$ and $α_{1} = α_{2} = 0.01$ (Gelman et al. 2013). The posterior distribution can be written as:

p (β (s), σ^{2} (s) Y, X, W (s)) \propto p (Y β (s), X, W, σ^{2} (s)) \times p (β (s) Σ_{β}) \times p (σ^{2} (s)) .

(14)

According to (14), we can use Markov chain Monte Carlo (MCMC, Gelman et al. 2013) to estimate $β (s)$ and $σ^{2} (s)$ .

Bayesian Variable Selection

We first consider the regression problem for one location. Following the procedure of George and McCulloch (1993), the spike and slab prior for $β_{j} (s)$ can be formulated as:

β_{j} (s) | γ_{j} \sim (1 - γ_{j}) N (0, τ_{j}^{2}) + γ_{j} N (0, c_{j}^{2} τ_{j}^{2}),

(15)

where

P (γ_{j} = 0) = 1 - P (γ_{j} = 1) = p_{j}

(16)

When $γ_{j} = 0$ , $β_{j} (s) \sim N (0, τ_{j}^{2})$ , and when $γ_{j} = 1$ , $β_{j} (s) \sim N (0, c_{j}^{2} τ_{j}^{2})$ . Our interpretation of this prior is: we set $τ_{j}$ small enough so that if $τ_{j} = 0$ , $β_{j} (s)$ would be so small that we can “safely” estimate it as 0; inversely, we set c_j large so that if $γ_{j} = 1$ , we include the $β_{j} (s)$ into our final model. For the prior of $γ_{j}$ , we set $γ_{j} \sim Bernoulli (0.5)$ , which is a non-informative choice.

Bayesian Bandwidth Selection

In GWR, it is important to choose a proper bandwidth for the weighting functions. In the Bayesian approach, a prior can be given to the bandwidth b so that the optimal bandwidth can be simultaneously obtained together with the estimation of other parameters. The prior also depends on which measure of distance is used. A more detailed discussion of distance measures is given in Section “Spatial Weighting Function and Distances.” Using similar ideas as in Boehm Vock et al. (2015), a prior for bandwidth can be set as:

b \sim Uniform (0, D),

where D is the upper limit for the support of the distribution of b. Without any prior knowledge, D can be chosen large enough so that we start from a noninformative prior for the bandwidth, i.e. we start from an approximate global model where observations are always weighed equally. There are also some other choices of prior distributions for the bandwidth, such as the gamma distribution or discrete uniform distribution. If prior information is available about the bandwidth, parameters for the prior distributions can be set to incorporate such information. Our proposed model can be summarized as follows:

Y | β (s), X, W (s), σ^{2} (s) \sim MVN (X β (s), σ^{2} (s) W^{- 1} (s))

(17)

β_{j} (s) | γ_{j}, τ_{j} \sim (1 - γ_{j}) N (0, τ_{j}^{2}) + γ_{j} N (0, c_{j}^{2} τ_{j}^{2})

(18)

τ_{j}^{2} \sim IGamma (α_{1}, α_{2})

(19)

γ_{j} \sim Bernoulli (0.5)

(20)

w_{i} (s) = f (d_{i} | b)

(21)

b \sim Uniform (0, D)

(22)

where “IGamma” denotes the inverse-Gamma distribution, and f is the weighting function introduced in (9), $i = 1, \dots, n$ and $j = 1, \dots, p$ . As we incorporate the prior of b into our model, conjugate posterior distribution for b cannot be obtained. Therefore, we use the Metropolis–Hastings Algorithm (MH; Gelman et al. 2013) to estimate the parameters.

Bayesian Model Assessment

In Section “Spatial Weighting Function and Distances,” we introduced several spatial weighting functions that can be used in the GWR. In order to select the weighting scheme that fits the data best, we apply the most commonly used tools, the Deviance Information Criterion (DIC; Spiegelhalter et al. 2002) and the Logarithm of the Pseudo-Marginal Likelihood (LPML; Ibrahim, Chen, and Sinha 2013), for model selection. The DIC is defined as:

DIC = Dev (\bar{θ}) + 2 p_{D},

(23)

where 0 and 0 represent the parameter of interest and the corresponding posterior mean. The term $Dev (\cdot)$ denotes the deviance function, while p_D is the effective number of parameters in the model, given by $p_{D} = \bar{Dev} (θ) - Dev (\bar{θ}) .$ For the GWR model in our paper, the following deviance function can be specified (Ma, Chen, and Hu 2020):

Dev (β (s), W (s), σ^{2} (s)) = - 2 log f (Y | β (s), X, W (s), σ^{2} (s))

= n log (2 π) + log (σ^{2} (s)) - log (|W (s)|) + {(Y - X β (s))}^{⊤} σ^{- 2} (s) W (s) (Y - X β (s)),

where n is the total number of the observations. Therefore, the DIC for the GWR model can be given as:

DIC = Dev (\bar{β} (s), \bar{W} (s), {\bar{σ}}^{2} (s)) + 2 p_{D} = 2 \bar{Dev} (β (s), W (s), σ^{2} (s)) - Dev (\bar{β} (s), \bar{W} (s), {\bar{σ}}^{2} (s))

where $\bar{β} (s)$ , $\bar{W} (s)$ and ${\bar{σ}}^{2} (s)$ are posterior estimates obtained from MCMC results. A smaller value of DIC indicates a better model. It can be regarded as the Bayesian equivalent of AIC, where the term p_D is the penalty term for model complexity, similar to the p, i.e. dimension of parameter space, in AIC. Similar to AIC, DIC also takes both fitness and model complexity into account simultaneously.

The LPML is constructed based on the Conditional Predictive Ordinate (CPO) values, which are estimates of the probability for observing Y_i given that all other responses have been observed. Let $D_{(- i)} = \{Y_{j} : j = 1, \dots, i - 1, i + 1, \dots, N\}$ denote the observed data with the ith subject response deleted. The CPO for the ith subject is defined as:

{CPO}_{i} = \int f (Y β (s), X, W (s), σ^{2} (s)) π (W (s), β (s), σ^{2} (s) D_{(- i)}) d (W (s), β (s), σ^{2} (s)),

(24)

where $π (W (s), β (s), σ^{2} (s) D_{(- i)}) = \frac{\prod_{j \neq i} f (Y_{j} | β (s), X, W (s), σ^{2} (s)) π (W (s), β (s), σ^{2} (s) | D_{(- i)})}{c (D_{(- i)})},$ and $c (D_{(- i)})$ denotes the normalizing constant. Within the Bayesian framework, a Monte Carlo estimate of the CPO can be obtained as:

{\hat{C P O_{i}}}^{- 1} = \frac{1}{T} \sum_{t = 1}^{T} \frac{1}{f (Y_{i} | X_{i}, β_{t} (s), W_{t} (s), σ_{t}^{2} (s))},

where T is the total number of Monte Carlo iterations. Then an estimate of the LPML is given by:

\hat{L P M L} = \sum_{i = 1}^{N} log \hat{C P O_{i}}

(25)

A model with a larger LPML value indicates that it is preferred.

Simulation Studies

In this section, we present the performance of the proposed estimation and variable selection techniques under scenarios where the covariate effects do not vary spatially, and where the covariate effects do vary spatially. We use the spatial structure of thirty selected provinces in China in our simulations. A map of these provinces with their names is presented in Figure 1(a). Specifically, Hainan province is an island, and is therefore not connected with any others, which makes its graph distance with any other province infinity. It is, however, very close to Guangdong province, and they bear a lot of resemblance in both culture and economic development. Therefore, we modified the adjacency matrix so that Hainan and Guangdong are adjacent. The graph distance matrix is calculated based on the modified adjacency matrix. A visualization of the graph distance matrix is presented in Figure 1(b).

Figure 1.

(a) A map of the thirty selected provinces of China, with their names annotated. (b) Visualization of graph distance matrix for thirty selected provinces in China. Darker color indicates larger graph distance.

Denote the average parameter estimates as ${\bar{\hat{β}}}_{ℓ, m}$ , calculated as

{\bar{\hat{β}}}_{ℓ, m} = \frac{1}{100} \sum_{r = 1}^{100} {\hat{β}}_{ℓ, m, r},

(26)

where ${\hat{β}}_{ℓ, m, r}$ denotes the parameter estimate for the mth coefficient of province $ℓ$ in the rth replicate. The parameter estimates are evaluated based on their bias, standard deviation, mean squared error, and coverage rate of the 95% highest posterior density (HPD) intervals in the following four ways:

mean absolute bias (MAB) = \frac{1}{30} \sum_{ℓ = 1}^{30} \frac{1}{100} \sum_{ℓ = 1}^{100} |{\hat{β}}_{ℓ, m, r} - β_{ℓ, m}|,

(27)

mean standard deviation (MSD) = \frac{1}{30} \sum_{ℓ = 1}^{30} \sqrt{\frac{1}{99} \sum_{ℓ = 1}^{00} {({\hat{β}}_{ℓ, m, r} - β_{ℓ, m})}^{2}},

(28)

mean of mean squared error (MMSE) = \frac{1}{30} \sum_{ℓ = 1}^{30} \frac{1}{100} \sum_{ℓ = 1}^{100} {({\hat{β}}_{ℓ, m, r} - β_{ℓ, m})}^{2},

(29)

mean coverage rate (MCR) = \frac{1}{30} \sum_{ℓ = 1}^{30} \frac{1}{100} \sum_{ℓ = 1}^{100} I (β_{ℓ, m} \in 95 % HPD interval),

(30)

where $β_{ℓ, m}$ is the true underlying parameter, and $I (\cdot)$ denotes the indicator function. These measures are first calculated for each individual province over replicates, and then averaged over provinces.

The variable selection approach is evaluated using both the accuracy rate for a single variable ACC and for the entire model (Model ACC), defined as:

{ACC}_{m} = \{\begin{matrix} \frac{1}{100} \sum_{r = 1}^{100} I (covariate m is selected in the final model) β_{ℓ, m} \neq 0, l = 1, \dots, 30 \\ \frac{1}{100} \sum_{r = 1}^{100} I (covariate m is not selected in the final model) β_{ℓ, m} = 0, l = 1, \dots, 30 \end{matrix},

and

M o d e l A C C = \frac{1}{100} \sum_{r = 1}^{100} I (the exact true underlying model is selected)

To compare with frequentist approach, the same datasets are also fitted using the classical frequentist GWR approach, where the bandwidth selection is made based on minimizing the summation of SSE across thirty provinces over a grid of candidate bandwidths. The details of frequentist bandwidth selection, as well as the parameter estimates obtained using the optimal bandwidth, are presented in Section 1 of the supplemental material. To demonstrate that the graph distance produces credible parameter estimates, and that at the same time it circumvents the additional effort of threshold selection in weighting kernels such as (4) and (10), simulation study is done for the same designs to be presented, with the great circle distance used. The results are reported in Section 2 of the supplemental material. In addition, considering that a total of 150 observations with thirty locations still make a small sample, an additional simulation study with more than 300 locations using the spatial structure of census tracts in Hartford, Litchfield, and Middletown counties in Connecticut has been conducted, and included in Section 3 of the supplemental material.

For both the following simulation studies and the supplemental simulation studies, the effective number of parameters for the frequentist GWR is also calculated as in Brunsdon, Fotheringham, and Charlton (2000). Note that the frequentist GWR is based on one bandwidth only, and only the full model that includes all covariates is fitted, therefore the frequentist p_D should be used as a reference, instead of a criterion for direct comparison.

Simulation Without Spatially Varying Coefficients

Under the scenario where there are no spatially varying coefficients, we generate data using the same set of parameters for all provinces. The independent continuous covariates are generated i.i.d. from the standard normal distribution $N (0, 1)$ , denoted as $X_{1}, X_{2}, \dots, X_{5}$ , and we use the matrix X to denote the covariate matrix with 5 columns, with the ith column being X_i. The response vector Y is generated as $X β + ϵ$ , where $ϵ \sim MVN (0, I)$ . Different choices of $β$ have been used corresponding to different underlying true models. The parameter D for bandwidth is set to be 100. Given that the maximum graph distance in the spatial structure of the 30 selected provinces is 6, a bandwidth of 100 induces a weighting schemes that, even if the distance between one certain province and another province whose parameter estimates we want to obtain, this province gets assigned a relative weight of $exp (- 6 / 100)$ , which is approximately 0.941 and thus approximates a global model where every observation is equally weighed. This ensures that the prior for bandwidth b is sufficiently noninformative. For each province, five observations are generated, resulting in 150 observations per replicate. A total of hundred replicates are performed. For each replicate, a chain of length 10,000 is run without thinning, where the first 2,000 samples are discarded as burn-in. Three parameter settings similar to in Shao (1997) were used, with $β^{⊤} = (2, 0, 0, 4, 8), (2, 2, 0, 4, 8)$ and $(2, 2, 3, 4, 8)$ , respectively. The mean of bandwidths selected in the 100 replicates was also calculated. The results are reported in Table 1.

Table 1.

Average Parameter Estimates, Performance of Parameter Estimates, and Variable Selection Results When There is No Spatial Variation in the Underlying True Parameters. “True” Denotes Whether a Covariate is in the True Model or Not, and ACC Stands for Variable Selection Accuracy Rate.

		$\bar{\hat{β}}$	MAB	MSD	MMSE	MCR	TRUE	ACC(%)	Model ACC(%)	b
Setting 1	$β_{1}$	1.998	0.055	0.074	0.005	0.942	1	100	100	70.112
	$β_{2}$	0.007	0.037	0.046	0.002	0.99	0	100
	$β_{3}$	0.004	0.04	0.049	0.002	0.999	0	100
	$β_{4}$	4.019	0.062	0.077	0.006	0.974	1	100
	$β_{5}$	7.99	0.066	0.085	0.007	0.934	1	100
Setting 2	$β_{1}$	1.998	0.055	0.073	0.005	0.942	1	100	100	70.213
	$β_{2}$	2.012	0.064	0.08	0.006	0.975	1	100
	$β_{3}$	0.003	0.04	0.049	0.002	1	0	100
	$β_{4}$	4.019	0.062	0.077	0.006	0.975	1	100
	$β_{5}$	7.991	0.066	0.085	0.007	0.934	1	100
Setting 3	$β_{1}$	1.997	0.055	0.073	0.005	0.939	1	100	100	70.058
	$β_{2}$	2.012	0.065	0.08	0.007	0.973	1	100
	$β_{3}$	3.005	0.068	0.084	0.007	0.954	1	100
	$β_{4}$	4.018	0.062	0.076	0.006	0.977	1	100
	$β_{5}$	7.99	0.067	0.085	0.007	0.936	1	100

It is rather clear that when there is no spatial variation, the bandwidth is selected to be large, which induces a weighting scheme that assigns close to uniform weight to both nearby provinces and distant provinces. Under all three settings, the variable selection accuracies are all 100% for all five covariates, and the three model selection accuracies are 100% as well. The average effective number of parameters under the frequentist GWR are 11.08, 10.27 and 12.65 under the three settings, while under the Bayesian GWR, the values are 178.00, 177.89 and 177.57, respectively.

Simulation with Spatially Varying Coefficients

For estimation and variable selection in the presence of spatially varying coefficients, we use similar simulation schemes as in Xue, Schifano, and Hu (2019). The graph distance matrix visualized in Figure 1 is transformed using multidimensional scaling (MDS; Cox and Cox 2000) and mapped into a Cartesian space. Denoting the transformed coordinates corresponding to province as $(x_{ℓ}^{c}, y_{ℓ}^{c})$ , the $β$ vector for province is set to

β_{p, ℓ} = \{\begin{matrix} 0 x_{p} n o t i n t h e t r u e m o d e l \\ β_{p} + 0.2 (x_{ℓ}^{c} + y_{ℓ}^{c}) x_{p} i n t h e t r u e m o d e l \end{matrix} .

(31)

The variation pattern is visualized in Figure 2. The estimation and variable selection results for the setting in (31) are presented in Table 2. It can be seen that in all three settings, the MAB, MSD and MMSE are all larger than when there are no spatially varying covariate effects. There have been decrease in the MCR for parameters corresponding to covariates that are in the true model. The variable selection procedure, however, remains quite robust, and in all replicates of simulation are the correct models selected. The average bandwidths selected for the three settings are around sixty seven. This indicates that in the presence of spatial variation, for some replicates of simulation, the bandwidths tend to be large as the gain in stabilizing the parameter estimates for each location dominates the incurred bias. This has also been observed for the classical frequentist GWR, as presented in Supplemental Table 2. The frequentist GWR has an average effective number of parameters of value 19.32, 20.49 and 20.19 under the three settings, and the Bayesian GWR has 179.38, 179.58 and 179.32 instead.

Figure 2.

Visualization of parameter surfaces when the parameters vary according to (31).

Table 2.

Performance of Parameter Estimates, and Variable Selection Results When a Simple Linear Variation Pattern is Present. “True” Denotes Whether a Covariate is in the True Model or Not, and ACC Stands for Variable Selection Accuracy Rate.

		MAB	MSD	MMSE	MCR	TRUE	ACC(%)	Model ACC(%)	b
Setting 1	$β_{1}$	0.107	0.08	0.017	0.778	1	100	100	68.447
	$β_{2}$	0.041	0.05	0.002	0.998	0	100
	$β_{3}$	0.04	0.051	0.003	0.98	0	100
	$β_{4}$	0.111	0.088	0.019	0.761	1	100
	$β_{5}$	0.111	0.087	0.019	0.761	1	100
Setting 2	$β_{1}$	0.107	0.08	0.017	0.786	1	100	100	67.677
	$β_{2}$	0.111	0.088	0.019	0.766	1	100
	$β_{3}$	0.04	0.05	0.003	0.98	0	100
	$β_{4}$	0.111	0.089	0.019	0.766	1	100
	$β_{5}$	0.111	0.086	0.018	0.764	1	100
Setting 3	$β_{1}$	0.107	0.081	0.018	0.786	1	100	100	66.92
	$β_{2}$	0.112	0.089	0.019	0.772	1	100
	$β_{3}$	0.11	0.084	0.018	0.777	1	100
	$β_{4}$	0.112	0.09	0.019	0.766	1	100
	$β_{5}$	0.111	0.087	0.019	0.767	1	100

To study the estimation and variable selection performance under a scenario where regional variation exists, we choose to use the four major economic regions of China proposed during the Eleventh Five-year plan: the west (0), northeast (1), central (2), and east (3) regions. A visualization of these four regions is given in Figure 3. Provinces within each economic region are assigned the same parameter value. The four $β^{⊤}$ ’s under the three simulation settings are given in Table 3.

Figure 3.

Map of selected provinces in China colored by their economic regions proposed during the Eleventh Five-year Plan.

Table 3.

The True Underlying Parameters Used for the Four Regions in Each Setting.

	Region 0	Region 1	Region 2	Region 3
Setting 1	(1.8, 0, 0, 4.2, 7)	(1.5, 0, 0, 3.8, 9)	(2.2, 0, 0, 4, 8.5)	(2, 0, 0, 4, 8)
Setting 2	(1.8, 1.8, 0, 4.2, 7)	(1.5, 1.5, 0, 3.8, 9)	(2.2, 2.2, 0, 4, 8.5)	(2, 2, 0, 4, 8)
Setting 3	(1.8, 1.8, 2.9, 4.2, 7)	(1.5, 1.5, 3.4, 3.8, 9)	(2.2, 2.2, 3.1, 4, 8.5)	(2, 2, 3, 4, 8)

The estimation and variable selection results are presented in Table 4. Again, compared to results in Supplemental Table 3, both approaches yield similar MAB, MSD, MMSE and MCR/MCP. Similar to previously observed, the variable selection in the Bayesian approach effectively reduces the MAB, MSD and MMSE parameter estimates for variables that are not in the true underlying model.

Table 4.

Performance of Parameter Estimates and Variable Selection Results When Regional Variation is Present. “True” Denotes Whether a Covariate is in the True Model or Not, and ACC Stands for Variable Selection Accuracy Rate.

		MAB	MSD	MMSE	MCR	TRUE	ACC(%)	Model ACC(%)	b
Setting 1	$β_{1}$	0.205	0.107	0.061	0.557	1	100	100	69.469
	$β_{2}$	0.041	0.049	0.002	1	0	100
	$β_{3}$	0.043	0.053	0.003	0.999	0	100
	$β_{4}$	0.136	0.111	0.028	0.782	1	100
	$β_{5}$	0.662	0.125	0.554	0.167	1	100
Setting 2	$β_{1}$	0.205	0.108	0.061	0.555	1	100	100	70.156
	$β_{2}$	0.204	0.105	0.061	0.578	1	100
	$β_{3}$	0.044	0.054	0.003	0.999	0	100
	$β_{4}$	0.136	0.114	0.029	0.791	1	100
	$β_{5}$	0.662	0.126	0.554	0.168	1	100
Setting 3	$β_{1}$	0.205	0.111	0.062	0.559	1	100	100	70.491
	$β_{2}$	0.205	0.106	0.061	0.589	1	100
	$β_{3}$	0.146	0.113	0.034	0.743	1	100
	$β_{4}$	0.138	0.116	0.029	0.79	1	100
	$β_{5}$	0.662	0.126	0.554	0.171	1	100

The bandwidths selected average to around 70. This could be due to the fact that a province now has a few neighbors with exactly the same true underlying coefficients, and therefore the weighting function tries to weigh observations in the neighboring provinces equally as the local one. The frequentist GWR performed similarly, and results are included Supplemental Table 3. The average effective number of parameters under the frequentist GWR are 17.27, 16.73 and 15.28 under the three settings, while under the Bayesian GWR, the values are 178.35, 178.25 and 178.38, respectively.

Real Data Analysis

The proposed Bayesian GWR model is used to analyze province-level macroeconomic data in 30 selected provinces of China from year 2012 to year 2016, i.e. we have 150 observations in total. The Gross Domestic Product (GDP, in billions of CNY) is used as the spatial response variable (Y). Five covariates, including the resident population in millions (X₁), the urban population in millions (X₂), the fixed asset investment in the whole society in billions of CNY (X₃), total export value in billions of USD (X₄), and total import value in billions of USD (X₅), are incorporated in the model. The five-year means of the variables for each province are shown in Figure 4. Following the common practice in economics to account for long-tails (see, e.g. Wooldridge 2015), we take the logarithm of GDP before model fitting. All five covariates are continuous, and are therefore standardized before model fitting.

Figure 4.

Five-year means for variables in selected provinces of China.

The proposed Bayesian GWR model (17)–(22) is fitted on this dataset. Priors $σ^{2} (s) \sim IGamma (1, 1)$ and $β_{0} (s) \sim N (0, 1)$ are given, and we set $τ_{j}^{2} = 0.001$ , $c_{j}^{2} = 10, 000$ the common practice in spike-and-slab model selection, and D = 100 so that we start from an approximately uniform weight over all provinces. The same graph distance matrix as in Section was used. The length of chains was selected to be 5,000, with the first 2000 as burn-in. The unity, exponential, and Gaussian weighting schemes were considered, and the DIC and LPML were used to select the best among the three for this particular dataset. The DIC and LPML values as well as the effective number of parameters (p_D; Spiegelhalter et al. 2002) for these three weighting schemes are shown in Table 5. It can be seen that the Gaussian weighting scheme yields the smallest DIC value and the largest LPML value, indicating that the model with a Gaussian weighting scheme is selected as the best model among the candidate models.

Table 5.

DIC and LPML Values for Different Weighting Schemes.

	Unity Scheme	Exponential Scheme	Gaussian Scheme
DIC	12,584.45	12,525.65	12,511.21
LPML	−6,875.99	−6,876.55	−6,871.79
p_D	214.47	180.68	190

Under the GWR model with a Gaussian weighting scheme, the posterior modes of the indicators $γ_{j}$ are $(1, 1, 0, 1, 0)$ , respectively, which indicates that the covariates X₁, X₂ and X₄ are selected, while covariates X₃ and X₅ can be excluded from the regression model. Specifically, in our model, the number of resident population, the number of urban population, and total export value can help explain the change of GDP in each province. The posterior estimation results of the parameters under the Bayesian GWR model with a Gaussian weighting scheme are presented in Figure 5. The geographical variation in the parameters is rather obvious. We can see that the number of resident population and the total export value have significant negative impact on the increase of GDP, while the number of urban population has significant positive impact. For the Gaussian weight function, the posterior estimate of the bandwidth b is 9.40, which indicates that the most distant provinces are assigned a relative weight of 0.665 in the estimation for one particular province. The impact of resident population on GDP is larger in north China than in southeast China. Comparing this to the population density map in China made by the Center for Geographic Analysis at Harvard University (worldmap.havard.edu/maps/11756), it can be seen that the influence is bigger in less populous provinces, which is in accordance with our intuition. The increasing effect of urban population on GDP is smaller in southeast China than in northwest China. Considering the relatively higher urbanization in south and east China (L. Wang et al. 2012), this can be explained by the “decreasing marginal effect” in economics. Export appears to be more important to provinces in west China than in eastern areas, which can also be explained by the fact that east China is relatively more developed, and have a more diversified source of GDP.

Figure 5.

Plot of parameter estimates for the three covariates in the final model in each selected province.

For comparison, classical frequentist GWR in (2) is also fitted on the dataset. We report the parameter estimates, together with plots, in the supplemental material. Particularly, the two covariates dropped by our Bayesian variable selection approach have the smallest absolute parameter values among all five covariates, which indicates that our proposed approach is indeed capable of picking out the most influential factors. The parameter estimates are also plotted on maps as in Figure 5. There are slight differences in the values of parameter estimates, which is partly due to the fact that we only have 150 observations (thirty provinces, five years). It is, however, worth noticing that the trend of variation is consistent between the frequentist and Bayesian approaches.

Discussion

We developed a likelihood-based Bayesian approach to estimate regression coefficients in conjunction with spike and slab variable selection for geographically sparse data. The selection of bandwidth is discussed for a wide choice of weighting schemes using popular Bayesian model selection criteria such as the DIC and the LPML under the GWR context. The proposed methods are implemented in nimble. In our simulation studies, when there is no spatially varying covariate effect, the bandwidth is selected to give all observations close to uniform weight in estimating the coefficients for each individual location, whereas when there is indeed spatially varying covariate effect, the bandwidth is selected to achieve a balance between introducing bias for each location by taking into consideration nearby observations, and having too unstable estimates by placing the majority of emphasis on local observations and weighing down all others too heavily. The parameters estimated for each location have decent coverage rate that are close to the nominal 95% level.

Compared to the great circle distance, with a natural threshold of 1 to define “close enough,” the graph distance yields weighting schemes that produce models with robust parameter estimation and variable selection performance. Based on comparisons with classical frequentist GWR from the simulation studies, it is interesting for us to notice that the bandwidth selection results of the Bayesian approach are different from that of the frequentist approach. A partial reason for this pattern is that the bandwidth selection of frequentist approach is based on minimizing the summation of the sum of square errors (SSE) for each location. Our Bayesian approach, however, tries to maximize the whole data likelihood. Therefore, the frequentist approach tends to select smaller bandwidth than Bayesian approach.

A few issues beyond the scope of this paper are worth further investigation. In this work, we are only concerned with estimation of parameter for linear regression. Extension of similar ideas to generalized linear models, and semi-parametric models such as the Cox model, are worth developing. In the second alternative simulation scenario with regional variation patterns, both the frequentist and Bayesian GWR try to weigh neighbors as high as possible, leading to large bandwidths. Under the frequentist framework, clustering of covariate effects has been done using hierarchical clustering on the parameter estimation, which is ad hoc. Another approach is the penalized methods in Li and Sang (2019). In the Bayesian paradigm, however, hierarchical modeling provides an integrated framework that incorporates the latent cluster configuration layer. Development of such a framework is worth investigating. Also, we are assuming that a covariate is either in the true model for all locations, or not in the true model for all locations. There are cases where a covariate is important for some locations, but is minimally impacting for other locations. Identifying such locations is devoted to future research. Detecting a relationship between two areas that do not share a boundary (Gao and Bradley 2019) other than using graph distance is also an interesting future work.

Supplemental Material

Supplemental Material, gwrnimble_tobeongithub - Geographically Weighted Regression Analysis for Spatial Economics Data: A Bayesian Recourse

Supplemental Material, gwrnimble_tobeongithub for Geographically Weighted Regression Analysis for Spatial Economics Data: A Bayesian Recourse by Zhihua Ma, Yishu Xue and Guanyu Hu in International Regional Science Review

Supplemental Material

Supplemental Material, gwr_supp - Geographically Weighted Regression Analysis for Spatial Economics Data: A Bayesian Recourse

Supplemental Material, gwr_supp for Geographically Weighted Regression Analysis for Spatial Economics Data: A Bayesian Recourse by Zhihua Ma, Yishu Xue and Guanyu Hu in International Regional Science Review

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Dr. Ma’s research was supported by Project of Educational Commission of Guangdong Province of China #2019WQNCX104.

ORCID iDs

Yishu Xue

Guanyu Hu

Supplemental Material

Supplemental material for this article is available online.

References

Bhattacharyya

Bickel

P. J.

. 2014. “Community Detection in Networks Using Graph Distance.” arXiv preprint arXiv:1401.3915. https://arxiv.org/abs/1401.3915.

Boehm Vock

L. F.

Reich

B. J.

Fuentes

Dominici

. 2015. “Spatial Variable Selection Methods for Investigating Acute Health Effects of Fine Particulate Matter Components.” Biometrics 71 (1): 167–77.

Boscardin

W. J.

Gelman

. 1993. “Bayesian Computation for Parametric Models of Heteroscedasticity in the Linear Model.” In Advances in EconoMetrics, Vol. 11, edited by Fomby

T. B.

Hill

R. C.

, 87–109. Bingley, UK: Emerald Group Publishing Ltd.

Brunsdon

Fotheringham

A. S.

Charlton

M. E.

. 1996. “Geographically Weighted Regression: A Method for Exploring Spatial Nonstationarity.” Geographical Analysis 28 (4): 281–98.

Brunsdon

Fotheringham

A. S.

Charlton

M. E.

. 1998. “Geographically Weighted Regression.” Journal of the Royal Statistical Society: Series D (The Statistician) 47 (3): 431–43.

Brunsdon

Fotheringham

A. S.

Charlton

M. E.

. 2000. Geographically Weighted Regression as a Statistical Model. Technical report . Newcastle, UK: Department of Geography, University of Newcastle-upon-Tyne.

Cox

T. F.

Cox

M. A.

. 2000. Multidimensional Scaling. New York: Chapman and Hall/CRC.

da Silva

A. R.

Fotheringham

A. S.

. 2016. “The Multiple Testing Issue in Geographically Weighted Regression.” Geographical Analysis 48 (3): 233–47.

de Valpine

Turek

Paciorek

C. J.

Anderson-Bergman

Lang

D. T.

Bodik

. 2017. “Programming with Models: Writing Statistical Algorithms for General Model Structures with NIMBLE.” Journal of Computational and Graphical Statistics 26 (2): 403–13.

10.

Fotheringham

A. S.

Brunsdon

Charlton

M. E.

. 2002. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Chichester, UK: Wiley.

11.

Fotheringham

A. S.

Charlton

M. E.

Brunsdon

. 1998. “Geographically Weighted Regression: A Natural Evolution of the Expansion Method for Spatial Data Analysis.” Environment and Planning A 30 (11): 1905–27.

12.

Gao

Bradley

J. R.

. 2019. “Bayesian Analysis of Areal Data with Unknown Adjacencies Using the Stochastic Edge Mixed Effects Model.” Spatial Statistics 31:100357.

13.

Gelman

Stern

H. S.

Carlin

J. B.

Dunson

D. B.

Vehtari

Rubin

D. B.

. 2013. Bayesian Data Analysis. New York: Chapman and Hall/CRC.

14.

George

E. I.

McCulloch

R. E.

. 1993. “Variable Selection via Gibbs Sampling.” Journal of the American Statistical Association 88 (423): 881–89.

15.

Huffer

. 2020. “Modified Kaplan–Meier Estimator and Nelson–Aalen Estimator with Geographical Weighting for Survival Data.” Geographical Analysis 52 (1): 28–48.

16.

Ibrahim

J. G.

Chen

M.-H.

Sinha

. 2013. Bayesian Survival Analysis. New York: Springer Science & Business Media.

17.

Ishwaran

Rao

J. S.

. 2005. “Spike and Slab Variable Selection: Frequentist and Bayesian Strategies.” The Annals of Statistics 33 (2): 730–73.

18.

LeSage

J. P.

2004. “A Family of Geographically Weighted Regression Models.” In Advances in Spatial Econometrics. Methodology, Tools and Applications, edited by Anselin

Florax

R. J. G. M.

Rey

S. J.

, 241–64. Berlin, Heidelberg: Springer Berlin Heidelberg.

19.

Sang

. 2019. “Spatial Homogeneity Pursuit of Regression Coefficients for Large Datasets.” Journal of the American Statistical Association 114 (527): 1050–62.

20.

Chen

M.-H.

. 2020. “Bayesian Hierarchical Spatial Regression Models for Spatial Data in the Presence of Missing Covariates with Applications.” Applied Stochastic Models in Business and Industry. https://doi.org/10.1002/asmb.2568.

21.

Mei

C.-L.

Wang

Zhang

W.-X.

. 2006. “Testing the Importance of the Explanatory Variables in a Mixed Geographically Weighted Regression Model.” Environment and Planning A 38 (3): 587–98.

22.

Müller

Szymanski

Knop

Trinajstić

. 1987. “An Algorithm for Construction of the Molecular Distance Matrix.” Journal of Computational Chemistry 8 (2): 170–73.

23.

Páez

Uchida

Miyamoto

. 2002a. “A General Framework for Estimation and Inference of Geographically Weighted Regression Models: 1. Location-specific Kernel Bandwidths and a Test for Locational Heterogeneity.” Environment and Planning A 34 (4): 733–54.

24.

Páez

Uchida

Miyamoto

. 2002b. “A General Framework for Estimation and Inference of Geographically Weighted Regression Models: 2. Spatial Association and Model Specification Tests.” Environment and Planning A 34 (5): 883–904.

25.

Shao

1997. “An Asymptotic Theory for Linear Model Selection.” Statistica Sinica 7 (2): 221–64.

26.

Spiegelhalter

D. J.

Best

N. G.

Carlin

B. P.

Van Der Linde

. 2002. “Bayesian Measures of Model Complexity and Fit.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 64 (4): 583–639.

27.

Tobler

W. R.

1970. “A Computer Movie Simulating Urban Growth in the Detroit Region.” Economic Geography 46 (sup1): 234–40.

28.

Wang

Ying

Cheng

Wang

Liang

Huang

Gong

. 2012, August. “China’s Urban Expansion from 1990 to 2010 Determined with Satellite Remote Sensing.” Chinese Science Bulletin 57 (22): 2802–12.

29.

Wang

Mei

C.-L.

Yan

X.-D.

. 2008. “Local Linear Estimation of Spatially Varying Coefficient Models: An Improvement on the Geographically Weighted Regression Technique.” Environment and Planning A 40 (4): 986–1005.

30.

Wooldridge

J. M.

2015. Introductory Econometrics: A Modern Approach. Toronto, Canada: Nelson Education.

31.

Xue

Schifano

E. D.

. 2019. “Geographically Weighted Cox Regression for Prostate Cancer Survival Data in Louisiana.” Geographical Analysis. https://doi.org/10.1111/gean.12223.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

3.51 MB

4.72 MB

3.51 MB