# Bayesian Estimation of Partially Linear Additive Spatial Autoregressive Models with P-Splines

Authors: Zhiyong Chen; Minghui Chen; Guodong Xing
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, we aim to develop a partially linear additive spatial autoregressive model (PLASARM), which is a generalization of the partially linear additive model and spatial autoregressive model. It can be used to simultaneously evaluate the linear and nonlinear effects of the covariates on the response for spatial data. To estimate the unknown parameters and approximate nonparametric functions by Bayesian P-splines, we develop a Bayesian Markov Chain Monte Carlo approach to estimate the PLASARM and design a Gibbs sampler to explore the joint posterior distributions of unknown parameters. Furthermore, we illustrate the performance of the proposed model and estimation method by a simulation study and analysis of Chinese housing price data.


---

## Body

## 1. Introduction

Spatial econometrics is a subfield of econometrics dealing with spatial interaction effects among geographical units Elhorst [1]; it has been widely applied in many research fields such as economics, geography, and environmental science. The traditional spatial parametric autoregressive models have received a lot of attention from theoretical econometricians and applied researchers, for example, Anselin [2], Case [4], Cressie [3], LeSage [4], LeSage [5], Anselin and Bera [6], Lee [7], Lee [8], Kakamu and Wago [9], LeSage and Pace [10], Kazar and Celik [11], and Piribauer and Crespo Cuaresma [12] among others. However, the spatial parametric autoregressive models are highly sensitive to model misspecification; they may not be adequate in many complex situations. Indeed, it is confirmed that some economic relationships in space exhibit highly nonlinear shapes [13–16]. Neglecting the underlying nonlinear relationship in spatial parametric autoregressive models frequently leads to inconsistent estimation of the parameters.Although spatial nonparametric autoregressive models can be used to improve the performance of the spatial parametric autoregressive models, it is unavoidable to bear the risk of misspecifying the link function. However, when the dimension of the explanatory variables is high, a fully nonparametric spatial autoregressive model is too hard to be explored because of the so-called “curse of dimensionality” [17]. Several nonparametric dimension-reduction techniques have been developed to overcome the problem, for instance, single-index model [18, 19], additive model [20], and varying-coefficient model [21], among others. Thus, semiparametric spatial autoregressive models have been proposed for dealing with spatial data. To capture the underlying relationships between the response variables and their associated covariates, semiparametric spatial autoregressive models have gained much attention in the literature of statistics and econometrics. For example, Su and Jin [22] discussed the quasi-likelihood estimator for the semiparametric partially linear spatial autoregressive model; Su [23] proposed semiparametric generalized method of moment (GMM) estimation of the semiparametric spatial autoregressive model; Sun et al. [24] developed a profile likelihood estimator for the semiparametric spatial dynamic model; Chen et al. [25] presented a two-step Bayesian approach for the semiparametric spatial autoregressive model; Wei and Sun [26] considered the semiparametric GMM estimation of a spatial model with space-varying coefficients; Hoshino [27] proposed a semiparametric series GMM estimator for the semiparametric spatial autoregressive model; Krisztin [28] studied a novel Bayesian semiparametric estimation for the penalized spline spatial autoregressive model; and Wei et al. [29] considered profile quasi-maximum likelihood method to estimate the partially linear varying-coefficient spatial autoregressive model. Sun and Wu [30] discussed the GMM to estimate the partially linear single-index spatial autoregressive model. Cheng et al. [31] developed the GMM estimation of the partially linear single-index spatial autoregressive model. Du et al. [32] investigated a GMM estimator for PLASARM. Hu et al. [33] presented a new partial functional linear spatial autoregressive model.In recent years, semiparametric models have become an important research topic between statistics and econometrics due to the explanation of parametric modeling and the flexibility of nonparametric modeling. Among the various semiparametric models, the most popular ones are perhaps the partially linear additive models. They contain both linear and nonlinear additive components. As they can not only provide more flexible models than the stringent linear models but also mitigate the “curse of dimensionality” phenomenon encountered in nonparametric regression. Therefore, the partially linear additive models provide a good balance between the flexibility of the additive models and the interpretation of the partially linear models. A lot of authors have developed various methods to analyze such models: local linear method [34], spline estimation [35–37], variable selection [38–41], quantile regression [42–45], and so on. But most of this work has not been applied to analyze spatial data. It is a distinctive challenge facing analysts of spatial data to characterize certain flexible functional forms that try to account for the latent nonlinearity of PLASARM. Combining the partially linear additive models with the spatial autoregressive models, we develop a Bayesian P-splines method and design a Gibbs sampler to explore the joint posterior distributions, along with a Markov chain Monte Carlo algorithm, to estimate the unknown parameters and approximate the nonparametric functions of PLASARM.The rest of the article proceeds as follows: In Section2, we present PLASARM for spatial dependent data and discuss its identifiability conditions and then obtain the likelihood functions by approximating the link function with a Bayesian P-spline method. In Section 3, we specify the prior distributions, derive the full conditional posteriors of the unknown parameters, and describe the detailed sampling algorithms in order to provide a Bayesian P-splines method for analysis. The applicability and practicality of the proposed model and estimation method are demonstrated through a simulation study and a real dataset in Section 4. In Section 5, we conclude the paper with a summary.


## 2. Model and Likelihood

### 2.1. Model

The partially linear additive spatial autoregressive model is defined as follows:(1)yi=ρ∑l=1nwilyl+xiTα+∑j=1pgjzij+εi,i=1,…,n,where yi is the i-th observation of response variable, xi=xi1,…,xiqT and zi=zi1,…,zipT are the i-th observations of q and p dimensional covariate vectors, respectively, wil is a specified constant spatial weight, gj⋅j=1,…,p are unknown univariate nonparametric functions, ρ is an unknown spatial parameter reflecting spatial autocorrelation between neighbors with stability condition ρ<1, α=α1,…,αqT is a q×1 vector of the unknown parameters, and εi’s are mutually independent and identically distributed normal random variables with zero mean and variance σ2. To ensure model identifiability of the nonparametric functions, it is often assumed that the condition ∑i=1ngjzij=0 for j=1,…,p.


### 2.2. Bayesian P-Splines and Likelihood

In analyzing the proposed PLASARMs defined by (1), modeling the smooth nonparametric function is an important issue. Because of the advantages of the Bayesian approach and the nice features of Bayesian P-splines [46], we will concentrate on applying the Bayesian P-splines approach. We intend to approximate unknown function gj⋅ in (1) by B-splines [47]. For j=1,…,p, the unknown function gj⋅ is a polynomial spline of degree mj with kj order interior knots ξj=ξj1,…,ξjkjT with aj<ξj1<⋯<ξjkj<bj, that is,(2)gjzj=∑l=1KjBjlzjβjl=BjTzjβj,zj∈aj,bj,where Kj=1+mj+kj is the number of splines determined by the number of knots, Bjzj=Bj1zj,…,BjKjzjT is a Kj×1 vector of spline basis that is determined by the knot vector ξj, βj=βj1,…,βjKjT is a Kj×1 vector of spline coefficients, and(3)aj=min1⩽i⩽nzijandbj=max1⩽i⩽nzij.In practice, a choice ofKj in the range of 10 to 30 provides the flexibility of fitting. To prevent overfitting due to using a relatively large number of knots, Eilers and Marx [48] proposed applying a difference penalty on coefficients of adjacent B-splines. Let 1n=1,…,1T and Bj be an n×Kj matrix with BjTzij as its i-th row. To achieve identification, we set ∑i=1n∑l=1KjBjlzijβjl=0, which is equivalent to 1nTBjβj=0. Denote Qj=1nTBj; then, the constraint becomes Qjβj=0.It follows from (2) that model (1) can be represented as(4)yi=ρ∑l=1nwilyl+xiTα+∑j=1pBjTzjβj+εi=ρ∑l=1nwilyl+xiTα+BTziβ+εi,i=1,…,n,where BTzi=B1Tzi1,…,BpTzipT and β=β1T,…,βpTT. Then, the matrix form of the model (1) can be written as(5)y=ρWy+xTα+BTzβ+ε,where y=y1,…,ynT, x=x1,…,xnT, ε=ε1,…,εnT, W=wil is an n×n specified constant spatial weight matrix, and BTz is an n×K1+⋯+Kp matrix with BTzi as its i-th row.The likelihood function corresponding to (5) is as follows:(6)Lα,β,σ2,ρ|y,x,z∝σ−nIn−ρWexp−12σ2y−ρWy−xTα−BTzβTy−ρWy−xTα−BTzβ=σ−nAρexp−12σ2Aρy−xTα−BTzβTAρy−xTα−BTzβ≐σ−nAρexp−12σ2Aρy−Bx,zθTAρy−Bx,zθ,where x=x1,…,xnT, z=z1,…,znT, K=∑j=1pKj, θ=αT,βTT is a K+q×1 vector of the regression coefficient, Bx,z=x,BTz is an n×K+q matrix, In is an identity matrix of order n, and Aρ=In−ρW.


