Selection of the Bandwidth Matrix in Spatial Varying Coefficient Models to Detect Anisotropic Regression Relationships

Abstract

The commonly used Geographically Weighted Regression (GWR) fitting method for a spatial varying coefficient model is to select a bandwidth h for the geographic location (u, v), and assign the same weight to the two dimensions. However, spatial data usually present anisotropy. The introduction of a two-dimensional bandwidth matrix not only gives weight from two dimensions separately, but also increases the direction of kernel smoothness. The adaptive bandwidth matrix is more flexible. Therefore, in this paper, a two dimensional bandwidth matrix is introduced into the spatial varying coefficient model for parameter estimation. Through simulation experiments, the results obtained under the adaptive bandwidth matrix are compared with those obtained under the global bandwidth matrix, indicating the effectiveness of introducing the adaptive bandwidth matrix.

1. Introduction

The spatial varying coefficient model proposed by Brunsdon et al. (1996) [1] is a widely used and significant model in non-parametric statistics. It is a non-parametric smoothing method with the distance function between observation points as the weight.

The work of Brunsdon et al. (1998) [2] and Fotheringham et al. (1998) [3] shows the role of the Geographically Weighted Regression (GWR) method in studying the nonstationarity of spatial data regression relationships, and the effectiveness of the method was demonstrated through actual data. Since the two-dimensional region is larger than the one-dimensional region, it is easy to produce boundary effects. Wang et al. (2008) [4] extended the local linear fitting method from the varying coefficient model to the spatial varying coefficient model, and proposed the local linear GWR method. Recently, the backfitting technique has been successfully applied to the GWR models for dealing with the multiscale issue by manipulating the bandwidth in the kernel function (Fotheringham et al. (2017) [5], Leong and Yue (2017) [6], Yu et al. (2020) [7]). The kernel smoothing method is an effective method for solving non-parametric problems, in which the selection of the bandwidth h is more important than the selection of the kernel function itself. For the selection of bandwidth, the cross-validation method (CV) and its variants are widely used in practice, such as the generalized cross-validation method (GCV), or the AICc method to select the optimal bandwidth. If the estimated curve is smooth, the global bandwidth can be selected; while if the estimated curve has a complex structure, the adaptive bandwidth can be considered. There are two important studies on adaptive bandwidth: one is the K-nearest neighbor method proposed by Loftsgaarden and Quesenberry (1965) [8]; the second is the adaptive kernel density estimation proposed by Breiman et al. (1977) [9]. Subsequently, ref. [10] proposed using the square root method to select the adaptive bandwidth, that is, for any dimension, the bandwidth in the form $h_i\propto f{\left(x_i\right)}^{-1/2}$ can be selected. Fan and Gijbels (1995) [11] used a data-driven method to select the bandwidth; based on the estimated curvature, they proposed a new method for estimating the approximate deviation and variance of local polynomial fitting. Tian (2019) [12] used the cutting-edge knowledge of statistics to research the four major contemporary non-parametric statistics in the selection of bandwidth.

Many studies are limited to one-dimensional cases, while non-parametric problems will have more important applications in higher dimensions. The bandwidth matrix obtained under high-dimensional data not only affects the smoothness, but also affects the smoothing direction. For high-dimensional data, there are mainly the Plug-in insertion method, SCV and other methods for selecting the bandwidth matrix (Gramacki (2018) [13], Chacón and Duong (2018) [14]). Among them, the Plug-in insertion method is more stable for selecting bandwidth matrix H.

Wand and Jones (1994) [15] used the Plug-in method for multi-dimensional data to select a constrained bandwidth matrix. Duong and Hazelton (2003) [16] used the Plug-in method to select a binary unconstrained bandwidth matrix. Subsequently, Chacón and Duong (2010) [17] improved it by selecting unconstrained bandwidth matrix for high-dimensional data. In addition, the Bayesian method is also applied to the field of bandwidth selection. There are also many ways to select high-dimensional bandwidths by using the prior information from the data.

Brewer (2000) [18] proposed a new method for adaptive bandwidth in univariate kernel density estimation based on likelihood cross-validation and Bayesian graphical model analysis. Zhang et al. (2006) [19] took the bandwidth matrix elements as the parameters, and the posterior density was obtained by the Likelihood Cross-Validation criterion. The optimal window width was estimated by MCMC algorithm. Zougab et al. (2014) [20] used Bayesian estimation method to select adaptive bandwidth matrix, and so forth. For multi-dimensional data, besides bandwidth matrix, the product kernel is usually used to solve high-dimensional problems. Belaid et al. (2018) [21] applied the adaptive bandwidth selected by the Bayesian method to the high-dimensional discrete kernel. Gramacki and Gramacki (2017) [22] used the FFT method to solve some problems of the bandwidth matrix.

Sain et al. (1994) [23] used the CV method to select multiple bandwidths by the form of the product kernel. Wu et al. (2007) [24] regard the local bandwidth as the product of the local bandwidth factor and the global bandwidth, calculate the local window width factor with the clustering analysis algorithm, and obtain the multiple adaptive bandwidth. However, the kernel form of the multivariate product is essentially similar to the diagonal matrix. Xu (2009) [25] extended the bandwidth h of the coefficient function in the varying coefficient model to multidimensional, and theoretically tried to give a diagonal bandwidth matrix, where the values of diagonal elements were consistent. Wand and Jones (1993) [26] believe that the unconstrained bandwidth matrix is more general than the diagonal bandwidth matrix, and is applicable to a wider range. Unconstrained matrices can include many more smooth directions.

