# Bayesian Adaptive Lasso for the Partial Functional Linear Spatial Autoregressive Model

Authors: Dengke Xu; Ruiqin Tian; Ying Lu
Journal: Journal of Mathematics
 (2022)

---

## Abstract

This study introduces a partial functional linear spatial autoregressive model which can explore the relationship between a scalar spatially dependent response variable and predictive variables containing both multiple scalar covariates and a functional covariate. With approximating to the functional coefficient by Karhunen–Loève representation, we propose a Bayesian adaptive Lasso method to simultaneously estimate unknown parameters and select important covariates in the model, which can be performed by combining the Gibbs sampler and the Metropolis–Hastings algorithm. Some simulation studies are conducted and the results show that the proposed Bayesian method behaves well.


---

## Body

## 1. Introduction

With the rapid development of electronic technology and the emergence of advanced measurement tools, a large amount of data can be collected and stored cheaply. In particular, scientists usually confront the data providing information about curves, surfaces, or anything else varying with continuous variables, such as economic index, nuclear magnetic resonance data, weather data, and stock data. This type of data structure is called functional data and has been received an increasing level of attention in various applied fields, including financial engineering, environmental science, medicine, brain imaging, and public health. Therefore, many useful functional regression models have been proposed, and there has been a large amount of work in functional data analysis; see [1–5], among others. Especially, it often occurs in data modeling that a response is related to predictive variables including scalar covariates and a functional covariate, which refers to the partial functional linear model (PFLM). Recently, a lot of research work has been contributed to PFLM. For example, Shin [6] established the asymptotic normality of the estimated regression coefficients and the convergence rate of the estimated functional coefficient of the proposed partial functional linear regression model. Based on the results of Shin [6], Lu et al. [7] introduced a functional partially quantile regression model and gave the large sample properties of the proposed estimators. Yu et al. [8] investigated the hypothesis test of the parametric component in partial functional linear regression. Based on the incomplete data, Zhou and Peng [4] discussed a functional partial linear model when some responses are missing at random. Yu et al. [9] studied a single-index partially functional linear regression model and obtained the convergence rates and asymptotic normality of the resulting estimators. To sum up, most of the above references are concerned that all individuals in functional regression models are mutually independent and the response variables do not have a spatial dependence structure. To our knowledge, relatively few studies have been carried out when response variables are spatially dependent variables in functional data analysis.It is well known that a large number of models and methods are used to deal with spatially dependent variables. One useful approach in dealing with spatial dependence is the spatial autoregressive model, which has received great interest of both statisticians and econometricians. For example, Liu et al. [10] investigated variable selection in the spatial autoregressive model with independent and identical distributed errors based on a penalized quasi-maximum likelihood method. Xie et al. [11] considered conducting variable selection in spatial autoregressive models with a diverging number of parameters. Based on the idea of B-spline approximation and instrumental variables, Du et al. [12] proposed an estimate method for the partially linear additive spatial autoregressive models and obtained the convergence rate and asymptotic normality of the resulting estimators. Su and Jin [13] developed statistical inference for partially linear spatial autoregressive models by using the kernel estimate method and quasi-maximum likelihood method. However, the abovementioned papers focused mainly on that the predictor variables are real-valued variables.On the contrary, due to the recent dramatic evolution in advanced computational technologies, many subject areas based on numerical approximations have developed rapidly, such as [14, 15]. One of the biggest beneficiaries of these fields is Bayesian statistics. Specifically, Bayesian models allow one to incorporate appropriate level of uncertainty; thus, Bayesian works for various statistical models have also received substantial attentions in recent years. For example, based on the variance modeling technique and spline approximation, Xu and Zhang [16] discussed Bayesian inference for semiparametric joint mean and variance models. Combining the Gibbs sampler with the Metropolis–Hastings algorithm, Tang et al. [17] investigated Bayesian estimation and Bayesian local influence analysis of the transformation linear mixed models that the random effects follow an unknown distribution. Pfarrhofer and Piribauer [18] proposed two shrinkage priors to make Bayesian variable selection for high-dimensional spatial autoregressive models. Wang and Tang [19] considered Bayesian inference on a quantile regression model in the presence of nonignorable missing covariates. Zhao et al. [20] studied the Bayesian composite quantile regression with adaptive group Lasso penalty. However, to the best of our knowledge, there is little work done for Bayesian analysis of the functional linear spatial autoregressive model. In particular, Bayesian variable selection for this model has almost never been studied. In addition, spatial dependence is a common data feature in functional analysis fields such as meteorology, econometrics, oceanography, and environmental science. Therefore, the research of the functional linear spatial autoregressive model has a strong application background. Hence, in this study, we introduce a Bayesian adaptive Lasso approach into the functional linear spatial autoregressive model on the basis of the hybrid algorithm combining the Gibbs sampler and Metropolis–Hastings algorithm, which is the attempt for this model.The outline of this study is as follows. In Section2, we introduce a partial functional linear spatial autoregressive model and then give the likelihood function by approximating to the functional coefficient using Karhunen–Loève representation. In Section 3, in order to provide a Bayesian adaptive Lasso procedure, we specify the prior distributions, derive the full conditional distributions of unknown parameters, and then describe the detailed sampling algorithms by combining the Gibbs sampler and the Metropolis–Hastings algorithm. To illustrate the proposed methodology, results obtained from some simulation studies are presented in Section 4. The study is concluded with a brief discussion in Section 5.