## 2.1. Model

The partially linear additive spatial autoregressive model is defined as follows:(1)yi=ρ∑l=1nwilyl+xiTα+∑j=1pgjzij+εi,i=1,…,n,where yi is the i-th observation of response variable, xi=xi1,…,xiqT and zi=zi1,…,zipT are the i-th observations of q and p dimensional covariate vectors, respectively, wil is a specified constant spatial weight, gj⋅j=1,…,p are unknown univariate nonparametric functions, ρ is an unknown spatial parameter reflecting spatial autocorrelation between neighbors with stability condition ρ<1, α=α1,…,αqT is a q×1 vector of the unknown parameters, and εi’s are mutually independent and identically distributed normal random variables with zero mean and variance σ2. To ensure model identifiability of the nonparametric functions, it is often assumed that the condition ∑i=1ngjzij=0 for j=1,…,p.


## 2.2. Bayesian P-Splines and Likelihood

In analyzing the proposed PLASARMs defined by (1), modeling the smooth nonparametric function is an important issue. Because of the advantages of the Bayesian approach and the nice features of Bayesian P-splines [46], we will concentrate on applying the Bayesian P-splines approach. We intend to approximate unknown function gj⋅ in (1) by B-splines [47]. For j=1,…,p, the unknown function gj⋅ is a polynomial spline of degree mj with kj order interior knots ξj=ξj1,…,ξjkjT with aj<ξj1<⋯<ξjkj<bj, that is,(2)gjzj=∑l=1KjBjlzjβjl=BjTzjβj,zj∈aj,bj,where Kj=1+mj+kj is the number of splines determined by the number of knots, Bjzj=Bj1zj,…,BjKjzjT is a Kj×1 vector of spline basis that is determined by the knot vector ξj, βj=βj1,…,βjKjT is a Kj×1 vector of spline coefficients, and(3)aj=min1⩽i⩽nzijandbj=max1⩽i⩽nzij.In practice, a choice ofKj in the range of 10 to 30 provides the flexibility of fitting. To prevent overfitting due to using a relatively large number of knots, Eilers and Marx [48] proposed applying a difference penalty on coefficients of adjacent B-splines. Let 1n=1,…,1T and Bj be an n×Kj matrix with BjTzij as its i-th row. To achieve identification, we set ∑i=1n∑l=1KjBjlzijβjl=0, which is equivalent to 1nTBjβj=0. Denote Qj=1nTBj; then, the constraint becomes Qjβj=0.It follows from (2) that model (1) can be represented as(4)yi=ρ∑l=1nwilyl+xiTα+∑j=1pBjTzjβj+εi=ρ∑l=1nwilyl+xiTα+BTziβ+εi,i=1,…,n,where BTzi=B1Tzi1,…,BpTzipT and β=β1T,…,βpTT. Then, the matrix form of the model (1) can be written as(5)y=ρWy+xTα+BTzβ+ε,where y=y1,…,ynT, x=x1,…,xnT, ε=ε1,…,εnT, W=wil is an n×n specified constant spatial weight matrix, and BTz is an n×K1+⋯+Kp matrix with BTzi as its i-th row.The likelihood function corresponding to (5) is as follows:(6)Lα,β,σ2,ρ|y,x,z∝σ−nIn−ρWexp−12σ2y−ρWy−xTα−BTzβTy−ρWy−xTα−BTzβ=σ−nAρexp−12σ2Aρy−xTα−BTzβTAρy−xTα−BTzβ≐σ−nAρexp−12σ2Aρy−Bx,zθTAρy−Bx,zθ,where x=x1,…,xnT, z=z1,…,znT, K=∑j=1pKj, θ=αT,βTT is a K+q×1 vector of the regression coefficient, Bx,z=x,BTz is an n×K+q matrix, In is an identity matrix of order n, and Aρ=In−ρW.


## 3. Bayesian Estimation

We develop a Bayesian P-splines method with the Gibbs sampler to estimate the proposed model in this section. We begin with the specification of the prior distributions, then the derivations of the full conditional posteriors of all of the unknown parameters, and the narration of the detailed sampling scheme.


### 3.1. Priors

In order to develop a Bayesian P-splines method to estimate the model and eschew the use of improper to prevent improper joint posterior, we specify appropriate prior distributions for all the unknown parameters which include regression coefficientsθ=αT,βTT, variance σ2 , and spatial parameter ρ.Firstly, we set a hierarchical prior forα, which consists of a conjugate normal prior,(7)πα|σ2,τ0∝2πτ0σ2−q/2exp−αTα2τ0σ2,and an inverse-gamma prior,(8)πτ0∝τ0−rτα0/2−1exp−sτα022τ0,where rτα0 and sτα02 are prespecified hyperparameters.Secondly, to identify the unspecified smooth functiongj⋅ for j=1,…,p, the identifiability constraint Qjβj=0 should be imposed on βj. Under this constraint, we choose the random walk prior,(9)πβj|σ2,τj∝2πτjσ2−Kj−d/2exp−βjTMβjβj2τjσ2IQjβj=0,for βj, where I⋅ is the indicator function, d is the order of the random walk, Mβj is the penalty matrix that equals Dd−1×⋯×D0TDd−1×⋯×D0 for d-order random walk prior, Dl is a Kj−l−1×Kj−l matrix of the form(10)Dl=−110⋯00⋱⋱⋱⋮⋮⋱⋱⋱00⋯0−11,l=0,…,d−1.To enhance the flexibility and robustness of the method, we further put an inverse-gamma prior,(11)πτj∝τj−rτβj0/2−1exp−sτβj022τj,on τj, where rτβj0 and sτβj02 are prespecified hyperparameters. In addition, the following conjugate prior distribution of σ2 is assigned as(12)πσ2∝σ2−r0/2−1exp−s022σ2,where r0 and s02 are prespecified hyperparameters. Throughout this article, we set r0=s02=1 to obtain a Cauchy distribution of σ2 and use rτα0=rτβj0=1 and sτα02=sτβj02=0.005 to obtain a highly dispersed inverse-gamma prior for τj, j=0,1,…,p.Finally, the spatial autocorrelation coefficientρ is given a uniform prior ρ∼Uλmin−1,λmax−1, where λmin and λmax are the minimum and maximum eigenvalues of the standardized spatial weight matrix W, respectively. In other words,(13)πρ∝1.The joint prior of all unknown parameters is given by(14)πρ,α,β,σ2,τ=πρπσ2πτ0πα|σ2,τ0∏j=1pπβj|σ2,τjπτj,where τ=τ0,τ1,…,τp. Note that we have treated hyperparameter τ as a parameter vector for computational convenience.


### 3.2. The Full Conditional Posterior Distributions of the Parameters