In many practical fields, such as geography, ecology, and meteorology, the data obtained usually contain geographic location information, and this kind of data is called spatial data. Spatial effects mainly come from spatial dependence (that is, spatial correlation) and spatial heterogeneity (that is, spatial non-stationarity). The spatial data are usually not in the same quantity distribution in all directions of geographical position, and the data distribution is not uniform, showing anisotropy. At present, there is no good way to deal with this phenomenon. The local linear geographically weight regression method is usually adopted to solve the fitting problem of the spatial variable coefficient model. In this method, the distance function between observation points is usually used as the weight value, and the same window width h is selected for the geographical position (u, v), that is, the data of two dimensions are given the same weight as isotropy. Considering the spatial anisotropy and dependence of the two variables (u, v), this paper introduces the two-dimensional Gaussian kernel function into the coefficient function estimation of the spatial variable coefficient model, and uses the binary unconstrained bandwidth matrix H to select different bandwidths from the two dimensions of geographical location, u and v. This can increase the direction of the smooth kernel, which can also explain the anisotropy. Furthermore, the adaptive bandwidth matrix is studied, and the coefficient estimates of the spatial variable coefficient models of the binary bandwidth matrix and the adaptive window width matrix are obtained respectively. Through simulation experiments, it is shown that the introduction of adaptive bandwidth matrix is more effective and flexible than the global bandwidth matrix.

2. Selection and Estimation of Bandwidth of Spatial Varying Coefficient Model

2.1. Spatial Varying Coefficient Model

The spatial varying coefficient model is to embed the data with the geographical location information into the varying coefficient model, which can not only avoid the dimensional disaster problem in fitting, but can also well reflect the spatial non-stationary and heterogeneous characteristics of the regression relationship. The specific form of the spatial varying coefficient model is as follows:

(1)$Y_i={\beta }_1(u_i,v_i)X_{i1}+\cdots +{\beta }_p(u_i,v_i)X_{ip}+{\epsilon }_i,i=1,2,\cdots ,n,$

where $(Y_i;X_{i1},X_{i2},\cdots ,X_{ip})$ represents the observed value of the dependent variable and independent variable at the geographic location$(u_i,v_i)$, $(i=1,2,\cdots ,n),{\beta }_j(u,v)(j=1,2,\cdots ,p)$, is a function of geographic location, and ${\epsilon }_i(i=1,2,\cdots ,n)$ is the $i.i.d.$ error term, which satisfies $E\left({\epsilon }_i\right)=0,Var\left({\epsilon }_i\right)={\sigma }^2$. The model allows a spatially varying interception by setting $x_{i1}=1,i=1,2,\cdots ,n$.

Brunsdo et al. [1] explicitly gave the specific form of the spatial varying coefficient model for the first time and proposed a method using the distance function between observation points as the weight, namely the geographically weighted regression method, which essentially belongs to the Nadaraya–Watson kernel estimation method. The boundary effects exists in this method, while the boundary of the two-dimensional geographic region is larger than that of the one-dimensional region, and the data of the boundary region is more prone to distortion. The local polynomial fitting method can automatically correct the boundary effect (Mei and Wang (2012) [27]). Wang et al. (2008) [4] extended the local linear fitting method from the varying coefficient model to the spatial varying coefficient model, and presented the locally linear GWR method. In order to avoid the boundary effects, the local linear GWR method is used for model fitting in this paper.

For a given location $(u_0,v_0)$ in the studied area, the estimates ${\widehat{\beta }}_j(u_0,v_0)$ of ${\beta }_j(u,v)(j=1,2,\cdots ,p)$ at $(u_0,v_0)$, and the matrix notation as follows:

$\widehat{\mathit{\beta }}(u_0,v_0)=({\widehat{\beta }}_1(u_0,v_0),{\widehat{\beta }}_2(u_0,v_0),\cdots ,{\widehat{\beta }}_p(u_0,v_0){)}^T,$

(2)$\widehat{\mathit{\beta }}(u_0,v_0)=({\mathbf{I}}_p,{\mathbf{0}}_p,{\mathbf{0}}_p){\left({\mathbf{X}}^T(u_0,v_0)\mathbf{W}(u_0,v_0)\mathbf{X}(u_0,v_0)\right)}^{-1}{\mathbf{X}}^T(u_0,v_0)\mathbf{W}(u_0,v_0)\mathbf{Y},$

where ${\mathbf{I}}_p$ is the p-order identity matrix and ${\mathbf{0}}_p$ is the p-order zero matrix, and $Y=(Y_1,Y_2,\cdots ,Y_n{)}^T$ is the observations vector of the response variable. The design matrix at point $(u_0,v_0)$ is:

$X(u_0,v_0)=\left(\begin{array}{ccccccccc}{X}_{11}& \cdots & X_{1p}& X_{11}(u_1-u_0)& \cdots & X_{1p}(u_1-u_0)& X_{11}(v_1-v_0)& \cdots & X_{1p}(v_1-v_0)\\ X_{21}& \cdots & X_{2p}& X_{21}(u_1-u_0)& \cdots & X_{2p}(u_1-u_0)& X_{21}(v_1-v_0)& \cdots & X_{2p}(v_1-v_0)\\ ⋮& & ⋮& ⋮& & ⋮& ⋮& & ⋮\\ X_{n1}& \cdots & X_{np}& X_{n1}(u_1-u_0)& \cdots & X_{np}(u_1-u_0)& X_{n1}(v_1-v_0)& \cdots & X_{np}(v_1-v_0)\end{array}\right),$

and $\mathbf{W}(u_0,v_0)=diag(w_1(u_0,v_0),w_2(u_0,v_0),$$\cdots ,w_n(u_0,v_0))$ is a diagonal weight matrix with $w_i(u_0,v_0)=K_h\left(d_{0i}\right)(i=1,2,\cdots ,n.)$, $K_h(·)$ is a kernel function. The kernel functions used in this paper are all Gaussian kernel functions. For the traditional GWR method, the specific form of the one-dimensional Gauss kernel function is as follows:

