Sparse Boosting for Additive Spatial Autoregressive Model with High Dimensionality

1. Introduction

Spatial data, also known as geospatial data, commonly appear in fields such as environmental health, economics, and epidemiology. This type of data can be represented by numerical values in a geographic coordinate system. The spatial autoregressive (SAR) model, which includes a spatially lagged term of the response variable to account for spatial dependence among spatial units, has become an important research topic in recent econometrics and statistics literature. Significant developments in the estimation and inference of this model can be found in the works of Cliff and Ord [1], Anselin [2], Robert [3], among others.

The linear SAR model, which extends the ordinary linear regression model by incorporating a spatially lagged term of the response variable, has been extensively studied [4,5,6,7,8]. However, it imposes assumptions that might be unrealistic in practice, rendering it inefficient if not correctly specified. To enhance model flexibility and adaptability, the nonparametric SAR model allows for investigating the relationship between the response variable and predictors without assuming a specific shape for the relationship. This model has recently gained much attention from econometricians and statisticians. One of the most powerful and useful models in spatial statistics is the partial linear SAR model, where the spatially lagged response variable and some predictors enter the model linearly, while the remaining predictors are incorporated nonparametrically. Specifically, Su and Jin [9] proposed a profile quasi-maximum likelihood estimation method for the model. Koch and Krisztin [10] proposed an estimation method based on B-splines and genetic algorithms. Chen et al. [11] proposed a two-step Bayesian approach based on kernel estimation and Bayesian methods for inference. Krisztin [12] proposed a novel Bayesian semiparametric estimation method combining penalized splines with Bayesian methods. Li and Mei [13] proposed a method to test linear constraints on the parameters of the partially linear SAR model.

Beyond the partial linear SAR model, Su [14] studied the SAR model, where the spatially lagged response variable enters the model linearly and all predictors enter nonparametrically, proposing semiparametric GMM estimation under weak moment conditions. Wei et al. [15] studied SAR models with varying coefficients to capture heterogeneous effects of covariates and spatial interaction. Du et al. [16] considered a class of partially linear additive SAR models and proposed a generalized method of moments estimator. Cheng and Chen [17] studied partially linear single-index SAR models and proposed profile maximum likelihood estimators. In this paper, we will consider the additive SAR model, which is more flexible than the SAR model examined in [14], as it includes only one nonparametric component, easing the curse of dimensionality issue. In particular, we will address high-dimensional data settings where the dimensionality of predictors may exceed the sample size.

Data science is an ever-expanding field. High-dimensional models, where the feature dimension grows exponentially or non-polynomially fast with sample size, have become a focus in statistical literature. Including irrelevant features in the model can lead to undesirable computational issues and unstable estimation. To address these challenges, variable screening and selection techniques have been developed. Among these developments, penalized approaches such as lasso [18], SCAD [19], MCP [20], and their various extensions [21,22,23] have been thoroughly studied to identify important covariates and estimate the coefficients of interest simultaneously, thereby improving the prediction accuracy and interpretability of statistical models. For example, following the idea of group lasso [23] and adaptive lasso [22], Wei et al. [15] proposed a local linear shrinkage estimator for the SAR model with varying coefficients for model selection. Nevertheless, these regularization methods have one major disadvantage: they all involve tuning parameters that must be chosen using computationally intensive methods such as cross-validation.

In recent decades, boosting methods have become an effective alternative tool for high-dimensional data settings to perform variable selection and model estimation. They offer several advantages, such as relatively smaller computational cost, lower risk of overfitting, and simpler adjustments to incorporate additional constraints. Boosting was initially conceptualized as a machine learning algorithm that constructs better base learners to minimize the loss function with every iteration. The original version of the boosting algorithm proposed by Schapire [24] and Freund [25] did not fully exploit the potential of base learners. This changed when Freund and Schapire [26] proposed the AdaBoost algorithm, which could adapt to base learners more effectively. Following the development of this algorithm, many other boosting algorithms have been formulated, with major variations in their loss functions. For instance, AdaBoost with the exponential loss, L2Boosting [27] with the squared error loss, SparseL2Boosting [28] with the penalized loss, and HingeBoost [29] with the weighted hinge loss. Other versions of boosting algorithms proposed recently include pAUCBoost [30], Twin Boosting [31], Twin HingeBoost [29], ER-Boost [32], GSBoosting [33], among others. As demonstrated by Yue et al. [34,35,36,37], sparse boosting achieves better variable selection performance than L2Boosting as well as some penalized methods such as lasso. However, the application of (sparse) boosting approaches in high-dimensional SAR models has not yet been investigated.