In this subsection, we present the full conditional posterior distributions of all parameters that are used in the Gibbs sampler algorithm and describe the detailed sampling method.It follows from likelihood function (6) and priors (14) that, given the remaining unknown quantities, the conditional posterior distribution of spatial autocorrelation coefficient ρ is(15)pρ|y,x,z,α,β,σ2,τ∝Aρexp−12σ2Aρy−xTα−BTzβTAρy−xTα−BTzβ.Because it is not a standard density function, we cannot directly do simulations from (15). We prefer the Metropolis–Hastings algorithm [49, 50] to overcome the difficulty: generate a candidate ρ∗ from a truncated Cauchy distribution with location ρ and scale σρ on interval λmin−1,λmax−1, where σρ acts as a tuning parameter, and accept ρ∗ with probability(16)min1,pρ∗|y,x,z,α,β,σ2,τpρ|y,x,z,α,β,σ2,τ×Cρ,given the factor(17)Cρ=arctanλmax−1−ρ/σρ−arctanλmin−1−ρ/σρarctanλmax−1−ρ∗/σρ−arctanλmin−1−ρ∗/σρ.It is easy to see from likelihood function (6) and priors (14) that, given the spatial autocorrelation coefficient ρ and the hyperparameter τ, the joint posterior of all the parameters is proportional to(18)pθ,σ2|y,x,z,ρ,τ∝σ−nexp−12σ2Aρy−xTα−BTzβTAρy−xTα−BTzβ×σ−r0−2exp−s022σ2×σ−qexp−αTα2τ0σ2×∏j=1pσ−Kj+dexp−βjTMβjβj2τjσ2IQjβj=0∝exp−12σ2Aρy−Bx,zθTAρy−Bx,zθ×exp−12σ2θTdiagτ−1θIQβ=0×σ−n−r0−q−K+pd−2exp−s022σ2∝σ2−n+r0−pd/2−1exp−S2+s022σ2×2πσ2−K+q/2Ξ1/2exp−12σ2θ−θ^TΞθ−θ^IQβ=0,where Q=Q1,…,Qp, diagτ−1=diagτ0−1Iq,τ1−1Mβ1,…,τp−1Mβp, Ξ=diagτ−1+BTx,zBx,z, θ^=Ξ−1BTx,zAρy , and S2=yTAρTAρy−θ^TΞθ^.We can see from (18) that the method of composition [51] can be used to draw σ2 from the conditional inverse-gamma posterior,(19)pσ2|y,x,z,θ,ρ,τ∝σ2−n+r0−pd/2−1exp−S2+s022σ2,and to generate θ from the conditional normal posterior,(20)pθ|y,x,z,ρ,σ2,τ∝2πσ2−K+q/2Ξ1/2exp−12σ2θ−θ^TΞθ−θ^IQβ=0.In addition, to achieve identification, the constraintQjβj=0 should be imposed on βj. According to Panagiotelis and Smith [52], sampling βj from (20) is equivalent to sampling βjnew from Nβj∗,Σj∗, where Σj∗=τj−1Mβj+BjTBj−1, then βjnew is transformed to βj by(21)βj=βjnew−Σj∗QjTQjΣj∗QjT−1Qjβjnew,j=1,…,p.It is obvious that the hyperparameterτj for j=0,1,…,p is mutually independent in the posterior. The full conditional posterior of τj is an inverse-gamma distribution with density(22)pτ0|α,σ2∝τ0−q+rτα0/2−1exp−sτα02+αTα/σ22τ0,(23)pτj|β,σ2∝τj−Kj−d+rτβj0/2−1exp−sτβj02+βjTMβjβj/σ22τj,j=1,…,p,which can be simulated directly from (22) and (23), respectively.Algorithm 1: The pseudocode of the MCMC sampling scheme.
Input: Samples: yi,xi,zii=1,…,n.Initialization: Initialize the MCMC algorithm in iteration t=0 with Θ0, where all unknown quantities are drawn from the respective prior distributions.MCMC iterations: For t=1,2,3,…, given the current state Θt−1 successively, resample Θt from pΘ|y,x,z. The detailed Gibbs sampler cycles are outlined as follows:(i)
Drawρt from pρ|y,x,z,θt−1,σ2t−1,τt−1;(ii)
Drawσ2t from pσ2|y,x,z,θt−1,ρt−1,τt−1;(iii)
Drawθt from pθ|y,x,z,ρt−1,σ2t−1,τt−1, and adjust βt according to (21);(iv)
Drawτt from pτ|θt−1,σ2t−1, which is replaced by the following two steps:(v)
Drawτ0t from pτ0|αt−1,σ2t−1;(vi)
Drawτjt from pτj|βt−1,σ2t−1 for j=1,…,p.Output: An MCMC sample from the joint posterior distribution: Θtt=1,2,3,…..


### 3.3. Sampling

Let the parameter vectorΘ=ρ,σ2,θ,τ. The Bayesian estimate of Θ is obtained by observations drawn from the full conditional posterior of all parameters with running some Markov chain Monte Carlo (MCMC) tools such as the Gibbs sampler and the Metropolis–Hastings algorithm. The pseudocode of the MCMC algorithm in our method is implemented in Algorithm 1.


## 3.1. Priors

In order to develop a Bayesian P-splines method to estimate the model and eschew the use of improper to prevent improper joint posterior, we specify appropriate prior distributions for all the unknown parameters which include regression coefficientsθ=αT,βTT, variance σ2 , and spatial parameter ρ.Firstly, we set a hierarchical prior forα, which consists of a conjugate normal prior,(7)πα|σ2,τ0∝2πτ0σ2−q/2exp−αTα2τ0σ2,and an inverse-gamma prior,(8)πτ0∝τ0−rτα0/2−1exp−sτα022τ0,where rτα0 and sτα02 are prespecified hyperparameters.Secondly, to identify the unspecified smooth functiongj⋅ for j=1,…,p, the identifiability constraint Qjβj=0 should be imposed on βj. Under this constraint, we choose the random walk prior,(9)πβj|σ2,τj∝2πτjσ2−Kj−d/2exp−βjTMβjβj2τjσ2IQjβj=0,for βj, where I⋅ is the indicator function, d is the order of the random walk, Mβj is the penalty matrix that equals Dd−1×⋯×D0TDd−1×⋯×D0 for d-order random walk prior, Dl is a Kj−l−1×Kj−l matrix of the form(10)Dl=−110⋯00⋱⋱⋱⋮⋮⋱⋱⋱00⋯0−11,l=0,…,d−1.To enhance the flexibility and robustness of the method, we further put an inverse-gamma prior,(11)πτj∝τj−rτβj0/2−1exp−sτβj022τj,on τj, where rτβj0 and sτβj02 are prespecified hyperparameters. In addition, the following conjugate prior distribution of σ2 is assigned as(12)πσ2∝σ2−r0/2−1exp−s022σ2,where r0 and s02 are prespecified hyperparameters. Throughout this article, we set r0=s02=1 to obtain a Cauchy distribution of σ2 and use rτα0=rτβj0=1 and sτα02=sτβj02=0.005 to obtain a highly dispersed inverse-gamma prior for τj, j=0,1,…,p.Finally, the spatial autocorrelation coefficientρ is given a uniform prior ρ∼Uλmin−1,λmax−1, where λmin and λmax are the minimum and maximum eigenvalues of the standardized spatial weight matrix W, respectively. In other words,(13)πρ∝1.The joint prior of all unknown parameters is given by(14)πρ,α,β,σ2,τ=πρπσ2πτ0πα|σ2,τ0∏j=1pπβj|σ2,τjπτj,where τ=τ0,τ1,…,τp. Note that we have treated hyperparameter τ as a parameter vector for computational convenience.


## 3.2. The Full Conditional Posterior Distributions of the Parameters