## 2. Model and Likelihood

### 2.1. Model

In this study, we consider the following partial functional linear spatial autoregressive model:(1)Yi=ρ∑j=1nwijYj+ZiTθ+∫01βtXitdt+εi,where Yi is a real-valued spatially dependent response variable, Zi is a p-dimensional random variables, and Xit is a zero mean, second-order stochastic process belonging to L20,1, i=1,…,n. In addition, wij is the i,jth element of a given n×n spatial weights matrix W with zero diagonal elements, i.e., wij=0, for all i=j. Besides, θ=θ1,θ2,…,θpT is a p-dimensional unknown parameters, βt is an unknown square integrable function on 0,1, and εis are mutually independent and identically distributed normal random variables with zero mean and variance σ2.Next, for convenience, we use matrix and vector notations to represent variables and models. DenoteY=Y1,Y2,…,YnT, Z=Z1,Z2,…,ZnT, Xt=X1t,X2t,…,XntT, and ε=ε1,ε2,…,εnT. Then, model (1) can be written as(2)Y=ρWY+Zθ+∫01βtXtdt+ε,where the random errors ε∼N0,σ2In and In is an n×n identity matrix.


### 2.2. Likelihood

The functional part of the model involves infinite dimensional variables and unknown parameters, so we first transform the infinite dimensional problem into a finite dimensional problem. Specifically, define the covariance function for the functional variableX and its empirical version, respectively, as Ks,t=CovXt,Xs and K^s,t=1/n∑i=1nXisXit. A linear operator is defined by the covariance function K, which maps a function f to Kf given by Kfu=∫Ku,vfvdv. It is assumed that the linear operator with kernel K is positive definite. By Mercer’s theorem, the spectral decompositions of the covariance function K and K^ can be expressed as Ks,t=∑j=1∞λ˜jϕjsϕjt and K^s,t=∑j=1∞λ˜^jϕ^jsϕ^jt, respectively, where λ˜1>λ˜2>⋯>0 and λ˜^1≥λ˜^2≥⋯≥0 are, respectively, the ordered eigenvalues of the linear operators associated with kernels K and K^ and ϕj and ϕ^j are the corresponding orthonormal eigenfunctions. Obviously, the sequences ϕj and ϕ^j each form an orthonormal basis in L20,1. Thus, according to the Karhunen–Loève representation, we have the following expansions:(3)Xt=∑i=1∞ξiϕit,βt=∑i=1∞γiϕit,where ξi are uncorrelated random variables with zero mean and variance Eξi2=λ˜i, γi=β,ϕi, and ⋅,⋅ represents inner product. More details can also be found in Ramsay and Silverman [21]. Combining (2) with (3), we rewrite model (2) as(4)Y=ρWY+Zθ+∑j=1∞γjϕj,X+ε.Thus, the model in (4) can be well approximated as follows;(5)Y≈ρWY+Zθ+∑j=1mγjϕj,X+ε,where the truncation parameter m≤n is the cut-off level and generally diverges with n. We replace ϕj with the corresponding estimator ϕ^j for j=1,…,m; thus, we can rewrite model (5) as(6)Y≈ρWY+Zθ+Uγ+ε,where U=ϕj,Xj=1,…,m and γ=γ1,…,γmT.Then, we can obtain the likelihood function of model (5):(7)Lθ,γ,ρ,σ2|Y,Z,X∝Aσ2−n2exp−eTe2σ2,where e=AY−Zθ−Uγ and A=In−ρW.


## 2.1. Model

In this study, we consider the following partial functional linear spatial autoregressive model:(1)Yi=ρ∑j=1nwijYj+ZiTθ+∫01βtXitdt+εi,where Yi is a real-valued spatially dependent response variable, Zi is a p-dimensional random variables, and Xit is a zero mean, second-order stochastic process belonging to L20,1, i=1,…,n. In addition, wij is the i,jth element of a given n×n spatial weights matrix W with zero diagonal elements, i.e., wij=0, for all i=j. Besides, θ=θ1,θ2,…,θpT is a p-dimensional unknown parameters, βt is an unknown square integrable function on 0,1, and εis are mutually independent and identically distributed normal random variables with zero mean and variance σ2.Next, for convenience, we use matrix and vector notations to represent variables and models. DenoteY=Y1,Y2,…,YnT, Z=Z1,Z2,…,ZnT, Xt=X1t,X2t,…,XntT, and ε=ε1,ε2,…,εnT. Then, model (1) can be written as(2)Y=ρWY+Zθ+∫01βtXtdt+ε,where the random errors ε∼N0,σ2In and In is an n×n identity matrix.