However, there currently exists a lack of exploration regarding the application of (sparse) boosting methods in the study of high-dimensional spatial autoregressive (SAR) models, which provides an opportunity for further investigation in this research. This paper focuses on additive spatial autoregressive models with high-dimensional covariates and proposes a novel multi-step sparse boosting algorithm aimed at model-based prediction and variable selection, rather than relying on traditional regularization methods. A significant advantage of this new method is that it eliminates the need for time-consuming tuning parameter selection, thereby enhancing the convenience of modeling. Our approach is designed specifically for high-dimensional additive spatial autoregressive models, improving the flexibility and adaptability of the model in complex situations. Extending existing boosting techniques to such complex model structures presents a considerable technical challenge and requires deeper exploration both theoretically and practically. The detailed research presented in this paper facilitates the first application of this method in the field, thereby providing a foundation for future research. Through simulation studies and real data examples, we demonstrate the superior performance of the proposed method compared to various alternative algorithms. Our results not only fill the existing research gap in the literature but also offer empirical support for the multi-step sparse boosting method as an effective tool for high-dimensional data analysis, especially in scenarios involving complex relationships among variables.

The rest of the paper is organized as follows. In Section 2, the additive autoregressive model is formulated, and a multi-step sparse boosting algorithm is proposed. In Section 3, simulation studies are conducted to demonstrate the validity of this multi-step method. In Section 4, the performance of the multi-step sparse boosting method is evaluated by analyzing Boston housing price data. Concluding remarks are given in Section 5.

2. Methodology

2.1. Model and Estimation

Consider the following additive spatial autoregressive model:

$$Y=\rho WY+\sum _{j=1}^p{g}_j\left(X_j\right)+ϵ,$$

(1)

where $Y={(y_1,\cdots ,y_n)}^T$ is the $n\times 1$ vector of the observations of the response variable with n being the number of spatial units, W is an $n\times n$ spatial weight matrix of known constants with zero diagonal elements and $WY$ is referred to as spatial lag of Y, $\rho$ is the spatial autoregressive parameter with $\left|\rho \right|<1$, $X_j={(X_{1j},\cdots ,X_{nj})}^T$ is the $n\times 1$ vector of the j-th regressor and $X=(X_1,\cdots ,X_p)$ is the $n\times p$ observed matrix of regressors. For simplicity, the covariates $X_1,\cdots ,X_p$ are assumed to be distributed on compact intervals $[a_1,b_1],\cdots ,[a_p,b_p]$ respectively. $g_j(.)s$ are unknown smooth functions on $[a_j,b_j]$ with the assumption that $E\left[g_j\left(X_j\right)\right]=0$ for identifiability purpose. $ϵ={({ϵ}_1,\cdots ,{ϵ}_n)}^T$ is an $n\times 1$ vector of i.i.d disturbances with zero mean and finite variance ${\sigma }^2$.