In this subsection, we present the full conditional posterior distributions of all parameters that are used in the Gibbs sampler algorithm and describe the detailed sampling method.It follows from likelihood function (6) and priors (14) that, given the remaining unknown quantities, the conditional posterior distribution of spatial autocorrelation coefficient ρ is(15)pρ|y,x,z,α,β,σ2,τ∝Aρexp−12σ2Aρy−xTα−BTzβTAρy−xTα−BTzβ.Because it is not a standard density function, we cannot directly do simulations from (15). We prefer the Metropolis–Hastings algorithm [49, 50] to overcome the difficulty: generate a candidate ρ∗ from a truncated Cauchy distribution with location ρ and scale σρ on interval λmin−1,λmax−1, where σρ acts as a tuning parameter, and accept ρ∗ with probability(16)min1,pρ∗|y,x,z,α,β,σ2,τpρ|y,x,z,α,β,σ2,τ×Cρ,given the factor(17)Cρ=arctanλmax−1−ρ/σρ−arctanλmin−1−ρ/σρarctanλmax−1−ρ∗/σρ−arctanλmin−1−ρ∗/σρ.It is easy to see from likelihood function (6) and priors (14) that, given the spatial autocorrelation coefficient ρ and the hyperparameter τ, the joint posterior of all the parameters is proportional to(18)pθ,σ2|y,x,z,ρ,τ∝σ−nexp−12σ2Aρy−xTα−BTzβTAρy−xTα−BTzβ×σ−r0−2exp−s022σ2×σ−qexp−αTα2τ0σ2×∏j=1pσ−Kj+dexp−βjTMβjβj2τjσ2IQjβj=0∝exp−12σ2Aρy−Bx,zθTAρy−Bx,zθ×exp−12σ2θTdiagτ−1θIQβ=0×σ−n−r0−q−K+pd−2exp−s022σ2∝σ2−n+r0−pd/2−1exp−S2+s022σ2×2πσ2−K+q/2Ξ1/2exp−12σ2θ−θ^TΞθ−θ^IQβ=0,where Q=Q1,…,Qp, diagτ−1=diagτ0−1Iq,τ1−1Mβ1,…,τp−1Mβp, Ξ=diagτ−1+BTx,zBx,z, θ^=Ξ−1BTx,zAρy , and S2=yTAρTAρy−θ^TΞθ^.We can see from (18) that the method of composition [51] can be used to draw σ2 from the conditional inverse-gamma posterior,(19)pσ2|y,x,z,θ,ρ,τ∝σ2−n+r0−pd/2−1exp−S2+s022σ2,and to generate θ from the conditional normal posterior,(20)pθ|y,x,z,ρ,σ2,τ∝2πσ2−K+q/2Ξ1/2exp−12σ2θ−θ^TΞθ−θ^IQβ=0.In addition, to achieve identification, the constraintQjβj=0 should be imposed on βj. According to Panagiotelis and Smith [52], sampling βj from (20) is equivalent to sampling βjnew from Nβj∗,Σj∗, where Σj∗=τj−1Mβj+BjTBj−1, then βjnew is transformed to βj by(21)βj=βjnew−Σj∗QjTQjΣj∗QjT−1Qjβjnew,j=1,…,p.It is obvious that the hyperparameterτj for j=0,1,…,p is mutually independent in the posterior. The full conditional posterior of τj is an inverse-gamma distribution with density(22)pτ0|α,σ2∝τ0−q+rτα0/2−1exp−sτα02+αTα/σ22τ0,(23)pτj|β,σ2∝τj−Kj−d+rτβj0/2−1exp−sτβj02+βjTMβjβj/σ22τj,j=1,…,p,which can be simulated directly from (22) and (23), respectively.Algorithm 1: The pseudocode of the MCMC sampling scheme.
Input: Samples: yi,xi,zii=1,…,n.Initialization: Initialize the MCMC algorithm in iteration t=0 with Θ0, where all unknown quantities are drawn from the respective prior distributions.MCMC iterations: For t=1,2,3,…, given the current state Θt−1 successively, resample Θt from pΘ|y,x,z. The detailed Gibbs sampler cycles are outlined as follows:(i)
Drawρt from pρ|y,x,z,θt−1,σ2t−1,τt−1;(ii)
Drawσ2t from pσ2|y,x,z,θt−1,ρt−1,τt−1;(iii)
Drawθt from pθ|y,x,z,ρt−1,σ2t−1,τt−1, and adjust βt according to (21);(iv)
Drawτt from pτ|θt−1,σ2t−1, which is replaced by the following two steps:(v)
Drawτ0t from pτ0|αt−1,σ2t−1;(vi)
Drawτjt from pτj|βt−1,σ2t−1 for j=1,…,p.Output: An MCMC sample from the joint posterior distribution: Θtt=1,2,3,…..


## 3.3. Sampling

Let the parameter vectorΘ=ρ,σ2,θ,τ. The Bayesian estimate of Θ is obtained by observations drawn from the full conditional posterior of all parameters with running some Markov chain Monte Carlo (MCMC) tools such as the Gibbs sampler and the Metropolis–Hastings algorithm. The pseudocode of the MCMC algorithm in our method is implemented in Algorithm 1.


## 4. Numerical Illustration

A simulation study demonstrates that the proposed model and methodology perform satisfactorily in the finite sample. We apply the proposed model to analyze a real dataset.


### 4.1. Simulation