## 2.2. Likelihood

The functional part of the model involves infinite dimensional variables and unknown parameters, so we first transform the infinite dimensional problem into a finite dimensional problem. Specifically, define the covariance function for the functional variableX and its empirical version, respectively, as Ks,t=CovXt,Xs and K^s,t=1/n∑i=1nXisXit. A linear operator is defined by the covariance function K, which maps a function f to Kf given by Kfu=∫Ku,vfvdv. It is assumed that the linear operator with kernel K is positive definite. By Mercer’s theorem, the spectral decompositions of the covariance function K and K^ can be expressed as Ks,t=∑j=1∞λ˜jϕjsϕjt and K^s,t=∑j=1∞λ˜^jϕ^jsϕ^jt, respectively, where λ˜1>λ˜2>⋯>0 and λ˜^1≥λ˜^2≥⋯≥0 are, respectively, the ordered eigenvalues of the linear operators associated with kernels K and K^ and ϕj and ϕ^j are the corresponding orthonormal eigenfunctions. Obviously, the sequences ϕj and ϕ^j each form an orthonormal basis in L20,1. Thus, according to the Karhunen–Loève representation, we have the following expansions:(3)Xt=∑i=1∞ξiϕit,βt=∑i=1∞γiϕit,where ξi are uncorrelated random variables with zero mean and variance Eξi2=λ˜i, γi=β,ϕi, and ⋅,⋅ represents inner product. More details can also be found in Ramsay and Silverman [21]. Combining (2) with (3), we rewrite model (2) as(4)Y=ρWY+Zθ+∑j=1∞γjϕj,X+ε.Thus, the model in (4) can be well approximated as follows;(5)Y≈ρWY+Zθ+∑j=1mγjϕj,X+ε,where the truncation parameter m≤n is the cut-off level and generally diverges with n. We replace ϕj with the corresponding estimator ϕ^j for j=1,…,m; thus, we can rewrite model (5) as(6)Y≈ρWY+Zθ+Uγ+ε,where U=ϕj,Xj=1,…,m and γ=γ1,…,γmT.Then, we can obtain the likelihood function of model (5):(7)Lθ,γ,ρ,σ2|Y,Z,X∝Aσ2−n2exp−eTe2σ2,where e=AY−Zθ−Uγ and A=In−ρW.


## 3. Bayesian Adaptive Lasso

### 3.1. Priors

In this study, a Bayesian approach is implemented to estimate the unknown parametersθ, γ, ρ, and σ2. For carrying out it, we should appoint a prior distribution for unknown parameters of the model. Firstly, we need to specify the prior of θ and consider a Bayesian adaptive Lasso approach in this study which is the Bayesian version of the adaptive Lasso method. Specifically, according to the idea of Bayesian adaptive Lasso, the hierarchical priors for θ are considered as follows:(8)θ|η12,…,ηp2∼Np0,Dηfηj2=λj22expλj2ηj22pη=∏j=1pfηj2Dη=diagη12,…,ηp2,η=η12,…,ηp2T.In other words,θj|ηj2∼N0,ηj2, for j=1,…,p. The above hierarchical priors lead to the following conditionally independent double-exponential prior density for θ:(9)pθ|λ1,λ2,…,λp=∏j=1pλjexp−λjθj,which shows that θj|λj∼DEλj=λjexp−λjθj, for j=1,…,p, where DEλj denotes the double-exponential density function. According to Bayes theorem, we can obtain the posterior distribution of θ:(10)pθ|Y,Z,X,γ,ρ,σ2∝pY|Z,X,θ,γ,ρ,σ2Pθ|λ1,λ2,…,λp∝exp−12σ2AY−Zθ−UγTAY−Zθ−Uγ−∑j=1pλjθj,which implies that the posterior estimate of θ is just the adaptive Lasso estimator of θ. Generally, a big λj corresponds to unimportant variables, while a small λj corresponds to important variables. Thus, we can make variable selection through the value of λj′s. In particular, if all values of λj′s are equal, the above procedure becomes the well-known Bayesian Lasso method. To enable Bayesian analysis, in this study, we use the hierarchical Bayesian approach to select the tuning parameter λj by regarding λj as a random variable. Note that the tuning parameter λj is larger than zero, so it is natural to consider a Gamma prior. Here, we put Gamma prior on λj2 rather than λj for the convenience of sampling. Specifically, the Gamma prior for λj2 is given as follows:(11)λj2∼Γc1,d1∝λj2c1−1exp−d1λj2,where c1 and d1 are predefined hyperparameters. In addition, the priors of parameters γ and ρ are chosen as γ∼Nγ0,Bγ and ρ∼U−1,1, where the hyperparameters γ0 and Bγ are assumed to be known vectors or matrix. Besides, σ2 is taken to be Inverse Gamma (denoted as IG ( c0, d0)) with the density function given by σ2∼IGc0,d0∝σ2−c0−1exp−d0/σ2, where c0 and d0 are known positive numbers. In this study, we focus on cases that the prior distributions for the parameters of the models are listed above such as the normal distribution, inverse Gamma, and uniform distribution. However, the proposed procedure can easily be modified to other specific prior distributions. Thus, the joint prior of all the unknown parameters is given by(12)πθ,γ,ρ,σ2=pγpρpθ|λ1,λ2,…,λp∏j=1ppλj.


