Linear or smooth? Enhanced model choice in boosting via deselection of base-learners

With the introduction of newer and more complex model classes (e.g., distributional regression, Rigby and Stasinopoulos, 2005; Kneib, 2013) and also more flexible types of effects, in the last decades the toolbox of data analysts has become richer in potential statistical modelling options. However, this increased flexibility also leads to the burden of having to decide on selecting the best possible option from a vast number of suitable solutions. This is also reflected by the introduction of P-splines (Eilers and Marx, 1996), based also on the work by Brian Marx—who we remember with this special issue. P-splines allow for a very flexible form of effects for continuous variables, but at the same time their flexibility is regulated by penalization in order to ensure smooth effects. This trade-off between flexibility and simplicity is the core of many discussions in statistical modelling between the involved researchers trying to find the best solution for the research question at hand. Generally, the problem of selecting a model as complex as necessary but also as simple as possible can be either tackled from a subject-matter perspective or in a data-driven automated manner. With the emergence of larger data sets, new and modern sources (e.g., proteomics, genomics), there is also an increased need for automated procedures to answer the questions of dimensionality, complexity reduction and model choice.

A natural solution for complexity reduction in the context of statistical modelling are penalized regression techniques. The most popular regularized regression techniques as the lasso (Tibshirani, 1996), ridge regression (Hoerl and Kennard, 1970) or elastic net (Zou and Hastie, 2005) typically aim for linear effects. Boosting, on the other hand, is an indirect approach to regularized regression that can easily incorporate non-linear and smooth effects (as introduced by Schmid and Hothorn, 2008). For the smooth component, the boosting framework relies on the iterative application of P-splines (Eilers and Marx, 1996) as base-learners, that were also extended to monotone or cyclical effects (Hofner et al., 2016). Boosting algorithms originally emerged from machine learning, but were later adapted to statistical modelling (Bühlmann and Hothorn, 2007) and are nowadays a versatile option to fit various model classes, including very complex ones (for an overview see Mayr et al., 2014a, b). For the remainder of this article, we will focus on statistical boosting algorithms in the sense of component-wise gradient boosting with regression models as base-learners (Bühlmann and Hothorn, 2007). It can be shown, that with linear base-learners, boosting with early stopping leads to very similar solutions as the lasso (Hepp et al., 2016).

In the context of dimensionality or complexity reduction, statistical boosting algorithms allow simultaneously: (a) to select the most influential variables from a potentially high-dimensional set of candidate variables; (b) to perform model choice by selecting the most appropriate type of effect via a decomposition into linear and non-linear parts; and (c) to estimate the final prediction model.

The notion of model choice in the context of statistical boosting was introduced by Kneib et al. (2009). The basic idea is to decompose the effect of a continuous variable in a linear component, represented by a linear base-learner, and a smooth component, represented by a non-linear base-learner. To avoid bias towards more complex base-learners, either when choosing between linear and smooth effects or in the presence of categorical variables, Hofner et al. (2011) developed a framework for unbiased model choice by incorporating penalized least squares base-learners.

Over the last years, however, practical experience suggests that there is still a tendency of statistical boosting to select too complex models—particularly in settings with a relatively large number of observations n in comparison to the number of potential predictor variables p (Staerk and Mayr, 2021, Strömer et al., 2022). As overfitting is less problematic in rather low-dimensional settings, the algorithm stops relatively late and has the tendency to include many base-learners. The effect of some of these base-learners in the final model might be only very small, as they were likely only updated once or twice, however this behaviour limits the interpretability of the prediction model without having a relevant impact on prediction performance. This carries over to the question of model choice.

Further regularization might be indicated if variable selection or model choice is key. An approach which can be applied to many regularization techniques is stability selection (Meinshausen and Bühlmann, 2010, Shah and Samworth, 2013), which was later successfully transferred to boosting (Hofner et al., 2015; Mayr et al., 2016; Thomas et al., 2018). This approach is computationally expensive and it remains unclear how to fit the final model.

Here, we propose to tackle this issue with a similar approach as described by Strömer et al. (2022) in the context of variable selection: We directly enforce model choice by deselecting more complex base-learners that do not contribute enough to the final model. In other words, if the smooth component of a continuous predictor was selected, but does not have a larger impact on the predictive risk than the linear one, we can also deselect it and re-fit the model only with the simpler linear component. We investigate the performance of this pragmatic and simple approach in a small simulation study and illustrate its merits by modelling thyroid hormone levels from a large Galician cohort based on different continuous and categorical predictor variables.

Our approach to model choice is based on the application of boosting which emerged from machine learning (Freund and Schapire, 1996) performing classification based on simple decision trees, but was later adapted to estimate regression models via gradient descent in function space (Friedman et al., 2000; Friedman, 2001). For given observations $(y_i\mathrm{,}{x}_i^T)$, $i=1,\dots \mathrm{,}n$, and loss function $\rho$, for example, the common $L_2$ loss $\rho (y\mathrm{,}f(x))=(y-f(x{))}^2$, the component-wise boosting algorithm aims to minimize an empirical risk

by estimating a statistical model $f$. While the loss function defines the type of regression setting (e.g., $L_2$ leading to least-squares mean regression, $L_1$ to median regression), the algorithm allows for additional (pre-)specification of base-learners $h_{\text{}j}(x_l),j=1,...,J$ for each dependent variable $x_l\mathrm{,}l=1,...,p$ according to their potential effect on the outcome $y$, for example, linear, non-linear, smooth or spatial base-learners. The negative gradient vector $u$ of the loss function $\rho (y\mathrm{,}f)$ is evaluated in iteration $m$ with the fitted model at the previous iteration ${\widehat{f}}^{[m-1]}(x)$