The smoothing splines technique is typically used to estimate the unknown functions, demonstrating steady performance in practice. In this paper, we will use B-spline basis functions to approximate the unknown coefficient functions $g_j(.),j=1,\cdots ,p$. Let $B(.)=(B_1(.),\cdots ,B_L(.))$ be an equally spaced B-spline basis, where L is the dimension of the basis. Similar to Huang et al. [38], let ${\overline{B}}_{lj}={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{B}_l\left(X_{nj}\right)$, ${\pi }_{lj}(.)=B_l(.)-{\overline{B}}_{lj}$, ${\pi }_j(.)=({\pi }_{1j}(.),\cdots ,{\pi }_{Lj}(.))$. Then under appropriate smoothness assumptions, $g_j(.)\approx {\sum }_{l=1}^L{\pi }_{lj}(.){\alpha }_{lj}={\pi }_j(.){\alpha }_j$ for $j=1,\cdots ,p$. Then the spline estimator of $g(.)={\sum }_{j=1}^p{g}_j(.)$ is $\widehat{g}(.)=\pi (.)\alpha$, where $\pi (.)=({\pi }_1(.),\cdots ,{\pi }_p(.))$, $\alpha ={({\alpha }_1^T,\cdots ,{\alpha }_p^T)}^T$. It is obvious that ${\sum }_{i=1}^n{\widehat{g}}_j\left(X_{ij}\right)=0$ for $j=1,\cdots ,p$.

Then the model (1) can be rewritten as follows

$$\begin{array}{cc}\hfill Y& \approx \rho WY+\sum _{j=1}^p{\pi }_j\left(X_j\right){\alpha }_j+ϵ\hfill \\ \hfill & =\rho WY+\pi \alpha +ϵ,\hfill \end{array}$$

(2)

where $\pi =({\pi }_1\left(X_1\right),\cdots ,{\pi }_p\left(X_d\right)$. Let $P=\pi {\left({\pi }^T\pi \right)}^{-1}{\pi }^T$ denote the projection matrix onto the space spanned by $\pi$. Similar to [39], partialing out the B-spline approximation, we obtain

$$(I-P)Y\approx (I-P)\rho WY+(I-P)ϵ.$$

(3)

In the above equation, a problem of endogeneity emerges because the spatially lagged value of Y is correlated with the stochastic disturbance. Suppose H is a relevant vector instrument to eliminate the endogeneity of $WY$. Let $D=WY$. Regressing D on instrument H via OLS produces the fitted value for the endogenous variable D:

$$\widehat{D}=H{\left(H^{T}H\right)}^{-1}{H}^{T}D.$$

(4)

Substituting D with its predicted value in model (3), it becomes

$$\begin{array}{cc}\hfill (I-P)Y& \approx (I-P)\rho H{\left(H^{T}H\right)}^{-1}{H}^{T}D+(I-P)ϵ\hfill \\ \hfill & =(I-P)\rho MD+(I-P)ϵ,\hfill \end{array}$$

(5)

where $M=H{\left(H^{T}H\right)}^{-1}{H}^T$. Then by OLS,

$$\widehat{\rho }={\left(D^{T}M(I-P)MD\right)}^{-1}{D}^{T}M(I-P)Y.$$

(6)

Substituting D with $\widehat{D}$ and $\rho$ with $\widehat{\rho }$ in model (2), it becomes

$$Y\approx \widehat{\rho }D+\pi \alpha +ϵ.$$

(7)

Then the least squares loss function is close to

$${(Y-\widehat{\rho }D-\pi \alpha )}^T(Y-\widehat{\rho }D-\pi \alpha ).$$

(8)

However, when dimensionality of $\pi$ is larger than sample size n, least square estimation fails. In this case, we will adopt sparse boosting approach to estimate $\alpha$. Denote $\widehat{\alpha }$ as the estimator of $\alpha$ obtained through sparse boosting, using the squared loss function (8) as the loss function. The detailed sparse boosting algorithm will be given in the next subsection.

Consequently, $g_j(.)$ can be estimated by ${\widehat{g}}_j(.)={\pi }_j(.){\widehat{\alpha }}_j$ for $j=1,\cdots ,p$. The variance ${\sigma }^2$ can be estimated by

$${\widehat{\sigma }}^2={\displaystyle \frac{1}{n}}{(Y-\widehat{\rho }WY-\pi \widehat{\alpha })}^T(Y-\widehat{\rho }WY-\pi \widehat{\alpha }).$$

(9)

Similar to [39], we use an analogous rationality for the construction of instrument variables. In the first step, we simply regress Y on pseudo regressors variables $WY$ and $\pi$. Then the least squares loss function is

$${(Y-\rho WY-\pi \alpha )}^T(Y-\rho WY-\pi \alpha )$$

(10)

Denote $(\tilde{\rho },\tilde{\alpha })$ as the estimators of $(\rho ,\alpha )$ obtained through sparse boosting, using the squared loss function (10) as the loss function. Then, we can use the following instrumental variable:

$$\tilde{H}=W{(I-\tilde{\rho }W)}^{-1}\pi \tilde{\alpha }.$$

(11)

In the second step, instrumental variable $\tilde{H}$ is used to obtain the estimators $\overline{\alpha }$ and $\overline{\rho }$, which are then used to construct the instrumental variable

$$H=W{(I-\overline{\rho }W)}^{-1}\pi \overline{\alpha }.$$

(12)

Finally, use the instrumental variables H to obtain the final estimators $\widehat{\rho }$ and $\widehat{\alpha }={({\widehat{\alpha }}_1^T,\cdots ,{\widehat{\alpha }}_p^T)}^T$. The function $g_j(.)$ can be estimated by ${\widehat{g}}_j(.)={\widehat{\pi }}_j(.){\widehat{\alpha }}_j$ for $j=1,\cdots ,p$ and the response Y can be estimated by $\widehat{Y}={(I-\widehat{\rho }W)}^{-1}\pi \widehat{\alpha }$.

2.2. Sparse Boosting Techniques

Sparse boosting can be viewed as iteratively pursuing gradient descending in function space using a penalized empirical risk function that integrates squared loss and the complexity of the boosting measure. Similar to Yue et al. ([34,35,36,37]), we will adopt the g-prior minimum description length (gMDL) [40], a combination of squared loss and the trace of boosting operator, as the penalized empirical risk function to estimate the update criterion in each iteration and the stopping criterion. We use it because it has a data driven penalty to avoid the selection of the tuning parameter. To facilitate presentation, suppose the vector $Y^{*}$ is regressed on the $p^{*}$-dimensional matrix $X^{*}=(X_1^{*},\cdots ,X_p^{*})$. Then the least squares loss function involved in sparse boosting is ${(Y^{*}-X^{*}{\theta }^{*})}^T(Y^{*}-X^{*}{\theta }^{*})$, where ${\theta }^{*}={({\left({\theta }_1^{*}\right)}^T,\cdots ,{\left({\theta }_p^{*}\right)}^T)}^T$. The gMDL takes the form:

$$\begin{array}{cc}\hfill \mathrm{gMDL}(Y^{*},RSS,\mathrm{trace}\left(\mathcal{B}\right))=& \log \left(F\right)+{\displaystyle \frac{\mathrm{trace}\left(\mathcal{B}\right)}{n}}\log \left({\displaystyle \frac{Y^{*T}{Y}^{*}-RSS}{\mathrm{trace}\left(\mathcal{B}\right)\times F}}\right),\hfill \\ \hfill & F={\displaystyle \frac{RSS}{n-\mathrm{trace}\left(\mathcal{B}\right)}},\hfill \end{array}$$

(13)

where $RSS$ is the residual sum of squares and $\mathcal{B}$ is the boosting operator. The model achieve shortest description of data will be chosen.

We present the sparse boosting approach more specifically. The initial value of ${\theta }^{*}$ is set to the zero vector, i.e. ${\theta }^{*\left[0\right]}=0$. In each of the kth iteration ($0<k\le K$, and K is the maximum number of iterations considered in the first step), we use the residual $R^{\left[k\right]}=Y^{*}-X^{*}{\theta }^{*[k-1]}$ from the current iteration to fit each of the j-th component $X_j^{*},j=1,\cdots ,p^{*}$. The fit, denoted by ${\widehat{\lambda }}_j^{\left[k\right]}$, is calculated by minimizing the squared loss function ${(R^{\left[k\right]}-X_j^{*}\lambda )}^T(R^{\left[k\right]}-X_j^{*}\lambda )$ with respect to $\lambda$. Therefore, the least squares estimate is ${\widehat{\lambda }}_j^{\left[k\right]}={\left[{\left(X_j^{*}\right)}^T\left(X_j^{*}\right)\right]}^{-1}{\left(X_j^{*}\right)}^T{R}^{\left[k\right]}$, the corresponding hat matrix is ${\mathcal{H}}_d=\left(X_j^{*}\right){\left[{\left(X_j^{*}\right)}^T\left(X_j^{*}\right)\right]}^{-1}{\left(X_j^{*}\right)}^T$ and the residual sum of squares is $RS{S}_j^{\left[k\right]}={(R^{\left[k\right]}-X_j^{*}{\widehat{\lambda }}_j^{\left[k\right]})}^T(R^{\left[k\right]}-X_j^{*}{\widehat{\lambda }}_j^{\left[k\right]})$. The chosen element ${\widehat{s}}_k$ is attained by:

$$\begin{array}{c}\hfill {\widehat{s}}_k={\mathrm{argmin}}_{1\le j\le p^{*}}\mathrm{gMDL}(Y^{*},RS{S}_j^{\left[k\right]},\mathrm{trace}\left({\mathcal{B}}_j^{\left[k\right]}\right)),\end{array}$$

(14)

where ${\mathcal{B}}_j^{\left[1\right]}={\mathcal{H}}_j$ and ${\mathcal{B}}_j^{\left[k\right]}=I-(I-{\mathcal{H}}_j)(I-\nu {\mathcal{H}}_{{\widehat{s}}_{k-1}}).\cdots .(I-\nu {\mathcal{H}}_{{\widehat{s}}_1})$ for $k>1$ is the first step boosting operator for choosing jth element in the kth iteration. Hence, there is an unique element $X_{{\widehat{s}}_k}^{*}$ to be selected at each iteration, and only the corresponding coefficient vector ${\theta }_{{\widehat{s}}_k}^{*\left[k\right]}$ changes, i.e., ${\theta }_{{\widehat{s}}_k}^{*\left[k\right]}={\theta }_{{\widehat{s}}_k}^{*[k-1]}+\nu {\widehat{\lambda }}_{{\widehat{s}}_k}^{\left[k\right]}$, where $\nu$ is the pre-specified step-size parameter. All the other ${\theta }_j^{*\left[k\right]}$ for $j\ne {\widehat{s}}_k$ keep unchanged. We repeat this procedure for K times and the number of iterations K can be estimated by

$$\widehat{K}={\mathrm{argmin}}_{1\le k\le K}\mathrm{gMDL}(Y^{*},RS{S}_{{\widehat{s}}_k}^{\left[k\right]},\mathrm{trace}\left({\mathcal{B}}^{\left[k\right]}\right)),$$

(15)

where ${\mathcal{B}}^{\left[k\right]}=I-(I-\nu {\mathcal{H}}_{{\widehat{s}}_k}).\cdots .(I-\nu {\mathcal{H}}_{{\widehat{s}}_1})$.

From the sparse boosting, we get the estimator of ${\theta }^{*}$ by ${\theta }^{*\left[\widehat{K}\right]}={({\left({\theta }_1^{*\left[\widehat{K}\right]}\right)}^T,\cdots ,{\left({\theta }_p^{*\left[\widehat{K}\right]}\right)}^T)}^T$.

For least squares loss function (8), $Y^{*}=Y-\widehat{\rho }D$ and $X^{*}=\pi$, and for least squares loss function (10), $Y^{*}=Y$ and $X^{*}=(WY,\pi )$.

Overall, as illustrated in Figure 1, the flowchart visually summarizes the methodologies and steps articulated in this paper, thereby enhancing the understanding of our proposed multi-step sparse boosting algorithm for high-dimensional additive spatial autoregressive models.

In Step I, the algorithm starts with parameter initialization, followed by iterations to identify optimal variables using the g-prior minimum description length (gMDL) criterion. Coefficients are updated until a stopping criterion is met, resulting in covariance matrix estimates.

Step II applies the sparse boosting algorithm again to refine variable selection and update coefficients. This approach enhances model accuracy and selection efficiency, making it well-suited for high-dimensional data analysis.

3. Simulation

In this section, we investigate the finite sample performance of the proposed methodology with Monte Carlo simulation studies. The data is generated from the following model:

$$Y=\rho WY+\sum _{j=1}^{1000}{g}_j\left(X_j\right)+ϵ,$$

(16)

where $X_1,\cdots ,X_{1000}$ are all i.i.d on $[-1,1]$. $g_1\left(x\right)=g_3\left(x\right)=2tan\left(x\right)$, $g_2\left(x\right)=g_4\left(x\right)=2sin\left(x\right)$ and $g_5\left(x\right)=\cdots g_{1000}\left(x\right)=0$. The error term $ϵ$ follows a normal distribution with mean 0 and variance ${\sigma }^2$. Similar to [41], the weight matrix is set to be $W=I_R\otimes E_m$, where $E_m=(l_m^T{l}_m-l_m)/(m-1)$, $l_m$ is the m-dimensional unit vector, and ⊗ is Kronecker product. Thus, the number of spatial units is $n=R\times m$. For comparison, we consider $(R,m)=(10,10)$ and $(20,20)$, corresponding to $n=100$ and 400 respectively. We evaluate three different values of $\rho =0.2,0.5,0.8$, representing weak to strong spatial dependence of the responses. Additionally, we consider ${\sigma }^2=0.5$ and 1. The degree of freedom for B-splines basis is set to be $\lfloor (N/\lfloor N/log(N)\rfloor \rfloor$, where $\lfloor x\rfloor$ is the integer part of x.

In face of ultra-high dimensional data, pre-screening can be adopted to reduce the dimensionality to a moderate size. In particular, we adopt marginal screening by $l^2$-norm [42] to screen out irrelevant covariates. More specifically, we regress Y on each covariate to construct the marginal model $Y={\rho }_{0}WY+g_{0j}\left(X_j\right)+{ϵ}_{0j},j=1,\cdots ,p$. The empirical $l^2$-norm of an estimated function ${\widehat{g}}_{0j}\left(X_j\right)={({\widehat{g}}_{0j}\left(X_{1j}\right),\cdots ,{\widehat{g}}_{0j}\left(X_{nj}\right))}^T$ is defined as:

$$\parallel {\widehat{g}}_{0j}\left(X_j\right){\parallel }^2={\displaystyle \frac{1}{n}}\left({\widehat{g}}_{0j}\left(X_{1j}\right),\cdots ,{\widehat{g}}_{0j}\left(X_{nj}\right)\right){\left({\widehat{g}}_{0j}\left(X_{1j}\right),\cdots ,{\widehat{g}}_{0j}\left(X_{nj}\right)\right)}^T.$$

(17)

A greater value of this measure suggest a stronger association between the covariate and the response. Following the recommendation by [42,43,44,45], the covariates with the largest $\lfloor {\displaystyle \frac{n}{log\left(n\right)}}\rfloor$$l^2$-norm are selected. Thus, for sample size 100 and 400, the first 21 and 66 variables in the ranked list are selected to conduct downstream analysis, respectively. Thereafter, we proceed to use our proposed multi-step sparse boosting approach to build up the final parsimonious model. For comparison, besides the proposed method using sparse boosting in each step, we also examine other approaches such as $l_2$-boosting, lasso regression and elastic net regression in each step.

In our implementation, we set the maximum number of boosting iterations to $K=500$ and the elastic net mixing parameter to $0.5$. Penalized methods were executed using the R package glmnet, with tuning parameters selected via 5-fold cross-validation. To assess the performance of our approach, we analyze results from 500 replications, reporting the following metrics:

• S: coverage probability that the top $\lfloor {\displaystyle \frac{n}{log\left(n\right)}}\rfloor$ covariates after screening includes all important covariates;
• TP: the median of true positives;
• FP: the median of false positives;
• Size: the median of model sizes;
• ISPE: the average of in-sample prediction errors defined as ${(Y-\widehat{Y})}^T((Y-\widehat{Y})/n$;
• RMISE: the average of root mean integrated squared errors defined as
• $\sqrt{{\displaystyle \frac{1}{n}}{\sum }_{j=1}^4{(g_j\left(X_j\right)-{\widehat{g}}_j\left(X_j\right))}^T(g_j\left(X_j\right)-{\widehat{g}}_j\left(X_j\right))}$;
• Bias($\rho$): the mean bias of $\rho$;
• Bias($\sigma$): the mean bias of $\sigma$.

The simulation results presented in Table 1 and Table 2, based on 500 replications, provide comprehensive insights into the performance of various methods for variable screening and selection under two conditions ($\sigma =0.5$ and $\sigma =1.0$).

At $\sigma =0.5$, M1 consistently achieves high coverage probabilities, reaching 0.99 when $\rho =0.8$, while maintaining a true positive count of 4 and significantly reducing the false positive rates compared to other methods. For example, M2 has a maximum of 12 false positives when $n=100$ and $\rho =0.2$, while M3 and M4 peak at 15 and 16 false positives, respectively. In terms of in-sample prediction error (ISPE), M1 demonstrates robust performance with values ranging from 1.123 to 3.127. M2, M3, and M4 show higher ISPE values, indicating a tendency toward overfitting, with M2 reaching 3.405 at $n=100$ and $\rho =0.8$, M3 reaching 3.443, and M4 reaching 3.64. M1 also excels in root mean integrated squared error (RMISE), with values between 0.534 and 0.754 at $n=100$ and $n=400$, while M2 ranges from 0.579 to 0.782. At $n=100$, M3 and M4 have even higher values, further demonstrating M1’s effectiveness in making accurate predictions.

At $\sigma =1.0$, similar trends are observed with the performance of M1 being even more pronounced. For instance, M1 maintains a coverage probability of 1.00 when $\rho =0.8$. The ISPE for M1 ranges from 2.280 to 5.771, while M4’s ISPE peaks at 6.681, highlighting M1’s advantage in prediction accuracy. Additionally, the false positive rates for M2, M3, and M4 increase, reaching more than 40 at $n=400$, further affirming M1’s superior performance. Despite maintaining a median true positive count of 4, M1 shows lower bias in parameter estimation, with Bias($\rho$) ranging from -0.107 to 0.002 compared to higher biases in the other methods.

In summary, the multi-step sparse boosting method (M1) outperforms alternative approaches in variable selection and parameter estimation, excelling in maintaining lower false positive rates and superior estimation accuracy. This establishes the multi-step sparse boosting method as a highly effective method for high-dimensional data analysis, particularly in scenarios involving complex relationships among variables.

4. Real Data Analysis

In this section, we apply the proposed method to the Boston housing price data, originally collected by Harrison and Rubinfield [46] and integrated by Gilley and Pace [47]. The dataset is available for download from http://ugrad.stat.ubc.ca/R/library/mlbench/html/BostonHousing.html. It comprises median house prices observed in 506 census tracts in the Boston area in 1970, alongside a set of variables presumed to influence house prices. Table 3 provides a detailed description of these variables:

The Boston housing price dataset is widely used in spatial econometrics and has been extensively studied in the literature. For the same data set, Xie et al. [48] adopted spatial autoregressive model, Li and Mei [13] considered partially linear spatial autoregressive model and Du et al. [16] used partially linear additive spatial autoregressive models. In the following, we will consider the additive spatial autoregressive model

$$\begin{array}{cc}\hfill {\mathrm{MEDV}}_i=& \rho \sum _{j=1}^{506}{w}_{ij}{\mathrm{MEDV}}_j+g_1\left(\mathrm{CRIM}\right)+g_2\left(\mathrm{ZN}\right)+g_3\left(\mathrm{B}\right)+g_4\left(\mathrm{RM}\right)+g_5\left(\mathrm{DIS}\right)+g_6\left(\mathrm{NOX}\right)+g_7\left(\mathrm{AGE}\right)\hfill \\ \hfill & +g_8\left(\mathrm{INDUS}\right)+g_9\left(\mathrm{RAD}\right)+g_{10}\left(\mathrm{PTRATIO}\right)+g_{11}\left(\mathrm{LSTAT}\right)+g_{12}\left(\mathrm{TAX}\right)+{ϵ}_i,i=1,\cdots ,506.\hfill \end{array}$$

(18)

where ${\mathrm{MEDV}}_i$ represents the median house value for census tract i, $\rho$ is the spatial autoregressive parameter, $w_{ij}$ denotes the spatial weight matrix with entries set as $\max (1-{\displaystyle \frac{d_{ij}}{d_0}},0)$, where $d_{ij}$ is the Euclidean distance based on longitude and latitude coordinates of any two houses, $d_0$ is the threshold distance (set to $0.05$ following Su and Yang [4]), resulting in a weight matrix with 19.1% nonzero elements. For comparative analysis, we consider the following methods explored in our simulation study: multi-step sparse boosting, multi-step $l_2$-boosting, multi-step lasso and multi-step elastic net. These methods will be used to estimate the functions $g_j(.)$ and $\rho$, aiming to identify the significant determinants of housing prices while accounting for spatial dependencies among census tracts.

The results of variable selection, estimation and prediction are summarized in Table 4. In addition to the number of selected variables, we also report the in-sample prediction error (ISPE) ${\sum }_{i=1}^n{(Y_i-{\widehat{Y}}_i)}^2$ and out-of-sample prediction error (OSPE) measured by 5-fold cross-validation.

From Table 4, we observe that, compared to traditional lasso and elastic net methods, the $l_2$ Boosting and Sparse $l_2$ Boosting methods for the varying-coefficient model exhibit smaller OSPE, while also producing relatively sparser models. The smaller ISPE may be due to overfitting in the traditional models which choosing all of the variables. These conventional models exhibit errors that increase by a factor of 6 to 7 during out-of-sample predictions, whereas our proposed model demonstrates greater stability with smaller fluctuations when comparing ISPE to OSPE. Our Sparse $L_2$ Boosting method for the varying-coefficient AFT model performs remarkably well in terms of sparsity, estimation, and prediction.

Additionally, we present a Venn diagram in Figure 2 showing the overlapping genes identified by all different methods. We observe that Sparse$L_2$Boosting and $l_2$Boosting select 2 variables (RM and LSTAT) which mostly affect the Boston housing prices in common, while lasso and elastic net methods identify all the 12 variables in common. However, Sparse$L_2$Boosting and $l_2$Boosting produce distinct lists of selected variables compared to the other methods, with only 2 variables selected by all methods.

We plot the estimated curves of varying coefficients for RM and LSTAT, which are the two selected variables by sparse boosting method, with their 95% confidence bands constructed by 500 bootstrap resamples in Figure 3. All of the functions are quite different from a straight line and this suggest that the varying-coefficient model is more appropriate to describe the covariate effects on boston housing prices in our data. Using more sophisticated semiparametric model specification may provide more accurate model for the high-dimensional analysis. Furthermore, we observe different functional forms for the two impact variables.

5. Concluding Remarks

This paper introduces a useful multi-step sparse boosting algorithm specifically designed for additive spatial autoregressive models with high-dimensional covariates, addressing critical aspects of variable selection and parameter estimation. Our methodology effectively facilitates model-based prediction without the burden of time-consuming tuning parameter selection, thereby streamlining the modeling process.

Using B-spline basis functions for approximating varying coefficients is a critical component of our approach, which allows for smooth and accurate representations even when the underlying assumptions regarding smoothness may be violated, enhancing the robustness of the model. Moreover, a noteworthy aspect of our finding is the two-step sparse boosting approach employing the generalized Minimum Description Length (gMDL) model selection criterion. The gMDL criterion utilizes a data-driven penalty for squared loss, promoting model parsimony while effectively filtering out irrelevant covariates.

The simulation results illustrate the efficacy of our proposed approach. Our multi-step sparse boosting method consistently exhibits fewer false positives than alternative methods, while true positive rates accurately reflect the actual number of relevant variables. In contrast, traditional methods utilizing boosting, lasso, or elastic net tend to suffer from overfitting, as evidenced by lower in-sample prediction errors (ISPE) and larger model sizes. Meanwhile, the real-world application of our methodology to the Boston housing price data further confirms its effectiveness and stability in identifying significant predictors of housing prices. By accurately selecting key variables, our approach enhances the model’s predictive performance and interpretability.

However, several limitations of our study must be acknowledged. Firstly, although B-spline basis functions offer significant advantages, they may not fully capture the intricate complexities of all data relationships, potentially limiting the model’s applicability to more diverse or nuanced datasets. Therefore, while our research encompasses thorough simulations and empirical testing on the housing price datasets, it is imperative to extend evaluations across a wider array of datasets.

This method can be extended to various fields, not limited to environmental studies, healthcare, and social sciences. Additionally, it is applicable to both nonparametric and semiparametric models, although our focus here has been on nonparametric approaches. This broader analysis is crucial to comprehensively validate the model’s versatility and effectiveness across various contexts within the disciplines of econometrics and statistics.