### 3.2. Posterior Inference

Letψ=θ,γ,ρ,σ2, and we need to estimate the unknown parameters of ψ. Based on the likelihood function (7) and priors (12), the joint posterior distribution pψ|Y,Z,X of the unknown parameters is given by(13)pψ|Y,Z,X∝Lθ,γ,ρ,σ2|Y,Z,Xπθ,γ,ρ,σ2.To make Bayesian inference based on the joint posterior distribution (13), we firstly derive the full conditional distributions of the unknown parameters and then construct the Gibbs sampler with the Metropolis–Hastings sampling algorithms to generate posterior samples from their full conditional distributions as follows.(i)
Full conditional distribution ofθ:(14)pθ|Y,Z,X,γ,σ2∝exp−12θ−μ˜θTΣ˜θ−1θ−μ˜θ,whereμ˜θ=Σ˜θZTAY−Uγ/σ2,Σ˜θ=Σξ−1+ZTZ/σ2−1, and Σξ=diagη12,η22,…,ηp2.(ii)
Full conditional distribution ofγ:(15)pγ|Y,Z,X,θ,σ2∝exp−12γ−μ˜γTΣ˜γ−1γ−μ˜γ,whereμ˜γ=Σ˜γUTAY−Zθ/σ2+Bγ−1γ0 and Σ˜γ=Bγ−1+UTU/σ2−1.(iii)
Full conditional distribution ofη2:(16)pηj−2|Y,Z,X,θ,σ2∼Inverse Gaussianμ′,λ′,whereμ′=λj/θj, and λ′=λj2 in the parameterization of inverse Gaussian density function is given by(17)fx|μ′,λ′=λ′2πx−3/2exp−λ′x−μ′22μ′2x,x>0.(iv)
Full conditional distribution ofλj2:(18)pλj2|Y,Z,θ,σ2∝λj2aλj−1exp−bλjλj2,whereaλj=c1+1 and bλj=ηj2/2+d1.(v)
Full conditional distribution ofσ2:(19)pσ2|Y,Z,X,θ,γ,ρ∝σ2−aσ2−1exp−bσ2σ2,whereaσ2=n/2+c0 and bσ2=AY−Zθ−UγTAY−Zθ−Uγ/2+d0.(vi)
Full conditional distribution ofρ:(20)pρ|Y,Z,X,θ,γ,σ2∝Aexp−12σ2AY−Zθ−UγTAY−Zθ−Uγ.Consequently, based on the conditional posterior distributions from (14) to (20), we can construct an efficient MCMC-based sampling algorithm for generating posterior samples which is summarized in Algorithm 1. It is easy to see that the conditional posterior distributions from (14) to (19) are some familiar distributions, such as the Inverse Gaussian, Inverse Gamma, Gamma, and normal distributions. Thus, it is fast and convenient to generate posterior samples from these standard distributions. However, the conditional posterior distribution of (20) has nonstandard density function and rather complicated; it is hard to directly draw observations based on this full conditional distribution. We prefer the well-known Metropolis–Hastings algorithm to overcome the difficulty. Firstly, we choose the normal distribution N0,σρ2 as the proposal distribution, where σρ2 is chosen such that the average acceptance rate is about between 0.25 and 0.45 (Gelman et al. [22]). Then, the Metropolis–Hastings algorithm is implemented as follows: at the l+1th iteration with the current value ρl, a new candidate ρ∗ is generated from Nρl,σρ2 and is accepted with probability:(21)min1,pρ∗|Y,Z,X,θ,γ,σ2pρl|Y,Z,X,θ,γ,σ2.Consequently, we could obtain posterior samplers from the algorithm in Algorithm 1 which have converged to the joint posterior distribution of (13) after the burn-in period. Then, a total number of M MCMC samples ψk=θk,γk,ρk,σ2k:k=1,2,…,M,M<J can be collected via the proposed hybrid algorithm. Thus, the Bayesian estimates of the parameters θ^,γ^,ρ^,σ^2 are defined, respectively, as follows:(22)θ^=1M∑k=1Mθk,γ^=1M∑k=1Mγk,ρ^=1M∑k=1Mρk,σ^2=1M∑k=1Mσ2k.Algorithm 1: An MCMC-based sampling algorithm for ψ=θ,γ,ρ,σ2.
Input: set up initial values ψ0=θ0,γ0,ρ0,σ20, and the number of iterations of the sampling algorithm is J.Output: a sequence of posterior samples ψ1,ψ2,…,ψJ.MCMC iterations: (the detailed Gibbs sampler cycles are outlined as follows)Fork⟵1 to J do:(1)
Sampleσ2|Y,Z,X,θ,γ,ρ∼IGaσ2,bσ2;(2)
Sampleθ|Y,Z,γ,ρ,σ2∼Nμ˜θ,Σ˜θ;(3)
Sampleγ|Y,Z,θ,ρ,σ2∼Nμ˜γ,Σ˜γ;(4)
Sampleλj2|Y,Z,X,θ,γ,ρ,σ2∼Γaλj,bλj;(5)
Sampleηj−2|Y,Z,X,θ,γ,ρ,σ2∼ Inverse Gaussian μ′,λ′;(6)
Sampleρ|Y,Z,X,θ,γ,σ2 from (20) based on the Metropolis–Hastings algorithm.End