and afterwards every base-learner is fitted to the negative gradient vector $u^{[m]}$. Starting with an offset value ${\widehat{f}}^{[0]}$, for example, the mean of $y$ for mean regression, the algorithm fits all specified base-learners to the negative gradient vector of the loss function $u^{[m]}$ (which corresponds to the model fit ${\widehat{f}}^{[m-1]}$ from the previous iteration—see equation (2.1)) and chooses only the best performing base-learner $h_{j^{*}}$. The model is then subsequently updated by ${\widehat{f}}^{[m]}={\widehat{f}}^{[m-1]}+\nu h_{j^{*}}$ in component $j^{*}$ with a fixed step length $0<\nu \le 1$ until the algorithm reaches a final stopping iteration $m_{\text{stop}}$. The step-length is typically chosen to be small, for example, $\nu =0.1$. The main tuning parameter is the stopping iteration $m_{\text{stop}}$, which controls the complexity of the final model. Higher numbers of $m_{\text{stop}}$ hence lead to larger and more complex models, while stopping the algorithm early leads to sparser and simpler models. Usually, the optimal stopping iteration is determined via cross-validation, subsampling or bootstrapping.

Note that in most cases, the number of base-learners $J$ equals the number of candidate variables $p$ (one base-learner for each variable) but it is possible to specify multiple base-learners per variable in the case of a model including both main effects and interactions, or for the decomposition of linear and non-linear effects as we will discuss in this work. However, one could also specify a single base-learner based on multiple variables, for example, in the case of spatial effects. As in each boosting iteration only the best-fitting base-learner is used, that is, either the effect of one base-learner $h_{\text{}j}$ is updated or a new one is incorporated in the model, one can perform automated variable selection when the algorithm is stopped before convergence. Base-learners and their corresponding predictor variables that have never been selected up to that point are effectively excluded from the final model.

We have proposed a simple and pragmatic approach to enhance model choice regarding the decision whether to include a continuous variable with a linear or smooth effect in a boosted statistical model. We have illustrated both in a simulation study as well as on thyroid hormone data that our approach can help researchers to decide on the type of effect and at least in the considered settings led to very promising results: The resulting models were sparser and incorporated only non-linear effects, which were really necessary, but led basically to the same fit and prediction accuracy.

There are, however, several points and limitations to consider: First, our approach actively changes the optimal model selected by the boosting procedure. As the boosting algorithm is typically tuned for predictive risk, it hence should be expected that deselecting components that were initially selected will yield some kind of loss in prediction accuracy. The initial deselection step is controlled by the $\tau$ parameter, which is set to 0.01 (as proposed by Strömer et al., 2022). Tuning of this parameter is non-trivial, as on average the best prediction accuracy should be achieved without any deselection. Second, our model choice approach is designed to favor the simple linear component. The smooth component is only chosen if it outperforms the linear part w.r.t. risk reduction; in any other case we stay with the linear effect. The argument for this is the general notion that in doubt one should always opt for the simplest model possible—but this obviously might not be favourable for all research questions. Also this choice (setting a threshold of 0.5) is not tuned, and could lead to a loss of prediction accuracy. Third, the deselection procedure leads to a longer runtime as the model effectively needs to be fitted twice. Fourth, our approach is focusing on data-driven, automated variable selection and model choice. In many modelling situations, it might be more reasonable to decide based on subject matter knowledge or practical implications (what should the model be used for) on the particular type of effect for a variable, not based on an automated procedure. Fifth, our approach is simple and pragmatic—but not based on theoretical insights. Also inherited from the iterative updating of the boosting algorithm, there is no closed form solution to determine standard errors for effect estimates. The construction of confidence intervals or even significance testing can hence only be performed via resampling techniques (Hofner et al., 2016; Mayr et al., 2017; Hepp et al., 2019).

With all that in mind, we still have reason to believe that this simple procedure might be a valuable option for practical data modelling. Automated model choice is an often cited and highlighted feature of statistical boosting, but the practical relevance over the last years seemed to be limited. Combined with the proposed pragmatic deselection procedure, the decomposition of effects in linear and non-linear components could become a true asset—not only for the simple Gaussian models presented here. The development of more complex and flexible model classes (Kneib et al., 2021) calls also for methods to reduce the complexity again in order to keep the model manageable for data analysts and interpretable for subject matter researchers.

Further research is warranted on how this procedure performs in multi-dimensional optimization problems where the same continuous variable could enter in different model components (like for location and scale, Mayr et al., 2012). Another field of future research might be to extend the procedure not only for selecting the type of effect, but also to allocate variables to different model components like in joint models for longitudinal and time-to-event data (Waldmann et al., 2017; Rappl et al., 2022).

This article would not have been possible without two fellow researchers and friends that are no longer with us and are deeply missed. The work of Professor Brian D. Marx ($†$ 25 November 2021) on P-splines laid the foundation for how we model smooth effects today. With his enthusiasm and support Brian also triggered our commitment to statistical modelling in general. Professor Carmen Cadarso Suárez ($†$ 3 June 2022) was still with us when we started to work on this project and participated in the first meetings. Only Carmen’s collaborative way to work and her talent to bring people together made our collaboration and the work on the Galician cohort possible.