Consider the following model:(24)yi=ρ∑l=1nwilyl+xiTα+∑j=12gjzij+εi,i=1,…,n,where the covariate vector xi=xi1,xi2T follows a bivariate normal distribution with mean 0 and covariance matrix Σ=1−0.5−0.51, zi=zi1,zi2T is a bivariate vector, and both zi1 and zi2 are independent; they follow uniform distributions on −1,1 and 0,1, respectively. The link functions g1z=sinπz and g2z=4z1−z2−1, εi∼N0,σ2. The true values of parameters are assumed as α=1,−1T and σ2=0.25, respectively. For comparison, we consider three different cases of spatial parameters ρ1=0.25, ρ2=0.5, and ρ3=0.75, which represent from weak to strong spatial dependence of the response. Two kinds of spatial weight matrices W are used: one is Rook weight matrix [2] with n units, and the other is Case weight matrix [53] with r districts and m members in each district. The sample size of Rook weight matrix and Case weight matrix is n=100,169,400,625 and r,m=20,5,20,8,80,5,80,8, respectively.According to the above design, we conduct each simulation with 500 replications. Forj=1,…,p, we applied a natural cubic spline in which the knots with Kj=22 are placed at an equally spaced interval of the predictor variables, used second-order random walk penalties, and assigned hyperparameters r0,s02,rτα0,sτα02,rτβj0,sτβj02=1,1,1,0.005,1,0.005 in our computation. The second-order random walk penalties are used for the Bayesian P-splines to approximate the unknown smooth functions. The initial states of the Markov chain are chosen as follows: the unknown quantities are drawn from the respective prior distributions. By incrementally increasing or decreasing, the tuning parameter σρ is used such that the resultant acceptable rate for the parameter is around 25%.We can generate 6000 sampled values and delete the first 2000 values as a burn-in period for each of the replications by the proposed Gibbs sampler. In order to check the convergence of the proposed MCMC algorithm, we run five times with different starting values by the proposed Gibbs sampler for an arbitrarily selected replication. The sampled traces of parts of the unknown quantities for replication withr,m=80,5 and ρ=0.5 are displayed in Figure 1, where the unknown quantities include model parameters and fitted function on grid points. It is evident that the five parallel sequences mix reasonably well. The other cases have similar performance. Based on the five parallel sequences, we further calculate the “potential scale reduction factor” R^ [54] for each model parameter and fitted function on 10 selected grid points. All the values of R^ against the iteration numbers are displayed in Figure 2 (the case of the spatial parameter ρ=0.5). We observe that 2000 burn-in iterations have converged as all the values of R^ were less than 1.2. Based on 4000 sampled values, we then compute the means, the standard errors, and the 2.5th and 97.5th percentiles of the posterior distributions for all unknown quantities.Figure 1
Trace plots of five parallel sequences corresponding to different starting values for parts of the unknown quantities. (a)ρ, (b) α1, (c) α2, (d) σ2, (e) g1z15, (f) g1z35, (g) g1z55, (h) g1z75, (i) g2z15, (j) g2z35, (k) g2z55, and (l) g2z75.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 2
The “potential scale reduction factor”R^ for simulation results. (a) Rook weight matrix (n= 400). (b) Case weight matrix (r= 80, m= 5).
(a)(b)We evaluate the performance of the estimated functions by the mean absolute deviation errors (MADE):(25)MADEj=1101∑i=1101g^jzji−gjzjiand MADE=1p∑j=1pMADEj,at 101 fixed grid points zjii=1101 that are equally spaced chosen from interval aj,bj. To assess the accuracy of the estimates for the remaining parameters, we calculate the corresponding means across 500 replications for the posterior mean (Mean), bias (Bias), standard error (SE), and 2.5th and 97.5th percentiles of the parameters (the 95% posterior credible intervals, 95% CI). We also compute the standard derivations (SD) of the estimated posterior means to compare them with the means of the estimated posterior standard errors.Figure3(a) displays the boxplots of the MADE values with sample size n=100. Based on the Rook weight matrix, the medians are MADE1=0.1023, MADE2=0.0898, and MADE=0.0975. Based on the Case weight matrix, the medians are MADE1=0.1019, MADE2=0.0896, and MADE=0.0974. Figure 3(b) shows the boxplots of the MADE values with sample size n=400. Based on the Rook weight matrix, the medians are MADE1=0.0507, MADE2=0.0460, and MADE=0.0488. Based on the Case weight matrix, the medians are MADE1=0.0505, MADE2=0.0461, and MADE=0.0488. It shows that both the Rook and Case weight matrices can obtain reasonable estimations in fitting nonparametric functions.Figure 3
The boxplots (a) show the mean absolute deviation errors with sample sizen=100. The boxplots (b) show the mean absolute deviation errors with sample size n=400 (the three panels on the left are based on the Rook weight matrix and the three panels on the right are based on the Case weight matrix). (a) MADE (ρ = 0.5, n= 100), (b) MADE (ρ = 0.5, n= 400).
(a)(b)Table1 summarizes the estimation results. It is observed that all of the means of the estimates are very close to the respective true values, and the average values of the standard errors are close to the corresponding standard derivations, which indicate that the parameter estimates and the standard errors are accurate. This indicates that the proposed estimators of the parametric model perform better with increasing sample size. By comparing the estimates for ρ under the same sample sizes, we find the estimates of ρ with the Case weight matrix are much better than those with the Rook weight matrix. The possible main reason lies in the fact that the Rook weight matrix is more complicated than the Case weight matrix. Finally, the bigger the sample size under the same spatial complexity, the more accurate the estimates are. We have also repeated the aforementioned experiences with different starting values, and the results are similar. This indicates that the proposed Gibbs sampler works well.Table 1
Simulation results of parameter estimation.
ParameternRook weight matrixr,mCase weight matrixMeanBiasSESD95% CIMeanBiasSESD95% CIρ=0.2501000.246−0.0040.0480.0470.152,0.340(20, 5)0.245−0.0050.0430.0440.160,0.329α1=1.0000.992−0.0080.0620.0610.870,1.1140.992−0.0080.0620.0600.870,1.114α2=−1.000−0.9990.0010.0620.063−1.122,−0.877−1.0000.0000.0620.063−1.122,−0.877σ2=0.2500.2510.0010.0380.0350.187,0.3350.2510.0010.0380.0350.187,0.335ρ=0.5000.494−0.0060.0410.0410.412,0.5740.495−0.0050.0310.0310.434,0.554α1=1.0000.993−0.0070.0630.0610.870,1.1150.993−0.0070.0630.0600.870,1.115α2=−1.000−1.0000.0000.0620.063−1.123,−0.877−1.0000.0000.0620.063−1.123,−0.878σ2=0.2500.2510.0010.0380.0350.187,0.3360.2510.0010.0380.0350.187,0.336ρ=0.7500.744−0.0060.0280.0290.687,0.7990.747−0.0030.0160.0170.715,0.778α1=1.0000.993−0.0070.0630.0620.870,1.1170.994−0.0060.0630.0600.870,1.117α2=−1.000−1.001−0.0010.0630.063−1.124,−0.878−1.001−0.0010.0630.064−1.125,−0.878σ2=0.2500.2510.0010.0380.0350.187,0.3370.2510.0010.0390.0350.187,0.337ρ=0.2501690.2500.0000.0360.0360.180,0.320(20, 8)0.240−0.0100.0430.0410.157,0.323α1=1.0000.994−0.0060.0460.0460.903,1.0851.0010.0010.0480.0490.907,1.095α2=−1.000−0.9980.0020.0460.046−1.089,−0.907−0.9950.0050.0480.049−1.088,−0.901σ2=0.2500.249−0.0010.0280.0270.200,0.3100.2500.0000.0290.0280.199,0.313ρ=0.5000.4980.0020.0300.0310.438,0.5580.4930.0070.0290.0280.435,0.550α1=1.0000.994−0.0060.0460.0460.903,1.0861.0010.0010.0480.0490.908,1.095α2=−1.000−0.9980.0020.0460.046−1.089,−0.907−0.9950.0050.0480.049−1.089,−0.901σ2=0.2500.249−0.0010.0280.0280.200,0.3100.2500.0000.0290.0280.199,0.313ρ=0.7500.748−0.0020.0210.0210.707,0.7880.746−0.0040.0150.0150.716,0.775α1=1.0000.995−0.0050.0470.0460.904,1.0871.0020.0020.0480.0490.908,1.096α2=−1.000−0.9990.0010.0470.046−1.090,−0.907−0.9960.0040.0480.049−1.090,−0.902σ2=0.2500.249−0.0010.0280.0280.200,0.3110.2500.0000.0290.0290.199,0.314ρ=0.2504000.249−0.0010.0230.0230.204,0.29480,50.248−0.0020.0200.0200.209,0.286α1=1.0001.0010.0010.0290.0280.943,1.0591.0010.0010.0290.0280.943,1.059α2=−1.000−0.9970.0030.0290.030−1.055,−0.940−0.9970.0030.0290.030−1.055,−0.940σ2=0.2500.2500.0000.0180.0180.217,0.2880.2500.0000.0180.180.217,0.288ρ=0.5000.499−0.0010.0200.0200.460,0.5370.498−0.0020.0140.0140.470,0.526α1=1.0001.0010.0010.0300.0300.943,1.0591.0010.0010.0300.0300.943,1.059α2=−1.000−0.9970.0030.0300.030−1.055,−0.940−0.9970.0030.0300.030−1.055,−0.940σ2=0.2500.2500.0000.0180.0180.217,0.2880.2500.0000.0180.00180.217,0.288ρ=0.7500.749−0.0010.0130.0140.722,0.7750.749−0.0010.0070.0070.734,0.763α1=1.0001.0010.0010.0300.0290.943,1.0601.0020.0020.0300.0290.943,1.060α2=−1.000−0.9980.0020.0300.029−1.056,−0.939−0.9980.0020.0300.030−1.056,−0.940σ2=0.2500.2500.0000.0180.0180.217,0.2880.2500.0000.0180.0180.217,0.289ρ=0.2506250.247−0.0030.0180.0200.211,0.28380,80.2470.0030.0200.0190.208,0.287α1=1.0001.0000.0000.0230.0220.954,1.0460.998−0.0020.0230.0240.953,1.044α2=−1.000−0.9990.0010.0230.023−1.045,−0.953−1.002−0.0020.0230.023−1.047,−0.957σ2=0.2500.2510.0010.0140.0140.224,0.2800.2500.0000.0140.0140.224,0.279ρ=0.5000.497−0.0030.0160.0170.467,0.5280.498−0.0020.0140.0130.471,0.525α1=1.0001.0010.0010.0240.0230.955,1.0470.999−0.0010.0230.0240.953,1.044α2=−1.000−0.9990.0010.0240.023−1.045,−0.953−1.002−0.0020.0230.023−1.047,−0.957σ2=0.2500.2510.0010.0140.0140.224,0.2800.2500.0000.0140.0140.224,0.279ρ=0.7500.748−0.0020.0110.0110.727,0.7690.749−0.0010.0070.0070.735,0.763α1=1.0001.0010.0010.0240.0230.955,1.0470.999−0.0010.0230.0240.953,1.044α2=−1.000−0.9990.0010.0240.023−1.046,−0.953−1.002−0.0020.0230.023−1.048,−0.957σ2=0.2500.2510.0010.0140.0140.224,0.2810.2500.0000.0140.0140.224,0.280The median and SD of MADE values of the link functions are reported in Table2. It shows that the median and SD of MADE values of the link functions decrease along with increasing the sample size, which indicates that the estimates of the unknown functions are convergent. The spatial correlation coefficients also have little effect on the nonparametric estimates.Table 2
The median and SD of MADE values for simulation results.
Para.Rookg1zg2zCaseg1zg2zρnMedianSDMedianSDr,mMedianSDMedianSD0.2501000.1020.0340.0900.03120,50.1020.0340.0900.0300.5000.1020.0340.0900.0310.1020.0340.0900.0310.7500.1020.0340.0890.0310.1020.0340.0890.0310.2501690.0760.0260.0700.02320,80.0790.0270.0730.0230.5000.0760.0260.0700.0230.0790.0270.0730.0230.7500.0760.0260.0700.0230.0790.0270.0730.0230.2504000.0510.0160.0460.01480,50.0510.0160.0460.0140.5000.0510.0160.0460.0140.0510.0160.0460.0140.7500.0510.0160.0460.0140.0510.0160.0460.0140.2506250.0420.0130.0360.01280,80.0400.0140.0370.0110.5000.0420.0130.0370.0130.0400.0140.0370.0110.7500.0420.0130.0360.0130.0400.0140.0370.011The estimated functions, together with the 95% pointwise posterior credible intervals, are depicted in Figure4 from a typical sample under the spatial dependence ρ=0.5 when the sample size is n=100 and n=400, respectively. The typical sample is selected in such a way that its MADE value is equal to the median in the 500 replications. This indicates that the proposed estimators of nonparametric functions are performed better with increasing sample size. It took about 30 seconds for n=100 and 60 seconds for n=400 to run each replication using a PC with Intel (R) Core (TM) i7-8750H 2.20 GHz, respectively. Computer code is available upon request from the authors.Figure 4
The estimated functions (dotted lines) and their 95% pointwise posterior credible intervals (dot-dashed lines) for a typical sample. The solid lines represent the true functions (the two panels on the left are based on the Rook weight matrix and the two panels on the right are based on the Case weight matrix).
(a)(b)(c)(d)(e)(f)(g)(h)


