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.
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 Section 2, 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:where is a real-valued spatially dependent response variable, is a -dimensional random variables, and is a zero mean, second-order stochastic process belonging to , . In addition, is the th element of a given spatial weights matrix with zero diagonal elements, i.e., , for all Besides, is a -dimensional unknown parameters, is an unknown square integrable function on , and s are mutually independent and identically distributed normal random variables with zero mean and variance .
Next, for convenience, we use matrix and vector notations to represent variables and models. Denote , , and . Then, model (1) can be written aswhere the random errors and is an 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 variable and its empirical version, respectively, as and A linear operator is defined by the covariance function , which maps a function to given by . It is assumed that the linear operator with kernel is positive definite. By Mercer’s theorem, the spectral decompositions of the covariance function and can be expressed as and , respectively, where and are, respectively, the ordered eigenvalues of the linear operators associated with kernels and and and are the corresponding orthonormal eigenfunctions. Obviously, the sequences and each form an orthonormal basis in . Thus, according to the Karhunen–Loève representation, we have the following expansions:where are uncorrelated random variables with zero mean and variance , , and represents inner product. More details can also be found in Ramsay and Silverman [21]. Combining (2) with (3), we rewrite model (2) as
Thus, the model in (4) can be well approximated as follows;where the truncation parameter is the cut-off level and generally diverges with . We replace with the corresponding estimator for ; thus, we can rewrite model (5) aswhere and .
Then, we can obtain the likelihood function of model (5):where and .
3. Bayesian Adaptive Lasso
3.1. Priors
In this study, a Bayesian approach is implemented to estimate the unknown parameters , , , and . 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:
In other words, , for . The above hierarchical priors lead to the following conditionally independent double-exponential prior density for :which shows that , for , where denotes the double-exponential density function. According to Bayes theorem, we can obtain the posterior distribution of :which implies that the posterior estimate of is just the adaptive Lasso estimator of . Generally, a big corresponds to unimportant variables, while a small corresponds to important variables. Thus, we can make variable selection through the value of . In particular, if all values of 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 by regarding as a random variable. Note that the tuning parameter is larger than zero, so it is natural to consider a Gamma prior. Here, we put Gamma prior on rather than for the convenience of sampling. Specifically, the Gamma prior for is given as follows:where and are predefined hyperparameters. In addition, the priors of parameters and are chosen as and , where the hyperparameters and are assumed to be known vectors or matrix. Besides, is taken to be Inverse Gamma (denoted as IG ( , )) with the density function given by , where and 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
3.2. Posterior Inference
Let , and we need to estimate the unknown parameters of . Based on the likelihood function (7) and priors (12), the joint posterior distribution of the unknown parameters is given by
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 : where , and .(ii)Full conditional distribution of : where and .(iii)Full conditional distribution of : where and in the parameterization of inverse Gaussian density function is given by(iv)Full conditional distribution of : where and (v)Full conditional distribution of : where and (vi)Full conditional distribution of :
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 as the proposal distribution, where 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 th iteration with the current value , a new candidate is generated from and is accepted with probability:
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 can be collected via the proposed hybrid algorithm. Thus, the Bayesian estimates of the parameters are defined, respectively, as follows:
| ||||||||||||||||||||||
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:where follows a multivariate normal distribution and , here , for To have a better comparison, we consider two different values and 1, which represent strong and weak signal-to-noise ratios. Furthermore, we take the spatial parameter , which describes different spatial dependencies and dependence directions. Similar to Xie et al. [11], the weight matrix is set to be , where 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 and , where are distributed as independent normal with mean zero and variance and . In the simulation, we consider the noninformative prior information type of hyperparameter values for unknown parameters : where is an -dimensional vector with all elements being 0. In addition, we choose to be 50 and 100 and to be 3; thus, the sample size equals to 150 and 300. Throughout our numerical studies, we choose the truncation parameter such that the first functional principal components’ scores can explain at least 90% of the total variability of the functional predictor .
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:where , are the grid points at which the functional coefficient estimator 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 directly, we plot the true value of functional coefficient against its estimate under different cases, and list these functional curve fitting results under different spatial parameters in Figures 1–4.
In Tables 1–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 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 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.
Data Availability
No data were used to support this study, and computer code is available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This work was supported by the Major Humanities and Social Sciences Tackling Project in Universities in Zhejiang Province (2018QN037).