Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
86,527 characters · 14 sections · 53 citation commands
Gradient Boosting for Spatial Regression Models with Autoregressive Disturbances
\email{[email removed]}
\affil*[1]{\orgdiv{Bielefeld University}, \orgname{Center for Mathematical Economics}, \orgaddress{\street{Universitätsstraße 25}, \city{Bielefeld}, \postcode{33615}, \state{NW}, \country{Germany}}}
\abstract{Researchers in urban and regional studies increasingly deal with spatial data that reflects geographic location and spatial relationships. As a framework for dealing with the unique nature of spatial data, various spatial regression models have been introduced. In this article, a novel model-based gradient boosting algorithm for spatial regression models with autoregressive disturbances is proposed. Due to the modular nature, the approach provides an alternative estimation procedure which is feasible even in high-dimensional settings where established quasi-maximum likelihood or generalized method of moments estimators do not yield unique solutions. The approach additionally enables data-driven variable and model selection in low- as well as high-dimensional settings. Since the bias-variance trade-off is also controlled in the algorithm, implicit regularization is imposed which improves prediction accuracy on out-of-sample spatial data. Detailed simulation studies regarding the performance of estimation, prediction and variable selection in low- and high-dimensional settings confirm proper functionality of the proposed methodology. To illustrative the functionality of the model-based gradient boosting algorithm, a case study is presented where the life expectancy in German districts is modeled incorporating a potential spatial dependence structure.}
In various real-world applications, applied researchers have to engage increasingly with data that exhibits distinct locational attributes, which establishes a connection between objects in space. Such spatial data usually includes a cross-sectional variable in which each observation corresponds to a spatial unit like a geographical location. Thus, dealing with spatial data poses certain challenges since classical statistical regression models heavily rely on the independence assumption between observations which simplifies the model but is commonly violated in the presence of spatial units. Therefore, so-called spatial regression models have been developed which explicitly incorporate the spatial dependence structure of the data lesage2009. Particularly, spatial dependencies can arise due to spatial autocorrelation which describes similarities between geographical locations in space. For instance, housing prices in cities might be more similar in neighboring locations of close proximity. Thus, under the assumption that the spatial dependencies only arise through the error terms, cliff1973 proposed based on whittle1954 to model the spatial autocorrelation in the disturbance process. Then, for any given variable of interest, the error term for each location depends on a weighted average of the disturbances in connected locations. Thus, spatial spillovers between observations are quantified through the spatial structure of the disturbances anselin1988.
Considering that the complexity of the linear predictor can potentially grow quickly in size depending on the modeling choices of the practitioner, approaches to model choice and variable selection have become increasingly more important. Particularly, many candidate spatial regression models with autoregressive disturbances can be chosen by the practitioner by the inclusion of additional spatially lagged independent variables in the linear predictor. However, model choice is non-trivial since modeling spatial dependencies can induce a non-nested model structure. If model estimation is based on the quasi-maximum likelihood (QML), a quantifiable value for model comparison can be obtained based on a likelihood ratio test lee2004, liu2019. Instead of utilizing specifically designed criteria, regularization techniques are popular alternatives for model choice and variable selection fahrmeier2013. Particularly, the least absolute shrinkage and selection operator tibshirani1996 is a popular example of such a regularization technique which has also been recently extended to spatial regression model with autoregressive disturbances cai2019,cai2020. Another option in classical linear regression models is the so-called model-based gradient boosting algorithm fahrmeier2013, hepp2016. Although originally proposed in the domain of machine learning for classification problems, gradient boosting has been extended for statistical regression models where it is known as component-wise, model-based or statistical gradient boosting. In principle, the algorithm is iterative in nature where the estimation problem reduces to fitting base-learners to the negative gradient of a prespecified loss function related to the statistical model of interest. In each iteration of the algorithm only the best performing base-learner is chosen from which a small fraction is added to the current linear predictor. Stopping the algorithm early allows for data-driven variable selection and yields interpretable results at each iteration mayr2014. Due to the modular and iterative nature of model-based gradient boosting, it remains a feasible approach even in high-dimensional data settings where the number of variables exceeds the number of observations bühlmann2006, bühlmann2007. Therefore, it is not surprising that model-based gradient boosting has been extended to a variety of regression models of increased complexity (see, for example, jobst2024, balzer2025).
Although sparse boosting for semi-parametric additive spatial autoregressive models has been recently explored yue2025, no prior work appears to have addressed a potential extension of model-based gradient boosting for (parametric) spatial regression models with autoregressive disturbances. To this end, the model-based gradient boosting algorithm in the mboost package bühlmann2007, hothorn2010, hofner2014, hofner2015 for fitting generalized linear, additive and interaction models of potential high-dimensionality in the programming language R R is extended to accommodate spatial regression models with autoregressive disturbances. To investigate proper functionality of the proposed model-based gradient boosting algorithm, in-depth simulation studies in low- and high-dimensional linear settings are conducted. The focus lies primarily on the evaluation of estimation, variable selection and prediction. To illustrate the potential application of model-based gradient boosting for spatial regression models with autoregessive disturbances, an application concerned with modeling the life expectancy in German districts is presented which heavily draws on the "Indikatoren und Karten zur Raum- und Stadtentwicklung" (INKAR) data base which includes a rich variety of variables BBSR2024.
The structure in this article is as follows: In Section (ref), the mathematical, theoretical framework of spatial regression models with autoregressive disturbances and model-based gradient boosting is introduced where both concepts are combined thereafter. Afterward, outcomes for the simulation studies are discussed in Section (ref). A description of the context, data situation and variables for the case study as well as the results for the application of model-based gradient boosting are presented in Section (ref). The article finishes then with a conclusion and a discussion in Section (ref).
Let $n \in \mathbb{N}$ denote the number of observations in a spatial dataset. For each $i \in \{1,\dots,n\}$, consider following spatial regression model with autoregressive disturbances
where $\bm{y}$ is a $n \times 1$ vector of observations, $\bm{X}$ is the $n \times p$ design matrix of $p \in \mathbb{N}$ exogenous variables, $\bm{\beta}$ are corresponding $p \times 1$ coefficients and $\bm{u}$ is the $n \times 1$ vector of disturbances. The spatial dependence in the data is assumed to enter the model in two ways. First, the disturbances are modeled as an autoregressive process which depend on a spatial autoregressive parameter $\lambda \in \{-1,1\}$, a spatial weight matrix $\bm{W}$ of size $n \times n$ that captures spatial connections between observations and $n \times 1$ idiosyncratic random innovations $\bm{\epsilon}$. Second, spatial lags of exogenous variables $\bm{W}\bm{X}$ and corresponding $p \times 1$ coefficients $\bm{\theta}$ are also included in modeling the spatial dependent variable of interest. Thus, the model in Equation (ref) is a generalization of special cases of spatial regression models and is labeled the spatial Durbin error model (SDEM). In the SDEM, spatial autocorrelation is accounted for in both the explanatory variables, as well as the error term. However, removing the spatially lagged independent variables results in the simpler spatial error model (SEM). In contrast, retaining the spatially lagged independent variables but removing the autoregressive nature of the disturbances results in the simpler spatial cross-regressive model (SLX). These special cases are discussed in detail in Appendices (ref) and (ref) anselin1988, lesage2009, halleck2015.
Additionally, regularity conditions have to be imposed which are useful in ensuring a proper estimation procedure. Assumption (ref) imposes homoskedasticity in the innovations. However, normality is formally not required. In Assumption (ref), the number of observations are linked to the spatial weight matrix which is satisfied if the $\bm{W}$ is row-normalized that is assumed to hold throughout this article. The stability condition $|\lambda| < 1$ ensures invertibility of $\bm{I} - \lambda \bm{W}$ and thus uniqueness of the autoregressive disturbances in terms of innovations. Similarly, Assumption (ref) ensures that the degree of spatial autocorrelation remains manageable lee2004.
The model in Equation (ref) can be written more compactly by combining all exogenous variables and the corresponding spatial lags into one design matrix $\bm{Z} =[\bm{X}, \bm{W} \bm{X}]$ and by stacking the corresponding coefficient vertically $\bm{\delta} = (\bm{\bm{\beta}, \bm{\theta}})^{\prime}$ as
where $\bm{\eta} = \bm{Z} \bm{\delta}$ is the so-called linear predictor. Although not the focus in this article, it is worth to note that the general model formulation also allows for the inclusion of additional spatial lags of exogenous variables, for example, $\bm{W}^2\bm{X}$ or even $\bm{W}^3\bm{X}$. To estimate the coefficients of the exogenous variables and the corresponding spatial lags of exogenous variables, it is convenient to transform the model into a single equation form as
Then Assumption (ref) has to be additionally imposed for a proper estimation procedure which states the usual full column rank condition of $\bm{Z}$.
If the autoregressive parameter $\lambda$ is known, then the loss function of Equation (ref) is the squared Mahalanobis distance of the residual vector
where the error covariance structure is induced by spatial dependence, and
The negative gradient vector is obtained as the derivative of the loss function with respect to the linear predictor as
which yields all necessary ingredients for the model-based gradient boosting algorithm kelejian1999, cai2019, cai2020.
Given the relevant expressions for the ingredients, model-based gradient boosting for the SDEM can be implemented. In principle, the ingredients play an important role in the iterations of the algorithm. Generally, the standard interpretation of boosting as the steepest descent in function space means that the algorithm reduces the empirical risk via the base-learners in an iterative fashion friedman2001. The base-learners refer to the functional form of the exogenous input variables. The algorithm begins with an empty model, and fits the specified base-learners to the negative gradient of the chosen loss function. Thus, proper functionality of the algorithm requires a pre-specified loss function which can be quite general. Afterward, the residual sum of squares is computed for each base-learner separately and the linear predictor is updated by a small fraction of the best performing base-learner. The algorithm then reevaluates the negative gradient and updates the linear predictor in an iterative manner until the specified number of boosting iterations are reached friedman2001, bühlmann2007, mayr2014. Algorithm (ref) adapts model-based gradient boosting for the SDEM and does not impose any limitations on the number of potential independent variables $q$. Indeed, the great advantage of model-based gradient boosting is the feasibility in high-dimensional settings where the number of variables is larger than the number of available observations. If the squared error is utilized as the loss function, model-based gradient boosting yields a consistent estimator in both low- and high-dimensional settings zhang2005, bühlmann2006.
Additionally, bühlmann2007 formally show that model-based gradient boosting with squared error loss converges to the ordinary least squares (OLS) solution if the number of boosting iterations $m_{\text{stop}}$ is chosen sufficiently large. As the main tuning parameter in the algorithm, $m_{\text{stop}}$ controls the so-called bias-variance trade-off. Thus, the accuracy of prediction can be improved and overfitting behavior mitigated. Due to the modular nature, the algorithm yields an interpretable solution at each iteration such that sparser models can be obtained by stopping the algorithm early instead of convergence. Since the algorithm only updates the linear predictor by one component at each iteration, variable selection and effects estimation shrinkage is also accounted for mayr2012.
The feasibility of Algorithm (ref) strongly relies on the assumption that the spatial autoregressive parameter $\lambda$ is apriori known. However, in real-world application settings, $\lambda$ is unknown which implies that $\bm{\Omega}(\lambda)$ occurring in the squared Mahalanobis distance and the negative gradient cannot be evaluated. Therefore, a three-step model-based gradient boosting procedure is proposed to enable the feasibility of Algorithm (ref) which relies on replacing the unknown $\lambda$ and $\bm{\Omega}(\lambda)$ by the estimated counterparts $\hat{\lambda}$ and $\bm{\Omega}(\hat{\lambda})$. In the first step, the model in Equation (ref) is written as
temporarily ignoring the potential autoregressive structure of the disturbances. The model in Equation (ref) can then be estimated using a variety of methods, as long as the resulting estimator $\bm{\tilde{\delta}}$ is consistent. For low-dimensional settings, a natural choice for the estimator is OLS. However, since the loss function associated with Equation (ref) is simply the squared error loss, model-based gradient boosting can also be employed. In high-dimensional settings, OLS does not yield unique solutions, so model-based gradient boosting which provides a consistent estimator is utilized instead.
In the second step, let $\bm{\tilde{u}} = \bm{y} -\bm{Z}\bm{\tilde{\delta}}$ denote the predictors of $\bm{u}$ based on a consistent estimator $\bm{\tilde{\delta}}$. Define $\bm{\bar{u}} = \bm{W}\bm{u}$, $\bm{\bar{\bar{u}}} = \bm{W}\bm{W}\bm{u}$ and the corresponding expressions based on the predictors as $\bm{\tilde{\bar{u}}} = \bm{W}\bm{\tilde{u}}$ and $\bm{\tilde{\bar{\bar{u}}}} = \bm{W}\bm{W}\bm{\tilde{u}}$. Adapt an identical notation pattern to $\bm{\bar{\epsilon}} = \bm{W} \bm{\epsilon}$. Then, if Assumptions (ref) to (ref) hold, following three moments can be obtained
where $\text{tr}(\cdot)$ denotes the trace of any matrix. Since the innovations can be written in terms of $\bm{\bar{u}}$ and $\bm{\bar{\bar{u}}}$ as $\bm{\epsilon} = \bm{u} - \lambda\bm{\bar{u}}$ and $\bm{\bar{\epsilon}} = \bm{\bar{u}} - \lambda\bm{\bar{\bar{u}}}$, a system of three equations can be obtained based on Equations (ref) and (ref)
The expressions for $\bm{\Gamma}(\cdot)$ and $\bm{\gamma}$ are given as
Replacing the moments in Equation (ref) by the corresponding sample moments yields
where
and $\bm{\nu}(\lambda, \sigma^2)$ is interpreted as a $3 \times 1$ residual vector. The estimators for $\lambda$ and $\sigma^2$ are obtained using nonlinear least squares and are denoted by $\hat{\lambda}$ and $\hat{\sigma}^2$. Based on Equation (ref), the non-linear least squares estimators are defined as
Let the Assumptions (ref) to (ref) hold. Then the nonlinear least squares estimators $\hat{\lambda}$ and $\hat{\sigma}^2$ are consistent estimators of $\lambda$ and $\sigma^2$ in the sense that $\hat{\lambda} \to_p \lambda$ and $\hat{\sigma}^2 \to_p \sigma^2$ for $n \to \infty$ sufficiently large.
In the third step, replace $\lambda$ in $\bm{\Omega}(\lambda)$ by the estimated counterparts $\hat{\lambda}$ yielding $\bm{\Omega}(\hat{\lambda})$. Thus, the squared Mahalanobis distance and negative gradient become
Finally replace the expressions in Algorithm (ref) by the corresponding expressions based on the estimated counterparts in Equations (ref) and (ref) which yields a feasible model-based gradient boosting algorithm for the SDEM kelejian1999.
In model-based gradient boosting, the standard approach for model and variable selection is through early stopping via the stopping criterion $m_{\text{stop}}$. Generally, $m_{\text{stop}}$ is usually chosen by means of k-fold cross-validation, subsampling or bootstrapping which have the tendency to include to many (non-informative) variables which a induces non-parsimonious final model mayr2012. Since the quality of the generalized moment estimators $\hat{\lambda}$ and $\hat{\sigma}^2$ depend on the quality of the predictors of $\bm{u}$, the inclusion of non-informative variables severely impacts the final estimators.
To mitigate the consequences and obtain sparser final models, the deselection algorithm proposed in stromer2022 is adapted to model-based gradient boosting for the SDEM which can be easily utilized in the first and third step of the procedure. In principle, the idea is to perform model-based gradient boosting and determine the optimal stopping iteration $m_{\text{opt}}$ using standard cross-validation techniques. Subsequently, variables that contribute the least to risk reduction are deselected, and model-based gradient boosting is reapplied using the remaining variables and the previously determined $m_{\text{opt}}$. More formally, let $\mathds{1}(\cdot)$ denote the indicator function, then the attributable risk reduction given by
is assumed to be a measure of importance for the $j$-th base-learner. Then $j^{{*}^{[m]}}$ is the component corresponding to the selected base-learner, $\left(r^{[m-1]} - r^{[m]}\right)$ is the risk reduction in iteration $m$ where $r^{[m-1]}$ and $r^{[m]}$ are corresponding risks. A base-learner and thereby the corresponding variable is deselected if
where $\tau \in (0,1)$ is a pre-specified threshold and $\left(r^{[0]} - r^{[m_{\text{stop}}]}\right)$ the total risk reduction. Thus, variables remain only in the model if the relative risk contribution is equal or larger than the threshold $\tau$. The choice of $\tau$ usually depends on the particular research situation at hand. However, a typical recommendation is to keep values low, that is, for example, $\tau = 0.01$ stromer2022. The complete procedure incorporating the deselection approach into model-based gradient boosting for the SDEM is given in Algorithm (ref).
To evaluate the performance of the proposed three-step model-based gradient boosting algorithm, simulation studies are conducted with the ingredients derived in Section (ref). Particularly, the performance of estimation, variable selection and prediction in low- as well as high-dimensional linear settings are evaluated. Additionally, an evaluation of the performance of the deselection algorithm is provided. The study design is motivated by the case study of modeling life expectancy in German districts in Section (ref). Specifically, the number of observations is fixed at $n = 400$. In contrast, the number of independent variables is varied between $q = 20$ and $q = 800$, indicating a low- $(n > q)$ and high-dimensional $(n < q)$ linear setting. The true data generating process is given by
where the variables are independently and identically drawn from the uniform distribution $\bm{X} \sim U(-2,2)$. The spatial autoregressive parameter is varied throughout the simulation study by $\lambda \in \{-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8\}$ and the innovations are normally distributed according to $\bm{\epsilon} \sim N(0,\sigma^2)$ with $\sigma^2 =1$. The spatial weight matrix $\bm{W}$ is generated based on a circular world in which each location is directly related to the five locations before and after, that is, $K = 5$. Additionally, $\bm{W}$ is row-normalized such that each row sums up to one. Simulation studies for varying number of related locations $K \in \{1,2,3,5,10,20\}$ is given in Appendix (ref). In the model-based gradient boosting algorithm, the corresponding base-learners are specified as simple linear regression models due to the nature of the data generating process. The learning rate is set to $s = 0.1$ since that is the usual practice (see, for example, schmid2008, mayr2012, hofner2014). The optimal stopping criterion $m_{\text{opt}}$ is found by minimizing the empirical risk via 25-fold subsampling. In each simulation setting, a total of $n_{\text{sim}} = 100$ repetitions are conducted. Additionally, different approaches for the consistent estimator $\bm{\tilde{\delta}}$ in the first step are considered, along with their impact on the final results. Specifically, the reported methods are first-step OLS (LS-GB), first-step gradient boosting (GB-GB), and first-step gradient boosting with deselection (DS-GB) where applicable.
Regarding the performance of variable selection and deselection, the criteria are chosen based on the confusion matrix. In particular, the reported variable selection criteria are the true positive rate (TPR), which is the proportion of correctly selected variables out of all true informative variables, the true negative rate (TNR), which is the proportion of correctly non-selected variables out of all true non-informative variables and the false discovery rate (FDR), which is the proportion of non-informative variables in the set of all selected variables stehmann1997.
The performance of estimation is evaluated by reporting the bias, the mean squared error (MSE) and the empirical standard error (ESE) for $\lambda$ defined as
For all proposed performance criteria, lower values are always preferred. Additionally, the effects shrinkage estimation for the independent variables is evaluated via visualization by boxplots to highlight the median, outliers and quartiles over 100 repetitions morris2019.
Furthermore, the prediction accuracy is evaluated based on an additional test data set. The test data $\bm{y_{\text{test}}}$ is generated according to the same data generating process as the train data with an identical number of observations $n_{\text{test}} = 400$. The chosen criteria are the quasi negative log-likelihood (NLL), the root mean squared error of prediction (RMSEP) and mean absolute error of prediction (MAEP) defined as
For all proposed performance criteria of prediction accuracy, lower values are always preferred.
The simulation study is conducted in the programming language R R. The QML and generalized method of moments (GMM) estimation of the SDEM in the low-dimensional linear setting is performed via the spatialreg package bivand2021, pebesma2023. The presented graphics are created with the tidyverse packages tidyvere2019. Model-based gradient boosting for generalized, additive and interaction models can be found in the mboost package bühlmann2007, hothorn2010, hofner2014, hofner2015. An implementation for model-based gradient boosting for the SDEM via the novel spatial error family incorporating the deselection algorithm and the R code for reproducibility of all simulation studies can be found in the GitHub repository \url{https://github.com/micbalz/SpatRegBoost}.
In Table (ref), the average selection rates over 100 repetitions for the low-dimensional linear settings can be seen. The results show a consistent TPR of $100\%$ across all spatial autoregressive parameters $\lambda$, indicating that all informative variables are selected on average. In contrast, for negative values of $\lambda$, a TNR of around $25\%$ is achieved, indicating that the method correctly avoids selecting only about four out of 16 non-informative variables on average, and incorrectly selects approximately 12 of them. As $\lambda$ becomes more positive, the TNR steadily increases, reaching its highest value of around $75\%$ for $\lambda = 0.8$, meaning that only about four out of 16 non-informative variables are incorrectly selected on average. However, a low TNR is accompanied by a high FDR across all values of $\lambda$. Although the TNR steadily increases as $\lambda$ becomes more positive, the corresponding decrease in FDR is less pronounced. Even at $\lambda = 0.8$, where the TNR is highest, the FDR remains substantial at approximately $45\%$. This indicates that nearly half of the selected variables are on average non-informative.
The performance of estimation for the linear effects in low-dimensional linear setting can be seen in Figure (ref). The estimated coefficients are presented by boxplots where the red horizontal line represents the true coefficient values according to the data generating process. Across all spatial autoregressive parameters $\lambda$, DS-GB manages to recover all informative variables. However, the coefficient estimates are slightly biased across all values of $\lambda$. In general, the magnitude of bias is stronger for spatial lags of the exogenous variables which is additionally increasing as $\lambda$ increases. However, a bias in the estimates for DS-GB is unsurprising and actually expected since the presented estimates for the linear effects are based on regularization via early stopping. In contrast, non-informative variables where rarely chosen with coefficient estimates clustered around $0$.
Additionally, the impact of the choice of the estimation method in the first step on the estimates of the $\lambda$ parameter are evaluated. The results are shown in Table (ref). Across all different spatial autoregressive parameters $\lambda$, the classical estimation strategies, namely QML and GMM struggle to reliable estimate $\lambda$ indicated by substantial downward biases. However, the bias seems to decrease as $\lambda$ increases and becomes positive yielding best values for $\lambda = 0.8$. Furthermore, QML outperforms corresponding GMM for positive $\lambda$ values. The results for LS-GB are identical to GMM, since the first step is based on OLS in both estimation strategies. The results change substantially for GB-GB and DS-GB. In principle, utilizing model-based gradient boosting in the first step to obtain consistent estimates of $\bm{\delta}$ leads to a decrease in the bias, MSE and ESE. The decrease in bias, MSE and ESE is further improved if model-based gradient boosting is combined with deselection in the first step. Moreover, the results show that GB-GB and DS-GB consistently outperform corresponding QML and GMM estimators.
Finally, the prediction performance on an independent test data set for the low-dimensionsal linear setting is evaluated. The results can be seen in Table (ref). In principle, the RMSEP, MAEP, and NLL are consistently higher across all values of the spatial autoregressive parameter $\lambda$ for QML and GMM, compared to the model-based gradient boosting algorithms. In contrast, the predictive performance of the model-based gradient boosting algorithms remains stable, with only minor differences at the decimal level.
In general, the presented simulation results show proper functionality of the model-based gradient boosting algorithm from the perspective of specificity and sensitivity in the low-dimensional linear setting across all considered spatial autoregressive parameters $\lambda$. This conclusion can be drawn based on the average selection rates where true informative variables are consistently selected across all 100 repetitions on average. Furthermore, the estimated linear effects show correct direction and algebraic sign. Nevertheless, biases are introduced due to the inherent regularization of model-based gradient boosting via early stopping. Although the average selection rates show high TNR and FDR, the actual impact in terms of estimated coefficients for non-informative variables remains very low, meaning the included non-informative variables are entering the final model only with low values. To provide intuition regarding the consistency and convergence of the proposed model-based gradient boosting algorithm, results are presented through a comparison with GMM coefficient estimates. These can be recovered by utilizing a sufficiently large stopping criterion, as demonstrated in Appendix (ref). Additionally, model-based gradient boosting outperforms the established QML and GMM estimators in estimating $\lambda$ when substantial noise is introduced through non-informative variables. As a result, the QML and GMM estimate exhibits a downward bias in the presence of such noise. This outcome is not unexpected, as spatial autocorrelation typically decreases when informative variables are added to the model. However, the simulation study also indicates that introducing non-informative noise can falsely reduce spatial autocorrelation, which has serious implications for model interpretation as it may lead to highly misleading conclusions. Utilizing model-based gradient boosting can effectively mitigate the consequences by not including non-informative variables, thereby improving quality of residuals in the first step and decreasing the bias across all values of $\lambda$. While model-based gradient boosting consistently outperforms classical QML and GMM estimation strategies in terms of predictive accuracy on unobserved test data, a notable drawback is its relatively FDR, indicating that the final model may still include many non-informative variables. However, its ability to exclude such variables improves drastically as $\lambda$ increases. Despite this improvement, the persistently high FDR suggests that model complexity remains an issue, even under stronger spatial dependence.
Similar to the low-dimensional linear setting, results are also evaluated for the high-dimensional linear setting. Here, only the GB-GB and DS-GB approaches are reported, as model-based gradient boosting remains the only feasible method when the number of variables exceeds the number of observations. Table (ref) presents the average selection rates. Across all values of the spatial autoregressive parameter $\lambda$, the TPR remains consistently at $100\%$, indicating that all informative variables are selected on average. In contrast to the low-dimensional linear setting, the TNR is also consistently high, demonstrating that most non-informative variables are successfully excluded. Notably, TNR improves with increasing $\lambda$, reaching approximately $94\%$ at $\lambda = 0.8$. Despite these promising findings, the FDR remains high suggesting that the final model still includes many non-informative variables. However, this outcome is not unexpected, given that only four informative variables are present and can easily be overshadowed by the large number of non-informative variables.
The estimation performance of the coefficients in the high-dimensional linear setting is illustrated in Figure (ref). Compared to the low-dimensional linear setting, the coefficient estimates exhibit a more pronounced shrinkage effect, particularly for the coefficients of the spatial lags of the exogenous variables. Furthermore, the shrinkage effect additionally intensifies as the spatial autoregressive parameter $\lambda$ increases and becomes more positive. Nevertheless, the algebraic signs and general direction of the coefficients remain consistently correct even in the presence of many non-informative variables.
Furthermore, the estimation performance for the spatial autoregressive parameter $\lambda$ in the high-dimensional linear setting can be seen in Table (ref). In principle, QML, GMM and LS-GB are not feasible anymore. Thus only GB-GB and DS-GB remain as feasible estimation strategies. The results for GB-GB indicate an increasing bias on the estimated $\hat{\lambda}$. This behavior can be explained by the increased shrinkage effect on the informative variables in combination with many non-informative variables since the generalized method of moment estimator heavily relies on the quality of the predictions of the residuals. In contrast, DS-GB almost always deselect all non-informative variables. Thus, even though the shrinkage effect is present, the additional noise from non-informative variables is removed which ensures a proper recovery of the true $\lambda$.
Finally, the predictive performance on independent test data in the high-dimensional linear setting is shown in Table (ref). As before, results are reported only for GB-GB and DS-GB, since QML, GMM, and LS-GB are infeasible in this setting. Although the performance itself is more varied, the absolute values of the performance criteria are comparable to those observed in the low-dimensional linear setting. As expected, the predictive accuracy of DS-GB improves with increasing values of $\lambda$ which corresponds to the reduction in bias observed for estimate $\hat{\lambda}$. However, a direct performance comparison with QML or GMM is not possible due to the infeasibility in the high-dimensional linear setting.
The results for the high-dimensional linear setting support the results of the low-dimensional linear settings and indicate proper functionality of the model-based gradient boosting algorithm from the perspective of specificity and sensitivity across all spatial autoregressive parameters $\lambda$. The TPR remains consistently at $100\%$ across all 100 repetitions on average. Although the estimated coefficients are stronger affected from shrinkage, the algebraic sign and general direction remain correct. The GB-GB exhibits a strong bias as $\lambda$ increases indicating that deselection in the first step is necessary to reliably estimate $\lambda$ in high-dimensional linear settings. Nevertheless, the great advantage is the feasibility of model-based gradient boosting in scenarios where the number of variables exceed the number of observations. In such situations, the established estimation strategies like QML and GMM fail entirely evidenced in the simulation study by missing values. For instance, the predictive performance as measured by the evaluation criteria, remains comparable in absolute terms to that of the low-dimensional linear setting. However, the results also reveal a consistently high FDR in the variable selection evaluation. Notably, no distinction is made between variables that are frequently selected and those that are selected only once, meaning all variables are equally weighted in the FDR computation. To address this limitation, the post-hoc deselection algorithm as proposed by stromer2022 is utilized to improve the FDR for model-based gradient boosting algorithm by removing variables that are likely selected by the algorithm merely by chance.
As evidenced in the results for the low- as well as high-dimensional linear setting, model-based gradient boosting for SDEM suffers from a "greedy" selection behavior which attributes to the inclusion of too many non-informative variables, thereby decreasing TNR and increasing FDR. The reason for the "greedy" nature lies in the fact that the solutions of model-based gradient boosting are optimal with respect to the $L_1$-arc-length. In fact, by taking any convex loss function, model-based gradient boosting is able to approximate the solution path of a strictly monotone $L_1$-regularized regression model. Therefore, model-based gradient boosting is unable to automatically deselect a variable once it has been added to the model hastie2007, hepp2016. To mitigate the consequences, the deselection algorithm for generalized linear, additive and interaction models is adapted to model-based gradient boosting and evaluated in the following stromer2022. The average selection rates after utilizing model-based gradient boosting with first step gradient boosting with deselection and an additional deselection step (DS-DS) are presented in Table (ref).
Indeed, the results clearly show that the TNR consistently remains at $100\%$ while the FDR decreases to $0\%$ on average. Thus, DS-DS successfully avoids selecting non-informative variables over all 100 repetitions, ensuring that only informative variables are included in the final model. This behavior is observed in low- as well as high-dimensional linear settings. Therefore, the deselection algorithm is able to mitigate the consequences of the inclusion of additional non-informative variables in boosted SDEM in both low- and high-dimensional linear settings. In fact, utilizing model-based gradient boosting in combination with deselection is strongly encouraged to ensure proper model and variable selection properties.
Life expectancy, defined as the average number of years a newborn in a given population is expected to live under current mortality conditions, is a statistical measure widely used to assess a variety of health, social, and economic outcomes marmot2005, cutler2006. Over the past several decades, the life expectancy has been steadily increasing across Europe indicating progress in healthcare, socioeconomic development, and public health interventions. Notably, countries in the European Union report an increase in longevity due to improved living standards, better medical care, and lower mortality from major diseases such as cardiovascular conditions and cancer. Nevertheless, significant regional disparities remain between and within countries. For instance, Germany, one of the largest and economically wealthiest nation in the European Union, has seen signs of stagnation in life expectancy in comparison to especially Northern European countries. Thus, the stagnation gives rise to important concerns about regional inequalities, healthcare system efficiency, lifestyle-related risks, and demographic changes oecd2024. In this case study, the goal is to model the life expectancy in German districts by utilizing spatial regression models with autoregressive disturbances to answer questions regarding underlying socio-economic and other determinants of life expectancy while incorporating district-level geographical disparities. In principle, the case study builds on recent previous research such as lampert2019, rau2020, siegel2022, jasilionis2023, marinetti2023, hoebel2025 which investigates current trends, socio-economic drivers of life expectancy, impact of the COVID-19 pandemic and patterns of mortality from both district-level and complete country perspectives.
To this end, a large-scale real-world data set, namely INKAR, is utilized. Managed by the Bundesinstitut für Bau-, Stadt- und Raumforschung, INKAR is an interactive (online) atlas about the living situation in Germany. It contains over 600 unique socio-demographic, socio-economic, and environmental indicators for distinct geographical locations allowing for evaluations of urban and rural disparities. In principle, INKAR is a panel data set relying on information starting from the year 1995. However, in this case study, the focus is on the cross-sectional life expectancy in German districts in the year 2019. Therefore, INKAR in the version of 2021 is utilized BBSR2024. The INKAR data set in the current as well as the 2021 version with the codebook for all indicators is freely available at \url{https://www.inkar.de/}.
To construct the data set suitable for this case study, the following preprocessing steps are necessary. First, the original number of districts available in INKAR version is 401. However, the data is merged with the German map data from the year 2024 available at \url{https://www.bkg.bund.de/}. The map data in the most recent version includes only 400 districts because the district "Eisenach" has been merged with a neighboring district in 2022 BKG2025. To ensure comparability between both data sets, "Eisenach" is thus removed from the INKAR 2021 data resulting in 400 districts of interest. Second, a spatial weight matrix has to be generated based on the specific spatial configuration of the German districts. Therefore, the centroids of the spatial polygons based on the shape file layer from the German map data are computed. To keep a close connection to the simulation study in Section (ref), a k-nearest neighbor structure is created utilizing $K = 10$, indicating that each district is assigned its ten geographically closest neighbors. Afterward, the resulting spatial weight matrix is row-normalized. Third, since the focus is primary on the variable selection evaluation, the indicators are allowed to be quite general. Thus, any domain knowledge is not imposed in the selection of variables. However, dubious indicators like the indicator number of a district as well as indicators describing similar phenomena are removed before model estimation. Furthermore, transformations for appropriate variables, that is, centering and scaling, are applied to ensure comparability across variables in terms of the scale.
The final data set is composed of 400 observations where each location corresponds to a district in Germany. These districts are divided again into 294 rural districts and 106 urban cities. A map of the life expectancy in German districts in the year 2019 can be seen in Figure (ref).
The map reveals distinct disparities in life expectancy across German districts, with clear spatial clustering among neighboring regions. Notably, the highest life expectancy is observed in southern Bavaria, near the borders with Austria and Switzerland. Additional areas of high life expectancy can be found in Lower Saxony, close to the Dutch border. In contrast, life expectancy in the new federal states (former East Germany) is generally lower than in the old federal states (former West Germany), indicating a clear east-west divide. Details on the data types and descriptions of the dependent and independent variables of interest are provided in Table (ref).
As evidenced by Figure (ref), ignoring spatial dependence in the estimation process may lead to a misrepresentation of the underlying determinants influencing life expectancy. Therefore, the objective is to model life expectancy in German districts for the year 2019 using spatial regression models with autoregressive disturbances. In particular, the SDEM is employed, as spatial dependence may arise not only in the disturbances but also through spatial lags of the exogenous variables. For example, the number of medical doctors in neighboring districts may influence life expectancy, especially given the relatively short distances between locations. Generally, the results are presented with a clear focus on model-based gradient boosting. Since the application settings is a low-dimensional linear setting, results are reported for the QML, GMM, LS-GB, GB-GB, DS-GB and DS-DS defined exactly as in the simulation study in Section (ref). All exogenous variables presented in Table (ref) as well as the corresponding first-order spatial lags and an intercept are included in the estimation yielding $65$ coefficients to estimate in a non-parsimonious model. Regarding the setup for the model-based gradient boosting algorithm, the narrative from the simulation study is followed and the learning rate is thus set to $s = 0.1$. Furthermore, for the optimization of the stopping criterion $m_{\text{stop}}$, the search is conducted by minimizing the empirical risk via 25-fold subsampling. For the deselection algorithm, $\tau$ is set to $0.01$. The estimated coefficients for all estimation strategies can be seen in Table (ref).
For brevity, the coefficients of the spatial lags of exogenous variables are not reported individually. Instead, the number of spatial lag variables (No. of $\bm{W}$) included in each model is presented. Regarding variable selection, the LS-GB approach reduces the total number of coefficients from $65$ to $53$ out of which $23$ can be attributed to the spatial lags of exogenous variables. A similar but improved outcome is observed for GB-GB where the number of variables is reduced to $41$ out of which $15$ correspond to the coefficients of the spatial lags of exogenous variables. The best results in terms of variables selection are achieved when model-based gradient boosting is combined with deselection. For DS-GB, the number of selected variables is reduced to $18$ with no coefficients of the spatial lags of exogenous variables included. DS-DS further improves performance of variable selection through post-hoc deselection, resulting in a final model with only $8$ variables. These results suggest that many of the variables listed in Table (ref) are non-informative. Notably, model-based gradient boosting combined with deselection effectively reduces the initially non-parsimonious SDEM to a simpler and more parsimonious SEM in a data-driven manner, without the need for formal likelihood-based model selection criteria.
Regarding the spatial autoregressive parameter, all estimation strategies consistently indicate positive spatial dependence between neighboring districts, although the strength varies. Consistent with the downward bias observed in the simulation study, both QML and GMM estimators tend to underestimate the spatial autoregressive parameter $\lambda$. This underestimation can be attributed to the inclusion of non-informative variables which are largely eliminated through model-based gradient boosting. Consequently, estimates of $\lambda$ obtained from model-based gradient boosting are generally higher, depending on the first step estimation method. Moreover, the algebraic sign and overall direction of the coefficients remain largely consistent across all estimation strategies. The primary differences lie in the magnitude of the coefficients, which varies due to the regularization effects introduced by early stopping in model-based gradient boosting. Based on the results for DS-DS, the most important variables explaining the life expectancy in German districts in 2019 are the average age of inhabitants, the unemployment rate, the debt quota, the proportion of employees with an academic degree, the household income per inhabitant, the rent prices and the number of care-dependent individuals. Since the independent variables are transformed by simple scaling and centering, the coefficients have an intuitive and simple interpretation. For instance, an increase in rent prices by one euro ceteris paribus increases the average life expectancy by $0.2109$ years on average. Conversely, holding all other variables constant, a $1$ percentage point increase in the share of private debtors is associated with an average $0.3923$ year decrease in average life expectancy. Additionally, $\sigma$ increases when model-based gradient boosting is utilized which is unsurprising, as fewer variables are retained in the final model, potentially leaving more unexplained variation in the innovations.
To summarize the findings, the positive spatial autoregressive parameter $\lambda$ indicates strong spatial dependence between districts, that is, the life expectancy in one district is strongly influenced by the life expectancy in neighboring districts. Thus, health outcomes cluster spatially as suggested in Figure (ref). Furthermore, an older average population is associated with lower life expectancy. Similarly, the number of care-dependent individuals negatively influence average life expectancy meaning that underlying poor health status decreases the average life expectancy. Districts with higher unemployment and debt quotas have on average lower average life expectancies indicating that financial hardship, fiscal stress or economic precarity is negatively correlated with longevity. In contrast, higher rent prices, higher household income and higher share of employees with an academic degree lead to a higher average life expectancy which shows that inhabitants in wealthier districts and better socioeconomic status live longer on average.
The key findings and main contributions of this article are: (a) Model-based gradient boosting is extended for spatial regression model with autoregressive disturbances, namely for the SDEM which unifies simpler models by considering spatial dependence in the linear predictor as well as the disturbances. (b) The algorithm is implemented in the mboost package bühlmann2007, hothorn2010, hofner2014, hofner2015 by providing a novel spatial error family. (c) Model-based gradient boosting for SDEM heavily relies on the knowledge about the spatial autoregressive parameter $\lambda$. However, in real-world application settings $\lambda$ is generally unknown making model-based gradient boosting infeasible. Therefore, a so-called feasible model-based gradient boosting algorithm for the SDEM is proposed which relies on replacing the unknown $\lambda$ in the ingredients by the corresponding estimated counterpart $\hat{\lambda}$ based on a generalized moment estimator proposed in kelejian1999. (d) In extensive simulation studies for low- as well as high-dimensional linear settings, the results show proper functionality of the proposed feasible model-based gradient boosting algorithm. Particularly, estimation is accompanied by high TPR values and coefficients are estimated with high accuracy although biases are introduced due to regularization via early stopping. Model-based gradient boosting also outperforms standard QML and GMM estimators in terms of bias, MSE and ESE when estimating the autoregressive parameter $\lambda$ in the presence of non-informative variables. Additionally, the predictive performance is always better for model-based gradient boosting in comparison to QML and GMM. (e) The great advantage of model-based gradient boosting is the feasibility in high-dimensional linear settings where the number of variables exceeds the number of observations. Due to the modular nature, model-based gradient boosting does not require refitting and allows for a direct interpretation of the effects of coefficients on the dependent variable of interest enabling a potential application for a wide range of real-world settings. (f) Additionally, the feasible model-based gradient boosting is applied in a real-world application setting where the life expectancy in German districts is modeled. Additional case studies revisiting classical spatial econometric data sets, namely Boston housing prices harrison1978 and Columbus crime rate anselin1988 are provided in Appendix (ref).
Naturally, limitations, improvements and extensions beyond the scope of this article have to be acknowledged. Although simulation studies have been provided, the scope of our simulation settings is of course by no means exhaustive. Further simulation studies for varying spatial weight matrices aim at partially address these limitations in Appendix (ref) but much more complex scenarios can be easily considered. For instance, additional correlated noise variables or higher spatial lags of exogenous variables may be included in the simulation settings. A major limitation of the proposed model-based gradient boosting algorithm is the restricted applicability. Specifically, the algorithm can only be applied in settings where the dependent variable does not appear as a spatial lag in the spatial regression model. Thus, the implemented model-based gradient boosting algorithm is not readily extendable to spatial autoregressive models which are promising candidates for a wide range of applications. Additionally, the algorithm assumes homoskedastic innovations which is a potentially unrealistic assumption in many real-world application settings. In fact, heteroskedasticity is a common feature of spatial datasets lesage2009. Revisiting the average selection rates reveals high FDR values indicating the inclusion of many non-informative variables of minor importance in the final model. Too mitigate the consequences, the utilization of a deselection algorithm is proposed. Alternatively, the application of the so-called stability selection could be utilized which enables the control over the amount of false positive meinshausen2010, shah2013, hofner2015. Furthermore, probing can adapted for the model-based gradient boosting for the SDEM which is currently available for univariate location models. The general idea is to obtain sparser models by stopping the algorithm as soon as the first randomly permuted version of a variable is added thomas2017. Finally, p-values for individual base-learners could be obtained by utilizing permutation techniques hepp2019.
Therefore, the next goal is to extend model-based gradient boosting for the SDEM by incorporating heteroskedastic innovations and adapting it to panel data models with spatially correlated error terms kapoor2007. Such an extension would enable a more nuanced analysis of life expectancy across German districts. For example, the INKAR dataset provides panel data making it possible to study the underlying determinants of life expectancy over time. Beyond the spatial econometrics context, the novel model-based gradient boosting could be extended and applied to network econometric models lee2010. Therefore, practitioners and applied statisticians working with spatial data are encouraged to utilize model-based gradient boosting for the SDEM as a valuable alternative for estimation, regularization, model and variable selection in future research.
\backmatter
\bmhead{Supplementary information} The data as well as the accompanying codebook on the life expectancy in German districts is publicly available via \url{https://www.inkar.de/} (DL-DE BY 2.0). All R-code for the implemented model-based gradient boosting algorithm along with the simulation studies is publicly available in the following GitHub repository \url{https://github.com/micbalz/SpatRegBoost}.
\bmhead{Acknowledgments} The work on this article was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within project 492988838.
\bmhead{Conflict of interests} The author declares that there is no conflict of interest.