### 4.2. A Real Data Example

As an illustration of our proposed model and method with a real example, we analyze housing prices data in Chinese cities. The data are obtained from the China City Statistical Yearbook 2013 [56] and the China Statistical Yearbook for Regional Economy 2013 [57]. They cover 286 cities at/above the prefecture level (except for the cities in Taiwan, Hong Kong, and Macau). In our real data analysis, we consider that the dependent variable y is each city’s average selling price of residential houses (denoted by HP). There are three covariates: (1) population density (PD); (2) per capita disposable income of urban households (INC); and (3) loan-to-GDP ratio (MON), which indicates the degree of monetary policy whether ease or tightness.This motivates us to consider the following partially linear additive spatial autoregressive model:(26)yi=ρ∑l=1nwilyl+αxi+∑j=12gjzij+εi,i=1,…,n,where the response variable yi=HPi, xi1=PDi, zi1=INCi, zi2=MONi. With regard to the choice of the weight matrix, according to the practice in Sun et al. [24], we use the Euclidean distance in terms of any two houses to calculate the spatial weight matrix W. We use longitude and latitude to represent the location, which is denoted as si=Loni,Lati. The spatial weight wil is(27)wil=exp−si−sl/∑k≠iexp−si−sk.We set as(28)PD=PD1/3−minPD1/3maxPD1/3−minPD1/3,where the transformation is used for transferring the asymmetric distribution of PD to nearly uniform distribution on 0,1. The other covariates variables of INC and MON are transformed so that the marginal distribution is approximately N0,1. These do not alter the model but facilitate implementation. For this dataset, we adopted natural cubic P-splines, used second-order random walk penalties, and assigned λ,r0,s02,rτα0,sτα02,rτβj0,sτβj02=2,1,1,1,0.005,1,0.005. The tuning parameter σρ is used such that the acceptable rate for updating ρ is around 25%.We run five times with different initial states by the proposed Gibbs sampler and generate 10000 sampled values following a burn-in of 20000 iterations in each replication. Figure5 plots the traces of parts of the unknown quantities. It can be seen that the five parallel sequences mix very well. Based on the five parallel sequences, the “potential scale reduction factor” R^ is calculated and plotted in Figure 6, from which we observe that the Markov chain has converged within 20000 burn-in iterations.Figure 5
Trace plots of five parallel sequences corresponding to different starting values for parts of the unknown quantities. (a)α, (b) g1z15, (c) g1z35, (d) g1z55, (e) g1z75, (f) g2z25, (g) g2z45, (h) g2z65, and (i) g2z85.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 6
The “potential scale reduction factor”R^ for the housing prices of Chinese cities.The estimated parameters together with their standard errors and 95% posterior credible intervals are reported in Table3. The estimate of the spatial parameter is ρ^=0.615 with the 0.034 standard deviation which indicates there exists a strong and positive spatial relationship among the response. We observe that the covariate PD has significant effects on the response. The regression coefficient of PD is α^=0.258>0, in which the population density has a positive effect on the housing price.Table 3
The parameter estimation for the housing prices of Chinese cities.
ParameterMeanSE95% CIα0.2580.0730.112,0.401ρ0.6150.0340.550,0.680σ20.01680.00150.0141,0.0199The estimated functions, together with the 95% pointwise posterior credible intervals, are depicted in Figure7, which look like two nonlinear functions with an upward trend. The curves show that g1z1 and g2z2 have a local maximum of 0.764 at around z1=3.335 and 0.272 at around z2=2.592, respectively. It provides evidence that the effects of the covariates INC and MON on the response are significant and nonlinear.Figure 7
The estimated functions (dot line) and their 95% pointwise posterior credible intervals (dot-dashed lines) for the housing prices of Chinese cities in model (26). (a) g1z1 and (b) g2z2.
(a)(b)


## 4.1. Simulation