## 3.1. Priors

In this study, a Bayesian approach is implemented to estimate the unknown parametersθ, γ, ρ, and σ2. For carrying out it, we should appoint a prior distribution for unknown parameters of the model. Firstly, we need to specify the prior of θ and consider a Bayesian adaptive Lasso approach in this study which is the Bayesian version of the adaptive Lasso method. Specifically, according to the idea of Bayesian adaptive Lasso, the hierarchical priors for θ are considered as follows:(8)θ|η12,…,ηp2∼Np0,Dηfηj2=λj22expλj2ηj22pη=∏j=1pfηj2Dη=diagη12,…,ηp2,η=η12,…,ηp2T.In other words,θj|ηj2∼N0,ηj2, for j=1,…,p. The above hierarchical priors lead to the following conditionally independent double-exponential prior density for θ:(9)pθ|λ1,λ2,…,λp=∏j=1pλjexp−λjθj,which shows that θj|λj∼DEλj=λjexp−λjθj, for j=1,…,p, where DEλj denotes the double-exponential density function. According to Bayes theorem, we can obtain the posterior distribution of θ:(10)pθ|Y,Z,X,γ,ρ,σ2∝pY|Z,X,θ,γ,ρ,σ2Pθ|λ1,λ2,…,λp∝exp−12σ2AY−Zθ−UγTAY−Zθ−Uγ−∑j=1pλjθj,which implies that the posterior estimate of θ is just the adaptive Lasso estimator of θ. Generally, a big λj corresponds to unimportant variables, while a small λj corresponds to important variables. Thus, we can make variable selection through the value of λj′s. In particular, if all values of λj′s are equal, the above procedure becomes the well-known Bayesian Lasso method. To enable Bayesian analysis, in this study, we use the hierarchical Bayesian approach to select the tuning parameter λj by regarding λj as a random variable. Note that the tuning parameter λj is larger than zero, so it is natural to consider a Gamma prior. Here, we put Gamma prior on λj2 rather than λj for the convenience of sampling. Specifically, the Gamma prior for λj2 is given as follows:(11)λj2∼Γc1,d1∝λj2c1−1exp−d1λj2,where c1 and d1 are predefined hyperparameters. In addition, the priors of parameters γ and ρ are chosen as γ∼Nγ0,Bγ and ρ∼U−1,1, where the hyperparameters γ0 and Bγ are assumed to be known vectors or matrix. Besides, σ2 is taken to be Inverse Gamma (denoted as IG ( c0, d0)) with the density function given by σ2∼IGc0,d0∝σ2−c0−1exp−d0/σ2, where c0 and d0 are known positive numbers. In this study, we focus on cases that the prior distributions for the parameters of the models are listed above such as the normal distribution, inverse Gamma, and uniform distribution. However, the proposed procedure can easily be modified to other specific prior distributions. Thus, the joint prior of all the unknown parameters is given by(12)πθ,γ,ρ,σ2=pγpρpθ|λ1,λ2,…,λp∏j=1ppλj.


## 3.2. Posterior Inference