(3)$K_h\left(d_{ij}\right)=\frac{1}{h}K\left(\frac{d_{ij}}{h}\right)=\frac{1}{h}\frac{1}{\sqrt{2\pi }}exp(-\frac{1}{2}{\left(\frac{d_{ij}}{h}\right)}^2),j=1,2,\cdots ,n$

where, $d_{ij}$ represents the Euclidean distance between $(u_i,v_i)$ and $(u_j,v_j)(i=1.2,\cdots ,n).$ Specifically, the set of weights at $(u_0,v_0)$ in Equation (2) are $w_i(u_0,v_0)=K\left(\frac{d_{0i}}{h}\right),i=1,2,\cdots ,n.$

The specific form of the two-dimensional Gauss kernel function is as follows:

(4)$K_{\mathbf{H}}\left(\mathbf{t}\right)={\left|\mathbf{H}\right|}^{-\frac{1}{2}}K\left({\mathbf{H}}^{-\frac{1}{2}}\mathbf{t}\right)={\left|\mathbf{H}\right|}^{-\frac{1}{2}}\frac{1}{\sqrt[d]{2\pi }}exp(-\frac{1}{2}{(\mathbf{t}-{\mathbf{t}}_i)}^T{\mathbf{H}}^{-1}(\mathbf{t}-{\mathbf{t}}_i)).$

The set of weights at $(u_1,v_1)$ is ${\omega }_i(u_1,v_1)=K\left({\mathbf{H}}^{-\frac{1}{2}}({\mathbf{t}}_i-{\mathbf{t}}_1)\right),i=1,2,\cdots ,n.$ The spatial varying coefficient model involves two variables, u and v, so $d=2,t_i$ is a binary vector about u and v, $\mathbf{H}$ is a binary bandwidth matrix about $(u,v)$, the distance used here is formally known as the Mahalanobis distance.

The estimate of the dependent variable $\mathbf{Y}$ is:

(5)$\widehat{\mathbf{Y}}=\mathbf{LY},$

where $\mathbf{L}$ is the $p\times p$ hat matrix. The hat matrix used for parameter estimation using the one-dimensional kernel function has the following form:

(6)${\mathbf{L}}_1=\left(\begin{array}{c}({\mathbf{X}}_1^T,{\mathbf{0}}_{(1\times 2p)}){\left({\mathbf{X}}^T(u_1,v_1)\mathbf{W}(u_1,v_1)\mathbf{X}(u_1,v_1)\right)}^{-1}{\mathbf{X}}^T(u_1,v_1)\mathbf{W}(u_1,v_1)\\ ({\mathbf{X}}_2^T,{\mathbf{0}}_{(1\times 2p)}){\left({\mathbf{X}}^T(u_2,v_2)\mathbf{W}(u_2,v_2)\mathbf{X}(u_2,v_2)\right)}^{-1}{\mathbf{X}}^T(u_2,v_2)\mathbf{W}(u_2,v_2)\\ ⋮\\ ({\mathbf{X}}_n^T,{\mathbf{0}}_{(1\times 2p)}){\left({\mathbf{X}}^T(u_n,v_n)\mathbf{W}(u_n,v_n)\mathbf{X}(u_n,v_n)\right)}^{-1}{\mathbf{X}}^T(u_n,v_n)\mathbf{W}(u_n,v_n).\end{array}\right)$

The hat matrix used in parameter estimation using two-dimensional kernel function has the following form:

(7)${\mathbf{L}}_2=\left(\begin{array}{c}({\mathbf{X}}_1^T,{\mathbf{0}}_{(1\times 2p)}){\left({\mathbf{X}}^T(u_1,v_1){\mathbf{W}}_{\mathbf{H}}(u_1,v_1)\mathbf{X}(u_1,v_1)\right)}^{-1}{\mathbf{X}}^T(u_1,v_1)\mathbf{W}(u_1,v_1)\\ ({\mathbf{X}}_2^T,{\mathbf{0}}_{(1\times 2p)}){\left({\mathbf{X}}^T(u_2,v_2){\mathbf{W}}_{\mathbf{H}}(u_2,v_2)\mathbf{X}(u_2,v_2)\right)}^{-1}{\mathbf{X}}^T(u_2,v_2)\mathbf{W}(u_2,v_2)\\ ⋮\\ ({\mathbf{X}}_n^T,{\mathbf{0}}_{(1\times 2p)}){\left({\mathbf{X}}^T(u_n,v_n){\mathbf{W}}_{\mathbf{H}}(u_n,v_n)\mathbf{X}(u_n,v_n)\right)}^{-1}{\mathbf{X}}^T(u_n,v_n)\mathbf{W}(u_n,v_n).\end{array}\right)$

2.2. Selection of Bandwidth Matrix

(1) The bandwidth matrix of the unary variable is $\mathbf{H}=h^2$, Here, the square of the window width is used to consider the problem. In this paper, the generalized cross validation method (GCV) is used to select the one-dimensional bandwidth in traditional GWR method:

(8)$GCV\left(h\right)=\frac{1}{n}\sum _{i=1}^n{\left(\frac{Y_i-\widehat{Y_i}\left(h\right)}{1-\frac{1}{n}tr\left(\mathbf{L}\left(h\right)\right)}\right)}^2=\frac{n}{{(n-tr\left(\mathbf{L}\left(h\right)\right))}^2}\sum _{i=1}^n{(Y_i-\widehat{Y_i}\left(h\right))}^2,$

where L is the hat matrix, $tr(·)$ represents the trace of the matrix, $Y_i$ represents the value of the ith dependent variable, $\widehat{\mathbf{Y}}\left(h\right)=\mathbf{L}\left(h\right)\mathbf{Y}$. To select the optimal window width $h_0$, in practical application, a range is first given in accordance with a certain step size, the corresponding GCV(h) value is calculated, and the $h_0$ that makes GCV(H) reach the minimum is chosen, that is, $GCV\left(h_0\right)=\underset{h>0}{min}GCV\left(h\right)$, then $h_0$ is the selected optimal bandwidth.