Consider the following model:(24)yi=ρ∑l=1nwilyl+xiTα+∑j=12gjzij+εi,i=1,…,n,where the covariate vector xi=xi1,xi2T follows a bivariate normal distribution with mean 0 and covariance matrix Σ=1−0.5−0.51, zi=zi1,zi2T is a bivariate vector, and both zi1 and zi2 are independent; they follow uniform distributions on −1,1 and 0,1, respectively. The link functions g1z=sinπz and g2z=4z1−z2−1, εi∼N0,σ2. The true values of parameters are assumed as α=1,−1T and σ2=0.25, respectively. For comparison, we consider three different cases of spatial parameters ρ1=0.25, ρ2=0.5, and ρ3=0.75, which represent from weak to strong spatial dependence of the response. Two kinds of spatial weight matrices W are used: one is Rook weight matrix [2] with n units, and the other is Case weight matrix [53] with r districts and m members in each district. The sample size of Rook weight matrix and Case weight matrix is n=100,169,400,625 and r,m=20,5,20,8,80,5,80,8, respectively.According to the above design, we conduct each simulation with 500 replications. Forj=1,…,p, we applied a natural cubic spline in which the knots with Kj=22 are placed at an equally spaced interval of the predictor variables, used second-order random walk penalties, and assigned hyperparameters r0,s02,rτα0,sτα02,rτβj0,sτβj02=1,1,1,0.005,1,0.005 in our computation. The second-order random walk penalties are used for the Bayesian P-splines to approximate the unknown smooth functions. The initial states of the Markov chain are chosen as follows: the unknown quantities are drawn from the respective prior distributions. By incrementally increasing or decreasing, the tuning parameter σρ is used such that the resultant acceptable rate for the parameter is around 25%.We can generate 6000 sampled values and delete the first 2000 values as a burn-in period for each of the replications by the proposed Gibbs sampler. In order to check the convergence of the proposed MCMC algorithm, we run five times with different starting values by the proposed Gibbs sampler for an arbitrarily selected replication. The sampled traces of parts of the unknown quantities for replication withr,m=80,5 and ρ=0.5 are displayed in Figure 1, where the unknown quantities include model parameters and fitted function on grid points. It is evident that the five parallel sequences mix reasonably well. The other cases have similar performance. Based on the five parallel sequences, we further calculate the “potential scale reduction factor” R^ [54] for each model parameter and fitted function on 10 selected grid points. All the values of R^ against the iteration numbers are displayed in Figure 2 (the case of the spatial parameter ρ=0.5). We observe that 2000 burn-in iterations have converged as all the values of R^ were less than 1.2. Based on 4000 sampled values, we then compute the means, the standard errors, and the 2.5th and 97.5th percentiles of the posterior distributions for all unknown quantities.Figure 1
Trace plots of five parallel sequences corresponding to different starting values for parts of the unknown quantities. (a)ρ, (b) α1, (c) α2, (d) σ2, (e) g1z15, (f) g1z35, (g) g1z55, (h) g1z75, (i) g2z15, (j) g2z35, (k) g2z55, and (l) g2z75.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Figure 2
The “potential scale reduction factor”R^ for simulation results. (a) Rook weight matrix (n= 400). (b) Case weight matrix (r= 80, m= 5).
(a)(b)We evaluate the performance of the estimated functions by the mean absolute deviation errors (MADE):(25)MADEj=1101∑i=1101g^jzji−gjzjiand MADE=1p∑j=1pMADEj,at 101 fixed grid points zjii=1101 that are equally spaced chosen from interval aj,bj. To assess the accuracy of the estimates for the remaining parameters, we calculate the corresponding means across 500 replications for the posterior mean (Mean), bias (Bias), standard error (SE), and 2.5th and 97.5th percentiles of the parameters (the 95% posterior credible intervals, 95% CI). We also compute the standard derivations (SD) of the estimated posterior means to compare them with the means of the estimated posterior standard errors.Figure3(a) displays the boxplots of the MADE values with sample size n=100. Based on the Rook weight matrix, the medians are MADE1=0.1023, MADE2=0.0898, and MADE=0.0975. Based on the Case weight matrix, the medians are MADE1=0.1019, MADE2=0.0896, and MADE=0.0974. Figure 3(b) shows the boxplots of the MADE values with sample size n=400. Based on the Rook weight matrix, the medians are MADE1=0.0507, MADE2=0.0460, and MADE=0.0488. Based on the Case weight matrix, the medians are MADE1=0.0505, MADE2=0.0461, and MADE=0.0488. It shows that both the Rook and Case weight matrices can obtain reasonable estimations in fitting nonparametric functions.Figure 3
The boxplots (a) show the mean absolute deviation errors with sample sizen=100. The boxplots (b) show the mean absolute deviation errors with sample size n=400 (the three panels on the left are based on the Rook weight matrix and the three panels on the right are based on the Case weight matrix). (a) MADE (ρ = 0.5, n= 100), (b) MADE (ρ = 0.5, n= 400).
(a)(b)Table1 summarizes the estimation results. It is observed that all of the means of the estimates are very close to the respective true values, and the average values of the standard errors are close to the corresponding standard derivations, which indicate that the parameter estimates and the standard errors are accurate. This indicates that the proposed estimators of the parametric model perform better with increasing sample size. By comparing the estimates for ρ under the same sample sizes, we find the estimates of ρ with the Case weight matrix are much better than those with the Rook weight matrix. The possible main reason lies in the fact that the Rook weight matrix is more complicated than the Case weight matrix. Finally, the bigger the sample size under the same spatial complexity, the more accurate the estimates are. We have also repeated the aforementioned experiences with different starting values, and the results are similar. This indicates that the proposed Gibbs sampler works well.Table 1
Simulation results of parameter estimation.
ParameternRook weight matrixr,mCase weight matrixMeanBiasSESD95% CIMeanBiasSESD95% CIρ=0.2501000.246−0.0040.0480.0470.152,0.340(20, 5)0.245−0.0050.0430.0440.160,0.329α1=1.0000.992−0.0080.0620.0610.870,1.1140.992−0.0080.0620.0600.870,1.114α2=−1.000−0.9990.0010.0620.063−1.122,−0.877−1.0000.0000.0620.063−1.122,−0.877σ2=0.2500.2510.0010.0380.0350.187,0.3350.2510.0010.0380.0350.187,0.335ρ=0.5000.494−0.0060.0410.0410.412,0.5740.495−0.0050.0310.0310.434,0.554α1=1.0000.993−0.0070.0630.0610.870,1.1150.993−0.0070.0630.0600.870,1.115α2=−1.000−1.0000.0000.0620.063−1.123,−0.877−1.0000.0000.0620.063−1.123,−0.878σ2=0.2500.2510.0010.0380.0350.187,0.3360.2510.0010.0380.0350.187,0.336ρ=0.7500.744−0.0060.0280.0290.687,0.7990.747−0.0030.0160.0170.715,0.778α1=1.0000.993−0.0070.0630.0620.870,1.1170.994−0.0060.0630.0600.870,1.117α2=−1.000−1.001−0.0010.0630.063−1.124,−0.878−1.001−0.0010.0630.064−1.125,−0.878σ2=0.2500.2510.0010.0380.0350.187,0.3370.2510.0010.0390.0350.187,0.337ρ=0.2501690.2500.0000.0360.0360.180,0.320(20, 8)0.240−0.0100.0430.0410.157,0.323α1=1.0000.994−0.0060.0460.0460.903,1.0851.0010.0010.0480.0490.907,1.095α2=−1.000−0.9980.0020.0460.046−1.089,−0.907−0.9950.0050.0480.049−1.088,−0.901σ2=0.2500.249−0.0010.0280.0270.200,0.3100.2500.0000.0290.0280.199,0.313ρ=0.5000.4980.0020.0300.0310.438,0.5580.4930.0070.0290.0280.435,0.550α1=1.0000.994−0.0060.0460.0460.903,1.0861.0010.0010.0480.0490.908,1.095α2=−1.000−0.9980.0020.0460.046−1.089,−0.907−0.9950.0050.0480.049−1.089,−0.901σ2=0.2500.249−0.0010.0280.0280.200,0.3100.2500.0000.0290.0280.199,0.313ρ=0.7500.748−0.0020.0210.0210.707,0.7880.746−0.0040.0150.0150.716,0.775α1=1.0000.995−0.0050.0470.0460.904,1.0871.0020.0020.0480.0490.908,1.096α2=−1.000−0.9990.0010.0470.046−1.090,−0.907−0.9960.0040.0480.049−1.090,−0.902σ2=0.2500.249−0.0010.0280.0280.200,0.3110.2500.0000.0290.0290.199,0.314ρ=0.2504000.249−0.0010.0230.0230.204,0.29480,50.248−0.0020.0200.0200.209,0.286α1=1.0001.0010.0010.0290.0280.943,1.0591.0010.0010.0290.0280.943,1.059α2=−1.000−0.9970.0030.0290.030−1.055,−0.940−0.9970.0030.0290.030−1.055,−0.940σ2=0.2500.2500.0000.0180.0180.217,0.2880.2500.0000.0180.180.217,0.288ρ=0.5000.499−0.0010.0200.0200.460,0.5370.498−0.0020.0140.0140.470,0.526α1=1.0001.0010.0010.0300.0300.943,1.0591.0010.0010.0300.0300.943,1.059α2=−1.000−0.9970.0030.0300.030−1.055,−0.940−0.9970.0030.0300.030−1.055,−0.940σ2=0.2500.2500.0000.0180.0180.217,0.2880.2500.0000.0180.00180.217,0.288ρ=0.7500.749−0.0010.0130.0140.722,0.7750.749−0.0010.0070.0070.734,0.763α1=1.0001.0010.0010.0300.0290.943,1.0601.0020.0020.0300.0290.943,1.060α2=−1.000−0.9980.0020.0300.029−1.056,−0.939−0.9980.0020.0300.030−1.056,−0.940σ2=0.2500.2500.0000.0180.0180.217,0.2880.2500.0000.0180.0180.217,0.289ρ=0.2506250.247−0.0030.0180.0200.211,0.28380,80.2470.0030.0200.0190.208,0.287α1=1.0001.0000.0000.0230.0220.954,1.0460.998−0.0020.0230.0240.953,1.044α2=−1.000−0.9990.0010.0230.023−1.045,−0.953−1.002−0.0020.0230.023−1.047,−0.957σ2=0.2500.2510.0010.0140.0140.224,0.2800.2500.0000.0140.0140.224,0.279ρ=0.5000.497−0.0030.0160.0170.467,0.5280.498−0.0020.0140.0130.471,0.525α1=1.0001.0010.0010.0240.0230.955,1.0470.999−0.0010.0230.0240.953,1.044α2=−1.000−0.9990.0010.0240.023−1.045,−0.953−1.002−0.0020.0230.023−1.047,−0.957σ2=0.2500.2510.0010.0140.0140.224,0.2800.2500.0000.0140.0140.224,0.279ρ=0.7500.748−0.0020.0110.0110.727,0.7690.749−0.0010.0070.0070.735,0.763α1=1.0001.0010.0010.0240.0230.955,1.0470.999−0.0010.0230.0240.953,1.044α2=−1.000−0.9990.0010.0240.023−1.046,−0.953−1.002−0.0020.0230.023−1.048,−0.957σ2=0.2500.2510.0010.0140.0140.224,0.2810.2500.0000.0140.0140.224,0.280The median and SD of MADE values of the link functions are reported in Table2. It shows that the median and SD of MADE values of the link functions decrease along with increasing the sample size, which indicates that the estimates of the unknown functions are convergent. The spatial correlation coefficients also have little effect on the nonparametric estimates.Table 2
The median and SD of MADE values for simulation results.
Para.Rookg1zg2zCaseg1zg2zρnMedianSDMedianSDr,mMedianSDMedianSD0.2501000.1020.0340.0900.03120,50.1020.0340.0900.0300.5000.1020.0340.0900.0310.1020.0340.0900.0310.7500.1020.0340.0890.0310.1020.0340.0890.0310.2501690.0760.0260.0700.02320,80.0790.0270.0730.0230.5000.0760.0260.0700.0230.0790.0270.0730.0230.7500.0760.0260.0700.0230.0790.0270.0730.0230.2504000.0510.0160.0460.01480,50.0510.0160.0460.0140.5000.0510.0160.0460.0140.0510.0160.0460.0140.7500.0510.0160.0460.0140.0510.0160.0460.0140.2506250.0420.0130.0360.01280,80.0400.0140.0370.0110.5000.0420.0130.0370.0130.0400.0140.0370.0110.7500.0420.0130.0360.0130.0400.0140.0370.011The estimated functions, together with the 95% pointwise posterior credible intervals, are depicted in Figure4 from a typical sample under the spatial dependence ρ=0.5 when the sample size is n=100 and n=400, respectively. The typical sample is selected in such a way that its MADE value is equal to the median in the 500 replications. This indicates that the proposed estimators of nonparametric functions are performed better with increasing sample size. It took about 30 seconds for n=100 and 60 seconds for n=400 to run each replication using a PC with Intel (R) Core (TM) i7-8750H 2.20 GHz, respectively. Computer code is available upon request from the authors.Figure 4
The estimated functions (dotted lines) and their 95% pointwise posterior credible intervals (dot-dashed lines) for a typical sample. The solid lines represent the true functions (the two panels on the left are based on the Rook weight matrix and the two panels on the right are based on the Case weight matrix).
(a)(b)(c)(d)(e)(f)(g)(h)