Letψ=θ,γ,ρ,σ2, and we need to estimate the unknown parameters of ψ. Based on the likelihood function (7) and priors (12), the joint posterior distribution pψ|Y,Z,X of the unknown parameters is given by(13)pψ|Y,Z,X∝Lθ,γ,ρ,σ2|Y,Z,Xπθ,γ,ρ,σ2.To make Bayesian inference based on the joint posterior distribution (13), we firstly derive the full conditional distributions of the unknown parameters and then construct the Gibbs sampler with the Metropolis–Hastings sampling algorithms to generate posterior samples from their full conditional distributions as follows.(i)
Full conditional distribution ofθ:(14)pθ|Y,Z,X,γ,σ2∝exp−12θ−μ˜θTΣ˜θ−1θ−μ˜θ,whereμ˜θ=Σ˜θZTAY−Uγ/σ2,Σ˜θ=Σξ−1+ZTZ/σ2−1, and Σξ=diagη12,η22,…,ηp2.(ii)
Full conditional distribution ofγ:(15)pγ|Y,Z,X,θ,σ2∝exp−12γ−μ˜γTΣ˜γ−1γ−μ˜γ,whereμ˜γ=Σ˜γUTAY−Zθ/σ2+Bγ−1γ0 and Σ˜γ=Bγ−1+UTU/σ2−1.(iii)
Full conditional distribution ofη2:(16)pηj−2|Y,Z,X,θ,σ2∼Inverse Gaussianμ′,λ′,whereμ′=λj/θj, and λ′=λj2 in the parameterization of inverse Gaussian density function is given by(17)fx|μ′,λ′=λ′2πx−3/2exp−λ′x−μ′22μ′2x,x>0.(iv)
Full conditional distribution ofλj2:(18)pλj2|Y,Z,θ,σ2∝λj2aλj−1exp−bλjλj2,whereaλj=c1+1 and bλj=ηj2/2+d1.(v)
Full conditional distribution ofσ2:(19)pσ2|Y,Z,X,θ,γ,ρ∝σ2−aσ2−1exp−bσ2σ2,whereaσ2=n/2+c0 and bσ2=AY−Zθ−UγTAY−Zθ−Uγ/2+d0.(vi)
Full conditional distribution ofρ:(20)pρ|Y,Z,X,θ,γ,σ2∝Aexp−12σ2AY−Zθ−UγTAY−Zθ−Uγ.Consequently, based on the conditional posterior distributions from (14) to (20), we can construct an efficient MCMC-based sampling algorithm for generating posterior samples which is summarized in Algorithm 1. It is easy to see that the conditional posterior distributions from (14) to (19) are some familiar distributions, such as the Inverse Gaussian, Inverse Gamma, Gamma, and normal distributions. Thus, it is fast and convenient to generate posterior samples from these standard distributions. However, the conditional posterior distribution of (20) has nonstandard density function and rather complicated; it is hard to directly draw observations based on this full conditional distribution. We prefer the well-known Metropolis–Hastings algorithm to overcome the difficulty. Firstly, we choose the normal distribution N0,σρ2 as the proposal distribution, where σρ2 is chosen such that the average acceptance rate is about between 0.25 and 0.45 (Gelman et al. [22]). Then, the Metropolis–Hastings algorithm is implemented as follows: at the l+1th iteration with the current value ρl, a new candidate ρ∗ is generated from Nρl,σρ2 and is accepted with probability:(21)min1,pρ∗|Y,Z,X,θ,γ,σ2pρl|Y,Z,X,θ,γ,σ2.Consequently, we could obtain posterior samplers from the algorithm in Algorithm 1 which have converged to the joint posterior distribution of (13) after the burn-in period. Then, a total number of M MCMC samples ψk=θk,γk,ρk,σ2k:k=1,2,…,M,M<J can be collected via the proposed hybrid algorithm. Thus, the Bayesian estimates of the parameters θ^,γ^,ρ^,σ^2 are defined, respectively, as follows:(22)θ^=1M∑k=1Mθk,γ^=1M∑k=1Mγk,ρ^=1M∑k=1Mρk,σ^2=1M∑k=1Mσ2k.Algorithm 1: An MCMC-based sampling algorithm for ψ=θ,γ,ρ,σ2.
Input: set up initial values ψ0=θ0,γ0,ρ0,σ20, and the number of iterations of the sampling algorithm is J.Output: a sequence of posterior samples ψ1,ψ2,…,ψJ.MCMC iterations: (the detailed Gibbs sampler cycles are outlined as follows)Fork⟵1 to J do:(1)
Sampleσ2|Y,Z,X,θ,γ,ρ∼IGaσ2,bσ2;(2)
Sampleθ|Y,Z,γ,ρ,σ2∼Nμ˜θ,Σ˜θ;(3)
Sampleγ|Y,Z,θ,ρ,σ2∼Nμ˜γ,Σ˜γ;(4)
Sampleλj2|Y,Z,X,θ,γ,ρ,σ2∼Γaλj,bλj;(5)
Sampleηj−2|Y,Z,X,θ,γ,ρ,σ2∼ Inverse Gaussian μ′,λ′;(6)
Sampleρ|Y,Z,X,θ,γ,σ2 from (20) based on the Metropolis–Hastings algorithm.End


## 4. Simulation Study