(2) The form of the binary unconstrained window width matrix is as follows: $\mathbf{H}=\left[\begin{array}{cc}{h_1}^2& h_{12}\\ h_{12}& {h_2}^2\end{array}\right]$, H is a symmetric positive definite matrix, $h_1,h_2>0,\left|h_{12}\right|<h_1{h}_2$. For the selection of the binary unconstrained bandwidth matrix H, the Plug-in method proposed by Duong and Hazelton (2003) [16] is used for reference in this paper. The Plug-in method mainly relies on replacing the unknown function by the estimator in approximation, so that the asymptotic expression AMISE (Asymptotic Mean Integrated Squared Error) of the estimated MISE (Mean Integrated Squared Error) can reach the minimum value. Here, t can be regarded as a set of data about u and v, and f(t) is the target density function. The specific expressions of MISE and AMISE are as follows:

(9)$MISE\widehat{f}(·;\mathbf{H})=\mathbb{E}ISE\widehat{f}(·;\mathbf{H})=\mathbb{E}{\int }_{R^2}{[\widehat{f}(\mathit{x},\mathbf{H})-f\left(\mathit{x}\right)]}^{2}d\mathit{x}$

(10)$AMISE\widehat{f}(·;\mathbf{H})=n^{-1}{\left|\mathbf{H}\right|}^{-1/2}R\left(K\right)+\frac{1}{4}{\mu }_2{\left(K\right)}^2\left(vec{h}^T\mathbf{H}\right){\mathbf{\Psi }}_4\left(vech\mathbf{H}\right)$

(11)${\mathbf{H}}_{AMISE}=\underset{\mathbf{H}\in \mathbf{H}}{argmin}AMISE\widehat{f}(·;\mathbf{H}),$

where $R\left(K\right)={\int }_{R^2}K{\left(\mathit{x}\right)}^{2}d\mathit{x}$, ${\mu }_2\left(K\right)\mathbf{I}={\int }_{R^2}\mathit{x}{\mathit{x}}^{T}K\left(\mathit{x}\right)d\mathit{x}$, ${\mu }_2\left(K\right)<\infty$, $vech$ represents a special form of column vector, for example, $vech\left(\begin{array}{cc}a& c\\ c& d\end{array}\right)$ = $\left[\begin{array}{c}a\\ c\\ d\end{array}\right]$.

Then ${\mathbf{\Psi }}_4$ is a 3 × 3 matrix, the specific form is as follows:

(12)${\mathbf{\Psi }}_4={\int }_{R^2}vech[2D^{2}f\left(\mathit{x}\right)-dg\left[D^{2}f\left(\mathit{x}\right)\right]]vec{h}^T[2D^{2}f\left(\mathit{x}\right)-dg\left[D^{2}f\left(\mathit{x}\right)\right]]d\mathit{x}.$

$D^{2}f\left(\mathit{x}\right)$ is the Hessian matrix of f, and $dg\left[A\right]$ denotes the diagonal matrix formed by replacing all off-diagonal entries of A by zeroes. ${\mathbf{\Psi }}_4$ are integrals of products of density derivatives of $f(·)$.

To obtain the optimal bandwidth matrix H, the key is to estimate ${\mathbf{\Psi }}_4$. The elements in ${\mathbf{\Psi }}_4$ are ${\psi }_{\mathit{r}}={\int }_{R^2}{f}^{\left(\mathit{r}\right)}\left(\mathit{x}\right)f\left(\mathit{x}\right)d\mathit{x}$, $\mathit{r}=(r_1,r_2)$, $\left|\mathit{r}\right|=r_1+r_2$, $r_1$, $r_2$ are non-negative integers, $f^{\left(\right(\mathbf{r}\left)\right)}\left(\mathbf{t}\right)$ is the r-order partial derivative with respect to $(t_1,t_2)$. Firstly, the original data are standardized, and the test bandwidth matrix G is used to replace H to estimate each element ${\mathbf{\Psi }}_{\mathit{r}}$ in ${\mathbf{\Psi }}_4$, separately:

(13)${\widehat{\mathbf{\Psi }}}_{\mathit{r}}\left(\mathbf{G}\right)=n^{-1}\sum _{i=1}^n{\widehat{f}}^{\left(\mathit{r}\right)}({\mathbf{T}}_i;\mathbf{G})=n^{-2}\sum _{i=1}^n\sum _{j=1}^n{K}_{\mathbf{G}}^{\left(\mathit{r}\right)}({\mathbf{T}}_i-{\mathbf{T}}_j).$

Calculate the variance and deviation of related G. Let $\mathbf{G}=g^{2}I$. Since g is also unknown, calculate the variance and deviation of g to obtain the AMSE related to g:

(14)$AMSE{\widehat{\psi }}_{\mathit{r}}\left(g\right)=2n^{-2}{g}^{-2\left|\mathit{r}\right|-2}{\mathit{\psi }}_{\mathbf{0}}R\left(K^{\left(\mathit{r}\right)}\right)+{[n^{-1}{g}^{-\left|\mathit{r}\right|-2}{K}^{\left(\mathit{r}\right)}\left(\mathbf{0}\right)+\frac{1}{2}{g}^2{\mu }_2\left(K\right)\sum _{i=1}^2{\psi }_{\mathit{r}+2{\mathit{e}}_i}]}^2.$

${\mathit{e}}_i$ represents the unit vector whose ith value is 1, and sum over all $AMSE{\widehat{\psi }}_{\mathit{r}}\left(g\right)$ to get $SAMES{\mathbf{\Psi }}_4$.