## 4.2. A Real Data Example

As an illustration of our proposed model and method with a real example, we analyze housing prices data in Chinese cities. The data are obtained from the China City Statistical Yearbook 2013 [56] and the China Statistical Yearbook for Regional Economy 2013 [57]. They cover 286 cities at/above the prefecture level (except for the cities in Taiwan, Hong Kong, and Macau). In our real data analysis, we consider that the dependent variable y is each city’s average selling price of residential houses (denoted by HP). There are three covariates: (1) population density (PD); (2) per capita disposable income of urban households (INC); and (3) loan-to-GDP ratio (MON), which indicates the degree of monetary policy whether ease or tightness.This motivates us to consider the following partially linear additive spatial autoregressive model:(26)yi=ρ∑l=1nwilyl+αxi+∑j=12gjzij+εi,i=1,…,n,where the response variable yi=HPi, xi1=PDi, zi1=INCi, zi2=MONi. With regard to the choice of the weight matrix, according to the practice in Sun et al. [24], we use the Euclidean distance in terms of any two houses to calculate the spatial weight matrix W. We use longitude and latitude to represent the location, which is denoted as si=Loni,Lati. The spatial weight wil is(27)wil=exp−si−sl/∑k≠iexp−si−sk.We set as(28)PD=PD1/3−minPD1/3maxPD1/3−minPD1/3,where the transformation is used for transferring the asymmetric distribution of PD to nearly uniform distribution on 0,1. The other covariates variables of INC and MON are transformed so that the marginal distribution is approximately N0,1. These do not alter the model but facilitate implementation. For this dataset, we adopted natural cubic P-splines, used second-order random walk penalties, and assigned λ,r0,s02,rτα0,sτα02,rτβj0,sτβj02=2,1,1,1,0.005,1,0.005. The tuning parameter σρ is used such that the acceptable rate for updating ρ is around 25%.We run five times with different initial states by the proposed Gibbs sampler and generate 10000 sampled values following a burn-in of 20000 iterations in each replication. Figure5 plots the traces of parts of the unknown quantities. It can be seen that the five parallel sequences mix very well. Based on the five parallel sequences, the “potential scale reduction factor” R^ is calculated and plotted in Figure 6, from which we observe that the Markov chain has converged within 20000 burn-in iterations.Figure 5
Trace plots of five parallel sequences corresponding to different starting values for parts of the unknown quantities. (a)α, (b) g1z15, (c) g1z35, (d) g1z55, (e) g1z75, (f) g2z25, (g) g2z45, (h) g2z65, and (i) g2z85.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 6
The “potential scale reduction factor”R^ for the housing prices of Chinese cities.The estimated parameters together with their standard errors and 95% posterior credible intervals are reported in Table3. The estimate of the spatial parameter is ρ^=0.615 with the 0.034 standard deviation which indicates there exists a strong and positive spatial relationship among the response. We observe that the covariate PD has significant effects on the response. The regression coefficient of PD is α^=0.258>0, in which the population density has a positive effect on the housing price.Table 3
The parameter estimation for the housing prices of Chinese cities.
ParameterMeanSE95% CIα0.2580.0730.112,0.401ρ0.6150.0340.550,0.680σ20.01680.00150.0141,0.0199The estimated functions, together with the 95% pointwise posterior credible intervals, are depicted in Figure7, which look like two nonlinear functions with an upward trend. The curves show that g1z1 and g2z2 have a local maximum of 0.764 at around z1=3.335 and 0.272 at around z2=2.592, respectively. It provides evidence that the effects of the covariates INC and MON on the response are significant and nonlinear.Figure 7
The estimated functions (dot line) and their 95% pointwise posterior credible intervals (dot-dashed lines) for the housing prices of Chinese cities in model (26). (a) g1z1 and (b) g2z2.
(a)(b)


## 5. Summary

Spatial data are often encountered in economics, geography, and environmental science and can be analyzed by the spatial autoregressive models. To reduce the high risk of misspecification of the traditional spatial autoregressive models and avoid some serious drawbacks of fully nonparametric models, we have developed a PLASARM for spatial data, which combines the partially linear additive model and spatial autoregressive model. In this paper, we have considered a Bayesian MCMC approach with P-splines to analyze the proposed model and designed a Gibbs sampler to explore the joint posterior distributions based on the Bayesian P-splines technique. A simulation study demonstrates that the proposed model and methodology perform satisfactorily in the finite sample. The results have shown that they are reliable for recovering the true parameters and nonparametric components. We have further applied the proposed model and methodology to analyze a real dataset.Various extensions can be considered in the future. We can pursue the case where the covariates are high dimension in the model. The proposed models can be generalized in several ways. Although we focus on a PLASARM to assess the effects of the covariates on the response, the other semiparametric models such as single-index spatial autoregressive models and varying-coefficient spatial autoregressive models can also be considered. In addition, it is straightforward to generalize the model and methodology to handle more complicated data types, for example, mixed discrete, continuous, and ordinal data [55]. We leave these topics for future research.



---

*Source: 1777469-2021-07-15.xml*