In this section, we investigate the performance of the proposed model and Bayesian estimation method via Monte Carlo simulation. The datasets are generated from the following partial functional linear spatial autoregressive model:(23)Y=ρWY+Zθ+∫01βtXtdt+ε,where Z follows a multivariate normal distribution N0,Σz and θ=1,−1,0,0,1,0,0,0, here Σzi,j=0.5i−j, for i,j=1,2,…,8. To have a better comparison, we consider two different values σ2=0.25 and 1, which represent strong and weak signal-to-noise ratios. Furthermore, we take the spatial parameter ρ=0.5,0,−0.5, which describes different spatial dependencies and dependence directions. Similar to Xie et al. [11], the weight matrix is set to be W=IR⊗Hq,Hq=lqlqT−Iq/q−1, where lq is an q-dimensional vector with all elements being 1 and ⊗ means Kronecker product. For the functional part, we take the same form as Shin [6]; in other words, the functional coefficient βt=2sinπt/2+32sin3πt/2 and Xt=∑j=150ξjϕjt, where ξj are distributed as independent normal with mean zero and variance λ˜j=j−0.5π−2 and ϕjt=2sinj−0.5πt. In the simulation, we consider the noninformative prior information type of hyperparameter values for unknown parameters θ,γ,ρ,σ2: γ0=0m,Bγ=Im,c0=d0=0.01,c1=1,andd1=0.1, where 0m is an m-dimensional vector with all elements being 0. In addition, we choose R to be 50 and 100 and q to be 3; thus, the sample size n=R×q equals to 150 and 300. Throughout our numerical studies, we choose the truncation parameter m such that the first m functional principal components’ scores can explain at least 90% of the total variability of the functional predictor Xt.Based on the generated dataset and the above simulation settings, we use the previously proposed MCMC-based sampling algorithm in Algorithm 1 to obtain the Bayesian estimates of unknown parameters based on 100 replications. To investigate the convergence of the proposed MCMC algorithm, we compute the EPSR values (estimated potential scale reduction) [23] for a number of test runs on the base of three parallel chains of observations via three different starting values. We find that, in all test runs, the EPSR values are close to 1 and less than 1.2 after discarding the first 3,000 as the burn-in period. So, we can collect the observations after 3000 iterations with M = 2000 in producing the Bayesian estimates for each replication and posterior summary of the parameters is listed in Tables 1–3. In addition, the accuracy of the functional coefficient estimator is checked by the square root of average square errors (RASE) which is defined as follows:(24)RASE=1N∑s=1Nβ^ts−βts21/2,where ts,s=1,2,…,N, are the grid points at which the functional coefficient estimator β^t is evaluated and N = 200 is used in our simulation. The simulation results of the mean and standard error (SE) of RASE for the functional component under different cases are reported in Table 4. Moreover, to see accuracy of the estimated function βt directly, we plot the true value of functional coefficient βt against its estimate under different cases, and list these functional curve fitting results under different spatial parameters in Figures 1–4.Table 1
Bayesian estimates of unknown parameters under different cases whenρ=0.5.
σ2R,qParaθ1θ2θ3θ4θ5θ6θ7θ8ρσ21(50, 3)Bias−0.0250.032−0.0050.012−0.0230.014−0.0080.0030.0000.016SD0.0970.1030.0900.0770.1100.0880.0850.0750.0340.140PP0.000.000.980.980.000.980.990.960.000.00(100, 3)Bias−0.0220.022−0.0160.001−0.0220.0160.001−0.009−0.0030.010SD0.0750.0780.0670.0690.0660.0610.0590.0570.0260.093PP0.000.000.970.950.000.991.000.940.000.000.25(50, 3)Bias−0.0060.009−0.0010.005−0.0070.006−0.0060.002−0.0010.004SD0.0460.0490.0510.0440.0550.0500.0500.0430.0200.034PP0.000.000.960.980.000.980.960.950.000.00(100, 3)Bias−0.0090.008−0.0080.000−0.0090.0080.001−0.005−0.0010.003SD0.0360.0390.0370.0380.0330.0330.0330.0310.0160.023PP0.000.000.960.950.000.990.970.940.000.00Table 2
Bayesian estimates of unknown parameters under different cases whenρ=0.
σ2R,qParaθ1θ2θ3θ4θ5θ6θ7θ8ρσ21(50, 3)Bias−0.0250.032−0.0050.012−0.0230.014−0.0080.0030.0000.016SD0.0960.1020.0900.0770.1080.0880.0850.0750.0530.136PP0.000.000.980.980.000.980.980.970.950.00(100, 3)Bias−0.0240.023−0.0160.001−0.0230.0160.001−0.009−0.0050.007SD0.0750.0770.0680.0690.0660.0600.0590.0570.0420.090PP0.000.000.980.960.000.991.000.940.920.000.25(50, 3)Bias−0.0060.009−0.0010.005−0.0070.007−0.0060.002−0.0010.004SD0.0450.0490.0510.0440.0540.0500.0500.0430.0310.034PP0.000.000.960.980.000.980.970.960.980.00(100, 3)Bias−0.0090.009−0.0080.000−0.0090.0080.001−0.005−0.0020.002SD0.0360.0380.0370.0380.0330.0330.0330.0310.0250.022PP0.000.000.960.950.000.990.990.940.920.00Table 3
Bayesian estimates of unknown parameters under different cases whenρ=−0.5.
σ2R,qParaθ1θ2θ3θ4θ5θ6θ7θ8ρσ21(50, 3)Bias−0.0250.031−0.0050.012−0.0230.015−0.0090.0030.0000.017SD0.0970.1030.0900.0770.1090.0890.0860.0760.0590.138PP0.000.000.970.980.000.980.980.970.000.00(100, 3)Bias−0.0250.025−0.0160.001−0.0240.0160.002−0.010−0.0070.004SD0.0770.0780.0670.0690.0670.0600.0590.0570.0470.092PP0.000.000.970.950.000.981.000.940.000.000.25(50, 3)Bias−0.0070.009−0.0010.005−0.0070.007−0.0060.002−0.0020.004SD0.0460.0500.0510.0440.0550.0500.0500.0430.0350.033PP0.000.000.950.980.000.980.960.950.000.00(100, 3)Bias−0.0100.009−0.0080.000−0.0100.0080.001−0.005−0.0020.002SD0.0380.0380.0370.0380.0340.0330.0330.0310.0280.022PP0.000.000.960.950.000.980.980.940.000.00Table 4
The mean and SE of RASE for the functional coefficientβt under different cases.
ρσ2R,qMeanSE0.51(50, 3)0.6420.314(100, 3)0.3940.1730.25(50, 3)0.3900.170(100, 3)0.2610.09101(50, 3)0.6420.313(100, 3)0.3940.1720.25(50, 3)0.3910.170(100, 3)0.2620.089−0.51(50, 3)0.6430.312(100, 3)0.3960.1730.25(50, 3)0.3920.170(100, 3)0.2640.088Figure 1
The estimated functional coefficient versus true value ofβt with different ρs (from left to right responding to ρ=0.5,0,−0.5) when R=50,q=3,andσ2=0.25.Figure 2
The estimated functional coefficient versus true value ofβt with different ρs (from left to right responding to ρ=0.5,0,−0.5) when R=100,q=3,andσ2=0.25.Figure 3
The estimated functional coefficient versus true value ofβt with different ρs (from left to right responding to ρ=0.5,0,−0.5) when R=50,q=3,andσ2=1.Figure 4
The estimated functional coefficient versus true value ofβt with different ρs (from left to right responding to ρ=0.5,0,−0.5) when R=100,q=3,andσ2=1.In Tables1–3, “Bias” means the difference between the true value and the average of the Bayesian estimates of the parameters based on 100 replications, “SD” means standard deviation of the Bayesian estimates, and “PP” denotes the proportion that parameter was identified to be zero in 100 replications in terms of the criterion that a parameter was identified to be 0 if its 95% confidence interval contains zero. By careful observation of Tables 1–4, several conclusions are obtained and summarized as follows. (1) Under all the considered settings, the performance of the proposed Bayesian method is quite satisfactory in the sense that their values of Bias and SD are reasonable small. Furthermore, the results of Bayesian estimation become better and better as the sample size increases. (2) As expected, the performance of Bayesian estimation is similar under different spatial parameters. (3) As the variance σ2 becomes smaller, the Bayesian estimation method performs better. (4) As we predicted, the proposed Bayesian method could identify the correct models in most cases because the proportion values in tables corresponding to the important covariates are all zero, and PP values corresponding to unimportant covariates are more than 90%. (5) For the functional component, the mean values and SE values of RASE become smaller and smaller as the sample size increases. This indicates that the functional coefficient βt is estimated better and better. In addition, examination of Figures 1–4 shows that the shapes of the estimated functional coefficient are close to the corresponding true line under all the considered settings, which agrees with what was discovered from Table 4. To sum up, all the above findings show that the preceding proposed Bayesian estimation procedure and MCMC-based sampling algorithm can well recover the true information in partial functional linear spatial autoregressive model.


## 5. Conclusion and Discussion

In this study, we propose a partial functional linear spatial autoregressive model which can explore the relationship between a scalar spatially dependent response variable and predictive variables containing both multiple scalar covariates and a functional covariate. Based on functional principal components analysis, a Bayesian adaptive Lasso approach is developed to analyse the model by combining the Gibbs sampler and Metropolis–Hastings algorithm. We conduct some simulation studies to show the efficiency of the proposed Bayesian method. The results indicate that the developed Bayesian variable selection method is highly efficient and computationally convenient.Furthermore, there are several interesting extensions that can be considered in the future. Specific extensions are considered as follows. (i) Compared with classical mean regression, quantile regression analysis should be more robust to nonnormal errors and outliers. Therefore, using quantile regression technology to analyse functional data with spatial dependence warrants a future investigation. (ii) Missing data analysis has always been a statistical hotspot. Therefore, it is also worth studying of the current model when covariates are missing under different missingness mechanism. We leave these topics for future research.



---

*Source: 1616068-2022-04-15.xml*