(15)$SAMES{\mathbf{\Psi }}_4=\sum _{\mathit{r}:\left|\mathit{r}\right|=4}AMSE{\widehat{\psi }}_{\mathit{r}}\left(g\right).$

To take the g which minimizes $SAMES{\mathbf{\Psi }}_4$. Here, the final bandwidth matrix H is generally obtained by the two-stage method. The selection of the optimal bandwidth matrix can be implemented in the “ks” package of the R language.

(3) As for the selection of the adaptive bandwidth matrix, the Bayesian method proposed by Zougab et al. (2014) [20] was adopted in this paper. Under the Bayesian framework, the prior distribution is selected as the conjugate prior distribution, and the explicit expression of the adaptive window width matrix is obtained by combining with the square loss function. Here, t can be viewed as a set of data about u and v. The specific calculation steps are as follows: the kernel estimation is carried out by leave-one out method estimation.

${\widehat{f}}_{-i}\left({\mathbf{t}}_i\right|\left\{{\mathbf{t}}_{-i}\right\},{\mathbf{H}}_i)=\frac{1}{n-1}\sum _{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n{K}_{{\mathbf{H}}_i}({\mathbf{t}}_i-{\mathbf{t}}_j).$

The kernel function generally uses Gaussian kernel:

$\widehat{f}\left(\mathbf{t}\right)=\frac{1}{n}\sum _{i=1}^n\frac{1}{{\left(2\pi \right)}^{d/2}{\left|{\mathbf{H}}_i\right|}^{1/2}}exp\left\{-\frac{1}{2}{(\mathbf{t}-{\mathbf{t}}_i)}^T{{\mathbf{H}}_i}^{-1}(\mathbf{t}-{\mathbf{t}}_i)\right\}.$

At this time, its prior distribution can take its conjugate prior distribution, that is, the inverse Wishard prior distribution,

$\pi \left({\mathbf{H}}_i\right)=\frac{{\left|\mathbf{Q}\right|}^{r/2}{\left|{\mathbf{H}}_i\right|}^{-(r+d+1)/2}}{2^{rd/2}{\Gamma }_d(r/2)}exp\left\{-\frac{1}{2}tr\left(\mathbf{Q}{{\mathbf{H}}_i}^{-1}\right)\right\},$

where ${\Gamma }_d\left(q\right)=2^{dq}{\pi }^{d(d-1)/4}{\prod }_{k=0}^{d-1}\Gamma (q-\frac{k}{2})$ is the classic Gamma function.

Under the square loss function, the explicit expression of a posteriori estimate can be given, which can avoid the consumption of a lot of operation time during sampling calculation. For each bandwidth matrix ${\mathbf{H}}_i$, the form of a posteriori probability is as follows:

$\pi \left({\mathbf{H}}_i\right|{\mathbf{t}}_i)=\frac{\widehat{f}\left({\mathbf{t}}_i\right|\left\{{\mathbf{t}}_{-i}\right\},{\mathbf{H}}_i)\pi \left({\mathbf{H}}_i\right)}{{\int }_{\mathbf{H}}\widehat{f}\left({\mathbf{t}}_i\right|\left\{{\mathbf{t}}_{-i}\right\},{\mathbf{H}}_i)\pi \left({\mathbf{H}}_i\right)d{\mathbf{H}}_i}.$

The posterior distribution $\pi \left({\mathbf{H}}_i\right|{\mathit{x}}_i)$ can be expressed as follows:

$\pi \left({\mathbf{H}}_i\right|{\mathbf{t}}_i)=\frac{{\sum }_{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n{\left|{\mathbf{H}}_i\right|}^{-(r+d+2)/2}exp\left\{-\frac{1}{2}tr\left[\left\{({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T+\mathbf{Q}\right\}{{\mathbf{H}}_i}^{-1}\right]\right\}}{2^{(r+1)d/2}{\Gamma }_d((r+1)/2){\sum }_{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n{|({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T+\mathbf{Q}|}^{-(r+1)/2}}.$

For each j term in the posterior distribution, where $j\ne i$ is multiplied by something with unit 1, where the quantity with unit 1 can be expressed as:

$m=\frac{|({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T{+\mathbf{Q}|}^{-(r+1)/2}}{|({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T{+\mathbf{Q}|}^{-(r+1)/2}}.$

Then, we can get:

$\begin{array}{cc}\hfill & \pi \left({\mathbf{H}}_i\right|{\mathbf{t}}_i)\hfill \\ \hfill =& \frac{{\sum }_{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n{\left|{\mathbf{H}}_i\right|}^{-(r+d+2)/2}exp\left\{-\frac{1}{2}tr\left[\left\{({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T+\mathbf{Q}\right\}{{\mathbf{H}}_i}^{-1}\right]\right\}\times m}{2^{(r+1)d/2}{\Gamma }_d((r+1)/2){\sum }_{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n{|({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T+\mathbf{Q}|}^{-(r+1)/2}}\hfill \\ \hfill =& \sum _{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n\frac{|({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T{+\mathbf{Q}|}^{-(r+1)/2}}{{\sum }_{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n{|({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T+\mathbf{Q}|}^{-(r+1)/2}}\hfill \\ \hfill & \times \iota \omega (r+1,({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T+\mathbf{Q})\hfill \end{array}$

(16)${\mathbf{H}}_i=\frac{1}{r-d}\sum _{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n\frac{|({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T{+\mathbf{Q}|}^{-(r+1)/2}}{{\sum }_{\genfrac{}{}{0pt}{}{j=1}{j\ne i}}^n{|({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T+\mathbf{Q}|}^{-(r+1)/2}}\times \left\{({\mathbf{t}}_i-{\mathbf{t}}_j){({\mathbf{t}}_i-{\mathbf{t}}_j)}^T+\mathbf{Q}\right\},$

where r is the degree of freedom, given $r=n^{2/(d+4)}+d$, Q is the covariance matrix of the sample. For the adaptive bandwidth matrix, the specific form of the two-dimensional Gauss kernel function is as follows:

(17)$K_{{\mathbf{H}}_i}\left(\mathbf{t}\right)=|{\mathbf{H}}_i{|}^{-\frac{1}{2}}K\left({\mathbf{H}}_i^{-\frac{1}{2}}\mathbf{t}\right)={\left|{\mathbf{H}}_i\right|}^{-\frac{1}{2}}\frac{1}{\sqrt[d]{2\pi }}exp(-\frac{1}{2}{(\mathbf{t}-{\mathbf{t}}_i)}^T{\mathbf{H}}_i^{-1}(\mathbf{t}-{\mathbf{t}}_i)).$

The set of weights at $(u_1,v_1)$ is ${\omega }_i(u_1,v_1)=K\left({{\mathbf{H}}_1}^{-\frac{1}{2}}({\mathbf{t}}_i-{\mathbf{t}}_1)\right),i=1,2,\cdots ,n,$ here $d=2,{\mathbf{t}}_i$ is a binary vector about $(u,v),{\mathbf{H}}_i$ is a binary bandwidth matrix about $(u,v)$. The corresponding hat matrix L has the following forms:

(18)${\mathbf{L}}_3=\left(\begin{array}{c}({\mathbf{X}}_1^T,{\mathbf{0}}_{(1\times 2p)}){\left({\mathbf{X}}^T(u_1,v_1){\mathbf{W}}_{{\mathbf{H}}_1}(u_1,v_1)\mathbf{X}(u_1,v_1)\right)}^{-1}{\mathbf{X}}^T(u_1,v_1)\mathbf{W}(u_1,v_1)\\ ({\mathbf{X}}_2^T,{\mathbf{0}}_{(1\times 2p)}){\left({\mathbf{X}}^T(u_2,v_2){\mathbf{W}}_{{\mathbf{H}}_2}(u_2,v_2)\mathbf{X}(u_2,v_2)\right)}^{-1}{\mathbf{X}}^T(u_2,v_2)\mathbf{W}(u_2,v_2)\\ ⋮\\ ({\mathbf{X}}_n^T,{\mathbf{0}}_{(1\times 2p)}){\left({\mathbf{X}}^T(u_n,v_n){\mathbf{W}}_{{\mathbf{H}}_n}(u_n,v_n)\mathbf{X}(u_n,v_n)\right)}^{-1}{\mathbf{X}}^T(u_n,v_n)\mathbf{W}(u_n,v_n).\end{array}\right)$

3. Simulation Experiment

In order to reflect the spatial scale difference and interdependence between $(u,v)$, both $u_1$ and $v_1$ are taken as n sample points on uniformly distributed of $U(10,20)$, and set the matrix $A=\left[\begin{array}{cc}1& 2\\ 3& 4\end{array}\right]$ and $[u,v]=[u_1,v_1]*A$. By rotating $(u_1,v_1)$, a set of geographic locations $(u,v)$ with relevant relationships is obtained. The rotated $(u,v)$ was used in the experiment to reflect the spatial scale difference and the dependent relationship between the geographical positions of the spatial data.

The model used to generate spatial data is as follows:

(19)$Y_i={\beta }_1(u_i,v_i)X_{i1}+{\beta }_2(u_i,v_i)X_{i2}+{ϵ}_i,i=1,2,\cdots ,n,$

where the independent variable $X_1$ is a column vector of all 1, that is, the observation values $X_{11},X_{21},\cdots ,X_{n1}$ are all 1, and the observation values $X_{12},X_{22},\cdots ,X_{n2}$ are independent values generated from the uniform distribution U(0,2), the error ${ϵ}_1,{ϵ}_2,\cdots ,{ϵ}_n$ are n independent and identically distributed random numbers that obey the normal distribution $N(0,2^2)$, where the sample size is taken as $n=100$, let ${\beta }_1(u,v)=\frac{1}{6}(u+v),{\beta }_2(u,v)=\frac{1}{3}u$, and the real surfaces are shown in Figure 1.

The coefficient functions in Equation (19) are calibrated by the Formula (2), and the optimal bandwidth in the kernel function is selected in three ways: the GCV-based method in Equation (8), the Plug-in-based matrix H method in Equation (11) and the adaptive bandwidth matrix H method in Equation (16). We get the GCV-based optimal bandwidth $h_0=35.09$, and the Plug-in-based optimal bandwidth matrix H is ${\mathbf{H}}_0=\left[\begin{array}{cc}7.9068& 17.1485\\ 17.1485& 38.2639\end{array}\right]$. The accuracy of parameter estimation is evaluated by the mean square error (MSE). Its definition is $MSE\left(\widehat{\theta }\right)=\frac{1}{n}{\sum }_{i=1}^n{(\theta -\widehat{\theta })}^2$; the smaller the value of the $MSE$, the better the fitting effect. In order to compare the performance of the above three estimation methods, keep other quantities fixed, re-sample data with a sample size of n from $N(0,2^2)$, do $N=200$ repeated experiments, calculate the $MSE$ of $\mathbf{Y}$, ${\beta }_1,{\beta }_2$, and take the average of the N replications. The numerical results are shown in the Table 1. Additionally, Figure 2 and Figure 3 display the estimated coefficients and the estimation residuals of $\widehat{{\beta }_1}(u,v)$ and $\widehat{{\beta }_2}(u,v)$ with the GCV-based bandwidth, the bandwidth matrix and the adaptive bandwidth matrix method, respectively.

It can be found from the simulation results in Table 1 that the MSE of regression coefficients estimation obtained by using the bandwidth matrix H and the adaptive bandwidth matrix are smaller than that obtained by using GCV-based bandwidth h for 200 replications. The parameter estimation MSE of ${\beta }_1(u,v)$ and ${\beta }_2(u,v)$ obtained by the adaptive bandwidth matrix is less than that of the parameter estimation of ${\beta }_1(u,v)$ and ${\beta }_2(u,v)$ obtained by the direct use of the bandwidth matrix. In other words, the estimation result obtained by the adaptive bandwidth matrix is better than that obtained by the use of the bandwidth matrix. That means the performance of two dimensions’ unconstrained bandwidth matrix for geographical location $(u,v)$ is better than that of selecting one bandwidth for two variables $(u,v)$ in the spatially varying coefficient model. In addition, when the bandwidth matrix H obtained by the Plug-in method is used for parameter estimation, the running time of the program is shorter than that obtained by the GCV-based method.

From the residual graph of the coefficient function in Figure 2 and Figure 3, the estimated residuals are dispersed by selecting the same bandwidth h for u and v, and the estimated residuals are relatively concentrated around zero by selecting the bandwidth matrix H, although there are a few points with large differences. The residuals boxplots of ${\beta }_1(u,v)$ in Figure 2 and ${\beta }_2(u,v)$ in Figure 3 show that the estimated residuals of $\beta (u,v)$ are more concentrated by using the adaptive bandwidth matrix than the global bandwidth matrix. The reason for the occurrence of large bias points may be related to the non-diagonal elements of the bandwidth matrix H, or an adaptive bandwidth matrix can be introduced. The introduction of the adaptive bandwidth matrix can reduce the outliers. It further proves the effectiveness of introducing a two-dimensional adaptive bandwidth matrix into the spatial varying coefficient model.

The above results show that, in the spatial varying coefficient model, considering the spatial anisotropy and correlation of spatial coordinates $(u,v)$, the estimated coefficient surfaces obtained by selecting bandwidth matrix H is more accurate than that obtained by selecting the same bandwidth h for the two spatial coordinates $(u,v)$. In addition, compared with GCV-based method, using the Plug-in-based method to select the bandwidth matrix can reduce the number of iterations in model fitting and improve computational efficiency.

4. Empirical Research

The Fractional Vegetation Cover (FVC) quantifies the denseness of vegetation and reflects the growth trend of vegetation. FVC plays an important role in reflecting the process of climate change and is an important indicator of ecological environmental change and is a useful parameter for many environmental and climate-related applications. In this research, the datasets of variables FVC, temperature and precipitation in Yili, Xinjiang were collected to study the spatial varying relationships between FVC and climate factors.

4.1. Data Sources

Yili Kazakh Autonomous Prefecture was the selected research area, located in northwest China’s Xinjiang Uygur Autonomous Region, surrounded by mountains, and vulnerable to influence of airflow around the Tianshan mountains, the Wusun mountain, Kequqin and the Nalati mountain. Inside the Yili river near the Yili river valley, there is a humid climate and abundant rainfall; it is arguably the wettest region in Xinjiang. In this paper, the data on air temperature, precipitation and FVC in Yili, Xinjiang in 2011 adopted by Feng et al. (2016) [28] were used, with a total of 835 corresponding data points.

4.2. Modelling

The specific form of the spatial varying coefficient model in this section is as follows:

(20)$Y_i={\beta }_1(u_i,v_i)X_{i1}+{\beta }_2(u_i,v_i)X_{i2}+{\beta }_3(u_i,v_i)X_{i3}+{ϵ}_i,$

where $Y_i$ represents the observation value of the ith vegetation coverage index, $X_{i1}$ is always equal to 1, $X_{i2}$ represents the observation value of the ith annual average precipitation, and $X_{i3}$ represents the observation value of the ith annual average temperature, ${\beta }_1(u_i,v_i)$ represents the benchmark vegetation coverage index, ${\beta }_2(u_i,v_i)$ represents the influence intensity of the annual average precipitation at the ith point on the vegetation coverage index, ${\beta }_3(u_i,v_i)$ represents the intensity of influence of the annual average temperature at the ith point on the vegetation coverage index, ${ϵ}_i$ represents the residual. For the estimation of the coefficient function, the local linear GWR method is used for model fitting and parameter estimation.

4.3. Results and Analysis

4.3.1. Results of Unary Bandwidth and Bandwidth Matrix

According to the method of selecting the bandwidth matrix introduced in Section 2.2, the bandwidth h = 12,938.16 and the bandwidth matrix $\mathbf{H}=\left(\begin{array}{cc}614198610& 2083252\\ 2083252& 253867658\end{array}\right)$, according to formulas in Equations (3) and (4), a one-dimensional kernel function and a two-dimensional kernel function, are constructed to fit the model respectively and to obtain the MSE of the estimated value of the dependent variable and the running time of the program. The results in Table 2 show that the MSE(Y) obtained by the two-dimensional kernel function with the bandwidth matrix is smaller, the estimation is more accurate, the program running time is greatly shortened, and the running efficiency is greatly improved.

The estimated values obtained from the bandwidth matrix H in Table 3 show that: the coefficient function ${\beta }_2(u,v)$ and the coefficient function ${\beta }_3(u,v)$ of annual mean temperature and precipitation both show obvious spatial and regional differences. There was a positive correlation between precipitation and FVC, and the regions with more precipitation had higher vegetation coverage. There was also a positive correlation between temperature and FVC.

The results in Figure 4 show that, in the spatial varying coefficient model, the estimated value obtained by selecting the binary bandwidth matrix H for the kernel function is roughly consistent with the real vegetation coverage index, and the residual is small. The vegetation coverage index in most areas of Yili is relatively high, while the vegetation coverage index in the border area and the northwest is slightly lower, which is consistent with the actual situation. The reasons for this phenomenon may be related to temperature, precipitation, terrain, and so forth. It should be noted that, in addition to intuitively exploring spatial varying regression relationships by summary tables or images of the estimated coefficient functions, the statistical inference methods can be used to test the significance of spatial non-stationary relationships, including the residual sum of square based statistics with an F-distribution (Leung et al., 2000 [29]; Mei et al., 2006 [30]) and the bootstrap methods (Wei and Qi 2012 [31]; Mei et al., 2016 [32]).

4.3.2. Results of Bandwidth Matrix and Adaptive Bandwidth Matrix

According to the adaptive bandwidth matrix H method discussed in Section 2, Equation (16), the estimated result of the coefficient function can be found in Table 4, and Figure 5 shows the prediction of the FVC and the residuals with the adaptive bandwidth matrix.

By observing the estimated values of the coefficient functions in Table 4, it is easy to see that the values of $\beta (u,v)$ are mostly positive numbers. For ${\beta }_2(u,v)$ and ${\beta }_3(u,v)$, most of the values are positive, that is, precipitation and temperature have a positive effect on the FVC. By comparing Table 3 and Table 4 we can find that, compared to the range of parameter estimates obtained with a fixed bandwidth matrix, the range of parameter estimates obtained by the adaptive bandwidth matrix is smaller.

5. Summary and Prospect

(1) When there is a scale problem and a correlation between two variables, that is, the smoothness between the two variables is different and there is anisotropy in space, it is more reasonable to consider the two variables in selecting the binary bandwidth matrix than to select the same bandwidth fitting value between the two variables. At the same time, as an intermediary between a one-dimensional univariate and a high-dimensional multivariate, a bivariate has a smoother direction than a univariate kernel, and is easy to visualize. This also makes it possible to generalize the spatio-temporal variable coefficient model to the unconstrained bandwidth matrix of the spatio-temporal variable coefficient model or an even higher dimensional bandwidth matrix.

(2) Compared with the GCV method, using the Plug-in method to select the bandwidth matrix can reduce the number of iterations in the process of selecting the optimal bandwidth h, and can shorten the iteration time of program operation and improve the efficiency of calculation.

(3) By observing the residuals of the parameter estimates, we find that the residual range of the parameter estimates obtained under the bandwidth is small, while most of the residual values of the parameter estimates obtained under the bandwidth matrix are very concentrated. It will have points with larger residual values, which may be related to non-diagonal elements. We can also try to add an adaptive bandwidth matrix, which is worth continuing to study.

Author Contributions

Conceptualization, X.H. and H.Z.; methodology, Y.L.; software, Y.L.; validation, X.H. and H.Z.; formal analysis, H.Z. and Y.L.; investigation, Y.L. and X.H.; resources, Q.S.; data curation, Q.S.; writing—original draft preparation, Y.L.; writing—review and editing, H.Z.; visualization, Y.L. and H.Z.; supervision, X.H. and H.Z.; project administration, H.J.; funding acquisition, X.H. All authors have read and agreed to the published version of the manuscript.

Funding

This study is supported by the National Natural Science Foundation of China grand number 11961065 and U1703237). Zhang’s work is also partially supported by the Natural Science Foundation of Xinjiang grant number 2019D01C045 and The Ministry of education of Humanities and Social Science project grant number 19YJA910007.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data in the paper can be requested by email to the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The real surface of the coefficient function β1(u,v) (left) and β2(u,v) (right).

Figure 2. The estimated surface of β1(u,v) (top row) and residuals boxplots (bottom row) with the GCV−based method (left panel), the bandwidth matrix H (middle panel) and the adaptive bandwidth matrix H method (right panel).

Figure 3. The estimated surface of β2(u,v) (top row) and residuals boxplots (bottom row) with the GCV−based method (left panel), the bandwidth matrix H (middle panel) and the adaptive bandwidth matrix H method (right panel).

Figure 4. Contour oftrue FVC in the research area (top row); The estimated contour and the residual of the FVC with bandwidth matrix H (bottom row).

Figure 5. The estimation of the FVC with the adaptive bandwidth matrix H (left) and the residuals (right).

Table 1

MSE results of the GCV-based bandwidth h, the bandwidth matrix H and the adaptive bandwidth matrix for N = 200 replications.

Bandwidth Type	GCV-Based	Matrix H-Based	Adaptive Matrix H-Based
MSE(Y)	3.0298	2.3928	0.6633
MSE(${\beta }_1$)	0.0048	0.0136	0.0011
MSE(${\beta }_2$)	0.0105	0.0137	0.0012
Running time	14.24 s	2.38 s	5.67 s

Table 2

Selected bandwidth and results of MSE.

Bandwidth Type	Bandwidth h	Bandwidth Matrix H
MSE(Y)	0.074	0.039
Running time	33.69 min	1.32 min

Table 3

Regression coefficients estimation with bandwidth matrix H.

Regression Coefficients	Minimum	1/4 Quantile	Median	3/4 Quantile	Maximum
${\beta }_1(u,v)$	−10.5274	−0.0085	0.8041	1.3285	12.5262
${\beta }_2(u,v)$	−0.0349	−0.0017	0.0003	0.0021	0.0358
${\beta }_3(u,v)$	−0.0788	0.0013	0.0155	0.0566	0.3834

Table 4

Regression coefficients estimation with the adaptive bandwidth matrix H.

Regression Coefficients	Minimum	1/4 Quantile	Median	3/4 Quantile	Maximum
${\beta }_1(u,v)$	−5.3196	0.3858	0.8049	1.1690	6.6727
${\beta }_2(u,v)$	−0.0207	−0.0011	0.0002	0.0013	0.0196
${\beta }_3(u,v)$	−0.0288	−0.0008	0.0141	0.0520	0.2599