# Stagewise Boosting Distributional Regression

Mattias Wetscher
University of Innsbruck &Johannes Seiler
University of Innsbruck  0000-0002-9982-3001 0000-0001-5714-9234 Reto Stauffer
University of Innsbruck &Nikolaus Umlauf
University of Innsbruck  0000-0002-3798-5507 0000-0003-2160-9803
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Mattias Wetscher, Johannes Seiler, Reto Stauffer, Nikolaus Umlauf Stagewise
Boosting Distributional Regression Forward stagewise regression is a simple
algorithm that can be used to estimate regularized models. The updating rule
adds a small constant to a regression coefficient in each iteration, such that
the underlying optimization problem is solved slowly with small improvements.
This is similar to gradient boosting, with the essential difference that the
step size is determined by the product of the gradient and a step length
parameter in the latter algorithm. One often overlooked challenge in gradient
boosting for distributional regression is the issue of a vanishing small
gradient, which practically halts the algorithm’s progress. We show that
gradient boosting in this case oftentimes results in suboptimal models,
especially for complex problems certain distributional parameters are never
updated due to the vanishing gradient. Therefore, we propose a stagewise
boosting-type algorithm for distributional regression, combining stagewise
regression ideas with gradient boosting. Additionally, we extend it with a
novel regularization method, correlation filtering, to provide additional
stability when the problem involves a large number of covariates. Furthermore,
the algorithm includes best-subset selection for parameters and can be applied
to big data problems by leveraging stochastic approximations of the updating
steps. Besides the advantage of processing large datasets, the stochastic
nature of the approximations can lead to better results, especially for
complex distributions, by reducing the risk of being trapped in a local
optimum. The performance of our proposed stagewise boosting distributional
regression approach is investigated in an extensive simulation study and by
estimating a full probabilistic model for lightning counts with data of more
than 9.1 million observations and 672 covariates. Distributional Regression,
Stagewise Regression, Gradient Boosting, Variable Selection, Correlation
Filtering, Batchwise Updating Distributional Regression, Stagewise Regression,
Gradient Boosting, Variable Selection, Correlation Filtering, Batchwise
Updating Mattias Wetscher, Johannes Seiler, Reto Stauffer, Nikolaus Umlauf
Department of Statistics
Faculty of Economics and Statistics
Universität Innsbruck
6020 Innsbruck, Austria
E-mail:
E-mail: ,
E-mail: ,
E-mail:
URL: https://nikum.org/,
URL: https://retostauffer.org/

##  1 Introduction

Modern regression models can not only provide estimates for the expectation,
but full probabilistic predictions, which are particularly important in
numerous applications, e.g., for the prediction of severe weather using a
complex count model for the number of lightning strikes as demonstrated in
this article (similar models where used in, Simon _et al._ , 2018, 2019). A
well-known model class for such probabilistic settings is the generalized
additive model for location, scale and shape (GAMLSS; Rigby and Stasinopoulos,
2005; Klein _et al._ , 2015b), where each parameter of the response
distribution is modeled by covariates. In order to estimate a well-calibrated
forecast, a stable estimation algorithm is needed for GAMLSS, which is also
able to perform variable selection at the same time. A common choice in such
situations is the gradient boosting algorithm for GAMLSS (Hofner _et al._ ,
2021), which can perform variable selection for parameters of the response
distribution in high dimensions, i.e., even if the problem contains more
covariates than observations.

Gradient boosting in the aforementioned framework iteratively improves the
regression coefficients, thereby fitting the covariates to the negative
gradient of the log-likelihood function with respect to the linear predictors
and only updating the best performing covariate, either for each parameter of
the distribution, also called _cyclical update_ , or only the best performing
parameter, called _non-cyclical update_ (Thomas _et al._ , 2018). Essentially,
this updating scheme is a coordinate descent step, where in each iteration
only the covariate with the largest (in absolute value) partial derivative
gets updated. It is performed at a learning rate such that the objective
function improves very slowly from iteration to iteration, allowing for stable
estimates with shrinkage.

Figure 1:  Illustration of the vanishing gradient problem in boosting when
estimating distributional regression models. The first and second row show an
updating procedure based on the gradient. Row three and four show the
stagewise counterpart where a constant updating step size $\epsilon=0.01$ is
used. The first and third row show the (marginal) log-likelihood functions
$\ell(\beta_{\sigma})$ of the
$\texttt{NO}(\mu=\beta_{\mu},\log(\sigma)=\beta_{\sigma})$ in the left column
and the $\texttt{NBI}(\log(\mu)=\beta_{\mu},\log(\sigma)=\beta_{\sigma})$ in
the right column. The second and fourth row show the corresponding gradient
information. The true parameters ($\beta_{\sigma}^{\star}$) are indicated by a
vertical gray dashed line.

A particular disadvantage of gradient boosting is the search for the optimal
stopping iteration of the algorithm, since too many iterations lead to
overfitting. This problem is usually solved by performing cross-validation or
a computational more efficient method, such as probing (Thomas _et al._ ,
2017), which expands the original data with independent “shadow variables”. If
such a shadow variable is selected, the updating algorithm stops. However, in
the context of distributional regression, this procedures can easily lead to
suboptimal models, especially when gradient boosting is applied to complex
distributions, since the updating scheme depends heavily on the gradient
information. This is the case when the iterations for some parameters quickly
fall into regions of the objective function that are particularly flat, i.e.,
when the gradients become vanishingly small. As a result, some parameters of
the distribution may never be updated for both procedures, and using cross-
validation, the number of boosting iterations would be extremely high to find
the optimum.

The problem of parameters that are never selected for updating is discussed in
Zhang _et al._ (2022), where an adaptive step length is proposed and
analytical optimal step length parameters are provided for the normal location
scale model. The selection of the optimal step length is based on minimizing
the negative log-likelihood for each base learner. This step length is then
multiplied by a penalty $\lambda$ (e.g., $\lambda=0.1$) to again force slow
improvements in the objective function. Despite these efforts, the utilization
of gradients in this technique may still result in the occurrence of vanishing
gradients in complex distributional regression settings, leading to a
selection problem with a possibly high proportion of false negatives.

To address the issues related to vanishing gradients in boosting
distributional regression models, we propose to adapt forward stagewise
regression (Tibshirani, 2015). This updating scheme is very similar to
gradient boosting, but only uses information of the sign of the gradient to
update coefficients at a small constant rate $\epsilon$ (e.g.,
$\epsilon=0.01$). Tibshirani (2015) shows the efficiency of the approach in a
general framework for convex optimization. Since the step length of the update
is kept constant, a stagewise algorithm tailored for distributional regression
does not suffer from the vanishing gradient, and a more balanced choice of
covariates for the update is possible if $\epsilon$ is sufficiently small.

The vanishing gradient problem in distribution regression is illustrated in
Figure 1, where we compare the convergence behavior of the parameters of a
gradient-based updating scheme and our proposed stagewise boosting updating
scheme for the normal distribution (NO($\mu$, $\sigma$)) and the negative
binomial distribution of type one (NBI($\mu$, $\sigma$)). Let us assume that
the algorithms for both the NO and the NBI distribution have already
approximately reached the optimum for the first parameter $\mu$, i.e., that no
further update is required for $\mu$. In this case, the algorithms will only
update the $\sigma$ parameter. This behavior becomes apparent in simple
numerical examples with the NBI distribution. Nonetheless, it is noteworthy
that employing the NO distribution can also result in similar convergence
patterns. As can be seen in Figure 1, the subsequent update iterations for the
second parameter $\sigma$ of the NO distribution converge very quickly to the
optimum $\beta^{\star}_{\sigma}$ with both methods, the gradient information
only becomes very small near the optimum. However, the subsequent gradient
boosting updates for parameter $\sigma$ of the NBI distribution do not reach
the optimum after 1000 iterations and are already slowed down after the first
few iterations due to the very small gradient, as indicated by the updating
points. In contrast, our proposed stagewise updating for distributional
regression does not contain the full information of the gradient and instead
uses a (semi-)constant step length so that it can converge to the true
parameter $\beta^{\star}_{\sigma}$ after a few iterations. This behavior of
gradient boosting algorithms for distributional regression implies the
following: In a larger regression setting with a number of covariates, some
coefficients may never be updated (i.e., remain zero) because distribution
parameters with larger gradients dominate parameters with vanishingly small
gradients in the updating process. In this case, users would have to run the
boosting algorithm for a very long time to detect relevant variables. In
practice, however, this problem goes unnoticed because users typically only
look at the evolution of the log-likelihood to see if the boosting algorithm
has “converged”. And since the gradients and update steps are very small, one
concludes that no further improvement in the log-likelihood can be achieved.
This behavior is well illustrated in Figure 1 by the colored update steps,
which show a sharp break (far away from the optimum $\beta_{\sigma}^{\star}$),
so to speak, resulting in a very flat (almost horizontal) line when looking at
the log-likelihood curve over iterations. Similarly, this behavior is usually
not seen when looking at individual coefficient paths. Finally, this effect
naturally leads to high false negative rates (weak parameters) and high false
positive rates (dominant parameters) in settings with a large number of
potential variables to choose from (see simulation Section 5).

In order to address the aforementioned challenges associated with the use of
gradient boosting for distributional regression, this paper presents a novel
approach to distributional regression by adapting and expanding upon the
forward stagewise algorithm by

  * •

best-subset selection of distributional parameters in each iteration with a
semi-constant step length,

  * •

simple selection of stopping iterations by a novel correlation filtering
approach,

  * •

a batchwise variant of the algorithm for big data, based on stochastic
approximations of the updating steps, mitigating the likelihood of being
trapped in local optima.

Furthermore, we present an extensive simulation study in which we test the
proposed methods in comparison to established state-of-the-art methods in
highly challenging scenarios. We use two two-parameter distributions, a normal
($\texttt{NO}(\mu,\sigma)$) and a gamma distribution
($\texttt{GA}(\mu,\sigma)$), and a three-parameter distribution, the zero-
adjusted negative binomial distribution of type I
($\texttt{ZANBI}(\mu,\sigma,\nu)$), with different numbers of uninformative
variables (30 and 100) for each distributional parameter and different degrees
of multicollinearity between the variables.

To emphasize the usefulness of the approach, our scalable version of the
algorithm is illustrated by a prediction problem of lightning counts using a
high-dimensional dataset with approximately $9.1$ million observations and
$672$ variables. The prediction is based on the ZANBI distribution, where each
of the three distribution parameters is modeled by the variables selected by
our algorithm.

The remainder of this paper is organized as follows: In Section 2 we give a
brief overview of distributional regression models. Gradient boosting for
distributional regression is described in Section 3. Stagewise regression in
general as well as our novel variants for distributional regression using
semi-constant step length, best-subset updating, correlation filtering and
batchwise estimation are presented in Section 4. In Section 5 we present the
results of the simulation study, in Section 6 we demonstrate the application
of our methods to the complex problem of modeling lightning counts.

##  2 Distributional Regression

The underlying idea of distributional regression (Rigby and Stasinopoulos,
2005; Klein _et al._ , 2015a; Umlauf and Kneib, 2018) is to model all
parameters of an arbitrary parametric response distribution (rather than just
the mean) through covariates. In the following, we briefly introduce the
framework and describe classic maximum likelihood estimation of the
parameters.

###  2.1 Model Specification

For a dataset with $i=1,\ldots,n$ observations, where $y_{i}$ denotes the
response (which may also be non-continuous or multivariate) and
$\mathbf{x}_{i}$ available covariate information, we assume conditional
independence of the individual response observations given covariates.
Specifically, let

$y_{i}\sim\mathcal{D}\left(\theta_{1}(\mathbf{x}_{i}),\dots,\theta_{K}(\mathbf{x}_{i})\right),$

where $\mathcal{D}(\cdot)$ represents a distribution with $K$ parameters
$\theta_{ik}\equiv\theta_{k}(\mathbf{x}_{i})$, $k=1,\ldots,K$, that are linked
to additive predictors $\eta_{ik}\equiv\eta_{k}(\mathbf{x}_{i})$ using known
monotonic and twice differentiable functions $h_{k}(\cdot)$,

$h_{k}(\theta_{ik})=\eta_{ik}=\mathbf{x}_{ik}^{\top}\boldsymbol{\beta}_{k}.$

Here, $\mathbf{x}_{ik}$ is a row of the $n\times J_{k}$ predictor specific
design matrix $\mathbf{X}_{k}$ and the parameters
$\boldsymbol{\beta}_{k}=(\boldsymbol{\beta}_{0k},\ldots,\boldsymbol{\beta}_{J_{k}k})^{\top}$
are regression coefficients. In matrix notation, the predictors are written as

$\boldsymbol{\eta}_{k}=(\eta_{1k},\ldots,\eta_{nk})^{\top}=\mathbf{X}_{k}\boldsymbol{\beta}_{k},$

which is the fully parametric case of the GAMLSS (Rigby and Stasinopoulos,
2005).

###  2.2 Classic Likelihood Estimation

Parameter estimation is based on maximizing the log-likelihood function

$\ell(\boldsymbol{\beta};\mathbf{y},\mathbf{X})=\sum_{i=1}^{n}\log\,d_{y}\left(y_{i};\theta_{i1}=h_{1}^{-1}(\eta_{i1}),\ldots,\theta_{iK}=h_{K}^{-1}(\eta_{iK})\right),$
(1)

where $d_{y}(\cdot)$ is the density function and
$\mathbf{y}=(y_{1},\ldots,y_{n})^{\top}$ the response vector,
$\boldsymbol{\beta}=(\boldsymbol{\beta}_{1}^{\top},\ldots,\boldsymbol{\beta}_{K}^{\top})^{\top}$
the stacked vector of regression coefficients to be estimated and
$\mathbf{X}=(\mathbf{X}_{1},\ldots,\mathbf{X}_{K})$ represents the full
covariate data matrix.

For the maximization of (1), Rigby and Stasinopoulos (2005) use a modified
backfitting algorithm based on iteratively reweighted least squares (IWLS;
Gamerman, 1997). Similarly, Umlauf _et al._ (2018) use IWLS based iterations
for both, maximum likelihood and full Bayesian estimation by Markov chain
Monte Carlo (MCMC) simulation.

An important aspect in estimating distributional regression models often
involves variable selection. Classical maximum likelihood methods, due to
their cyclic updating nature, are not suitable for selection and
regularization. In contrast, boosting updates inherently perform variable
selection through early stopping, rendering it more suitable for this purpose.
Another shrinkage method is presented by Groll _et al._ (2019) based on
$L_{1}$-type regularization in the context of distributional regression for
metric covariates and both group and fused-LASSO for categorical variables.
The major drawback here, however, is that the algorithm allows little
flexibility and is additionally numerically very expensive as the shrinkage
parameters must be determined using a grid search. This makes the method
unattractive for many applications, e.g., models for very large data like in
our application (see Section 6). Therefore, boosting-type algorithms are
currently considered the leading methods for variable selection in the context
of distributional regression. This is largely attributed to their exceptional
numerical stability, which makes them particularly effective for this purpose,
as described in more detail in the following section.

##  3 Gradient Boosting

Gradient boosting is a widely used supervised machine learning technique that
is particularly effective in constructing predictive models, including
distributional regression models. It is a type of boosting algorithm that aims
to improve the prediction of the model by iteratively minimizing the residuals
or a general loss function in very small steps. In contrast to the traditional
backfitting step, which updates every coefficient in each iteration, gradient
boosting for distributional regression follows a different approach by only
updating the best-performing regression coefficient overall (non-cyclical
update) or the best-performing regression coefficient for each distributional
parameter (cyclical update) in each iteration (Hofner _et al._ , 2021; Thomas
_et al._ , 2018)111Hofner _et al._ , 2021 and Thomas _et al._ , 2018 formulate
gradient boosting with more general effects, i.e., baselearners and in the
context of minimizing a loss function. Instead of the latter we present it in
the context of maximizing a log-likelihood.. The nature of the iterative
updating also makes gradient boosting a powerful technique for modeling
complex data even with only limited sample size.

###  3.1 Cyclical Gradient Boosting

For the following assume that each covariate in $\mathbf{x}_{i}$ is
standardized. To illustrate the concept of cyclical gradient boosting (Hofner
_et al._ , 2016), let us consider a linear model given by

$y_{i}=\mathbf{x}_{i}^{\top}\boldsymbol{\beta}+\varepsilon_{i}\quad\text{and}\quad\varepsilon_{i}\sim\mathcal{N}(0,\sigma^{2}),$

which can be seen as a special case of the more general distributional
regression model:

$y_{i}\sim\mathcal{D}\left(y_{i};\theta_{i1}=\eta_{i1},\theta_{i2}=\exp(\eta_{i2})\right),$

where $\mathcal{D}(\cdot)$ represents the two-parameter normal distribution
$\mathcal{N}(\mu,\sigma^{2})$, and $\theta_{i1}$ and $\theta_{i2}$ correspond
to the mean parameter $\mu$ and variance parameter $\sigma^{2}$, respectively.
In the homoskedastic case, we have $\exp(\eta_{i2})=\sigma^{2}\,\,\forall\,i$.
In the full distributional regression setting, we model
$\exp(\eta_{i2})=\sigma_{i}^{2}$, hence, we specify the linear predictors
$\eta_{i1}=\mathbf{x}_{i1}^{\top}\boldsymbol{\beta}_{1}$ and
$\eta_{i2}=\mathbf{x}_{i2}^{\top}\boldsymbol{\beta}_{2}$.

In gradient boosting for distributional regression the intercepts, for the
linear model $\beta_{01}$ and $\beta_{02}$, are initialized by its
corresponding maximum likelihood estimates, which is usually relatively
straightforward to compute. All other coefficients are initialized as zero.

Then, starting from the initial values

$\boldsymbol{\theta}_{k}^{[0]}=h_{k}^{-1}\left(\boldsymbol{\eta}_{k}^{[0]}\right)=h_{k}^{-1}\left(\mathbf{X}_{k}\boldsymbol{\beta}_{k}^{[0]}\right),$

in each iteration $t=1,\ldots,T$, of the boosting algorithm, the coefficients
$\boldsymbol{\beta}_{1}^{[t]}$ and $\boldsymbol{\beta}_{2}^{[t]}$ are improved
sequentially and slowly. Starting from $\boldsymbol{\beta}_{1}^{[t]}$, fit
each column of the predictor-specific model matrix $\mathbf{X}_{1}$, denoted
by $({\mathbf{X}_{1}})_{\cdot j_{1}}$, to the gradient vector

$\mathbf{g}_{1}=\left(\frac{\partial\,\log\,d_{y}(y_{i};\theta_{i1}^{[t-1]},\theta_{i2}^{[t-1]})}{\partial\eta_{1}}\right)_{i=1,\dots,n}$

and select the coefficient to be updated according to the residual sum of
squares (RSS) criterion

$j_{1}^{*}=\underset{{j_{1}=1,\dots,J_{1}}}{\mathrm{argmin}}\left\|\mathbf{g}_{1}-({\mathbf{X}_{1}})_{\cdot
j_{1}}c_{j_{1}}\right\|_{2}^{2},$

where $c_{j_{1}}$ represents the coefficient of the least squares fit of
column $({\mathbf{X}_{1}})_{\cdot j_{1}}$ with the current gradient vector
$\mathbf{g}_{1}$. The updating step then is

$\boldsymbol{\beta}_{1}^{[t]}=\boldsymbol{\beta}_{1}^{[t-1]}+\epsilon\cdot
c_{j_{1}^{*}}\cdot\mathbf{e}_{j_{1}^{*}},$

where $\epsilon$ is the step length control parameter (e.g., set to
$\epsilon=0.1$) and $\mathbf{e}_{j_{1}^{*}}$ is an index vector consisting of
zeros, except for the position $j_{1}$, where it contains a one. This updating
step moves the coefficients slowly in the direction of $\mathbf{g}_{1}$, which
is similar to a coordinate descent step and computationally efficient, since
least squares estimation is fast and numerically stable. Secondly, the process
is repeated with $\boldsymbol{\beta}_{2}^{[t]}$, but the updated version
$\boldsymbol{\beta}_{1}^{[t]}$ is used, i.e., each column of
$({\mathbf{X}_{2}})_{\cdot j_{2}}$ is fitted to the gradient vector

$\mathbf{g}_{2}=\left(\frac{\partial\,\log\,d_{y}(y_{i};\theta_{i1}^{[t]},\theta_{i2}^{[t-1]})}{\partial\eta_{2}}\right)_{i=1,\dots,n}.$

Then, $j_{2}^{*}$ is determined with the respective RSS-criterion and the
update

$\boldsymbol{\beta}_{2}^{[t]}=\boldsymbol{\beta}_{2}^{[t-1]}+\epsilon\cdot
c_{j_{2}^{*}}\cdot\mathbf{e}_{j_{2}^{*}},$

is carried out.

In more generality, with $K$ distibutional paramerters, the updating is
sketched as follows:

$\displaystyle\mathbf{g}_{1}$
$\displaystyle=\left(\frac{\partial\,\log\,d_{y}(y_{i},\mathbf{x}_{i},(\boldsymbol{\beta}_{1}^{[t-1],\top},\ldots,\boldsymbol{\beta}_{K}^{[t-1],\top})^{\top})}{\partial{\eta_{1}}}\right)_{i=1,\dots,n}$
$\displaystyle\underset{\text{update}}{\rightarrow}$
$\displaystyle\boldsymbol{\beta}_{1}^{[t]}$ $\displaystyle\mathbf{g}_{2}$
$\displaystyle=\left(\frac{\partial\,\log\,d_{y}(y_{i},\mathbf{x}_{i},(\boldsymbol{\beta}_{1}^{[t],\top},\boldsymbol{\beta}_{2}^{[t-1],\top},\ldots,\boldsymbol{\beta}_{K}^{[t-1],\top})^{\top})}{\partial{\eta_{2}}}\right)_{i=1,\dots,n}$
$\displaystyle\underset{\text{update}}{\rightarrow}$
$\displaystyle\boldsymbol{\beta}_{2}^{[t]}$ $\displaystyle\vdots$
$\displaystyle\mathbf{g}_{K}$
$\displaystyle=\left(\frac{\partial\,\log\,d_{y}(y_{i},\mathbf{x}_{i},(\boldsymbol{\beta}_{1}^{[t],\top},\ldots,\boldsymbol{\beta}_{K-1}^{[t],\top},\boldsymbol{\beta}_{K}^{[t-1],\top})^{\top})}{\partial{\eta_{K}}}\right)_{i=1,\dots,n}$
$\displaystyle\underset{\text{update}}{\rightarrow}$
$\displaystyle\boldsymbol{\beta}_{K}^{[t]}.$

###  3.2 Non-Cyclical Gradient Boosting

Instead of cyclically updating all parameters, a variation of the gradient
boosting algorithm focuses on updating only the best model term. This is done
by selecting among the tentative updates

$\boldsymbol{\beta}_{k}^{[t]}=\boldsymbol{\beta}_{k}^{[t-1]}+\epsilon\cdot
c_{j_{k}^{*}}\cdot\mathbf{e}_{j_{k}^{*}},$ (2)

with

$j_{k}^{*}=\underset{{j_{k}=1,\dots,J_{k}}}{\mathrm{argmin}}\left\|\mathbf{g}_{k}-({\mathbf{X}_{k}})_{\cdot
j_{k}}c_{j_{k}}\right\|_{2}^{2},$

only the one that maximizes the log-likelihood most

$k^{*}=\underset{{k=1,\dots,K}}{\mathrm{argmax}}~{}\ell\left(\left(\boldsymbol{\beta}_{1}^{[t-1]\top},\ldots,\boldsymbol{\beta}_{k}^{[t]\top},\dots\boldsymbol{\beta}_{K}^{[t-1]\top}\right)^{\top};\mathbf{y},\mathbf{X}\right).$

The tentative updates are computed analogous to the cyclical updating but the
respective column of the predictor-specific model matrices $\mathbf{X}_{k}$
are fitted to a not updated version of the gradient,

$\mathbf{g}_{k}=\left(\frac{\partial\,\log\,d_{y}(y_{i};\theta_{i1}^{[t-1]},\ldots,\theta_{iK}^{[t-1]})}{\partial\eta_{k}}\right)_{i=1,\dots,n}$
(3)

The final update in iteration $t$ is only carried out in the $k^{*}$-th
distributional parameter and for the cases $k\neq k^{*}$, the coefficients
remain unchanged, i.e.,
$\boldsymbol{\beta}_{k}^{[t]}=\boldsymbol{\beta}_{k}^{[t-1]}$. In contrast to
cyclic updating, this so-called non-cyclic updating method has proven
advantages and is thoroughly investigated in Thomas _et al._ (2018).

###  3.3 Remarks

The advantages of gradually improving the parameters of the distributional
regression model are manifold, one of which is the ability to apply
regularization and variable selection tech- niques to the evolution of the
coefficients. In the gradient boosting setting, some popular tech- niques
include cross-validation (Hastie _et al._ , 2009), stability selection
(Meinshausen and Bühlmann, 2010; Thomas _et al._ , 2018), and variable
deselection (Strömer _et al._ , 2022). Cross-validation is a common
regularization method used to find an early stopping for the algorithm.
Stability selection is another method that involves randomly dividing the data
into multiple subsets, boosting the regression model on each subset, and then
selecting variables that are consistently selected across different subsets.
However, both methods have the disadvantage of requiring the model to be
calculated multiple times on different subsets, which can drastically increase
computation time. Variable deselection is a variable selection algorithm that
first grows a full model via gradient boosting and then deselects unimportant
variables with low impact on risk reduction. We use these methods in the
simulation in Section 5 as benchmark methods in combination with the non-
cyclical gradient boosting.

For the cyclical updating, cross-validation aims to find an early stopping for
each distributional parameter. However, this method can be computationally
expensive as a multidimensional grid of stopping iterations needs to be
considered, which makes it unsuitable for many (big data) applications.
Therefore, we only compare non-cyclical updating with our new methods in the
simulation.

One of the main challenges in dealing with complex distributions is the
selective saturation of some distributional parameters that correspond to
small gradient vectors $\mathbf{g}_{k^{*}}$. These parameters may be selected
late in the updating process or not at all, while other parameters may already
be saturated with selected variables. As a result, achieving a balanced
selection of variables for different parameters may prove challenging,
potentially resulting in reduced predictive performance and, more importantly,
incorrect conclusions regarding influencing factors. Furthermore, the log-
likelihood function for complex distributions may have multiple local maxima,
which increases the risk of getting stuck in such maxima due to the vanishing
gradient problem. To address these issues, we propose an adaptation of the
general stagewise regression algorithm (Tibshirani, 2015) for distributional
regression. In the next section, we provide a detailed description of our
approach and discuss its improvements.

##  4 Stagewise Boosting

In this section, we adapt the general stagewise regression presented in
Tibshirani (2015) to the distributional regression setting and present several
notable improvements. Please note that in the following we use non-cyclical
updating (see Section 3.2) throughout.

###  4.1 Stagewise Boosting Distributional Regression

We start with adapting the forward stagewise algorithm in Tibshirani (2015)
for a distributional gradient boosting setting similar to gradient boosting
(3.1) as follows. Assume all covariates are standardized, we then replace the
tentative updates in (3.1) with,

$\boldsymbol{\beta}_{k}^{[t]}=\boldsymbol{\beta}_{k}^{[t-1]}+\epsilon\cdot\mathrm{sign}\left(({\mathbf{X}_{k}})_{\cdot
j_{k}^{*}}^{\top}\mathbf{g}_{k}\right)\cdot\mathbf{e}_{j_{k}^{*}},$ (4)

where $\mathbf{g}_{k}$ is defined in (3). Furthermore, we replace the RSS-
criterion with the inner product criterion (IP), meaning the index $j_{k}^{*}$
is selected as the maximizing index in absolut value among $1,\dots,J_{k}$ for
the inner product $({\mathbf{X}_{k}})_{\cdot j_{k}^{*}}^{\top}\mathbf{g}_{k}$,
i.e.,

$j_{k}^{*}=\underset{{j=1,\dots,J_{k}}}{\mathrm{argmax}}\,\left|({\mathbf{X}_{k}})_{\cdot
j}^{\top}\mathbf{g}_{k}\right|.$

Using the linearity of the predictors $\eta_{k}$ and the chain rule, one can
derive the following equality,

$\displaystyle({\mathbf{X}_{k}})_{\cdot j_{k}^{*}}^{\top}\mathbf{g}_{k}$
$\displaystyle=\sum_{i=1}^{n}({\mathbf{X}_{k}})_{ij_{k}^{*}}\frac{\partial\,\log\,d_{y}(y_{i};\theta_{i1}^{[t-1]},\ldots,\theta_{iK}^{[t-1]})}{\partial\eta_{k}}$
$\displaystyle=\sum_{i=1}^{n}\frac{\partial\eta_{ik}}{\partial\beta_{j_{k}^{*}}}\frac{\partial\,\log\,d_{y}(y_{i};\theta_{i1}^{[t-1]},\ldots,\theta_{iK}^{[t-1]})}{\partial\eta_{k}}=\frac{\partial\,\ell(\boldsymbol{\beta}^{[t-1]};\mathbf{y},\mathbf{X})}{\partial\beta_{j_{k}^{*}k}}.$

This means the inner product $({\mathbf{X}_{k}})_{\cdot
j_{k}^{*}}^{\top}\mathbf{g}_{k}$ coincides with the derivative of the log-
likelihood with respect to $\beta_{j_{k}^{*}k}$.

Since we have standardized variables, the correlation $c_{j_{k}k}$ is given by

$c_{j_{k}k}=\frac{1}{n-1}(\mathbf{X}_{k})_{\cdot j_{k}}^{\top}\mathbf{g}_{k}.$

Furthermore, the residual sum of squares (RSS) can be expressed in terms of
the correlation as follows:

$\left\|\mathbf{g}_{k}-(\mathbf{X}_{k})_{\cdot
j_{k}}c_{j_{k}}\right\|_{2}^{2}=\mathbf{g}_{k}^{\top}\mathbf{g}_{k}\left(1-c_{j_{k}k}^{2}\right).$

This shows that the RSS and the inner product (IP) criteria are equivalent.

Thus, the only difference to gradient boosting is that the update step length
is constant, denoted by $\epsilon$. Since a constant step length is used, the
stagewise boosting update rule is not susceptible to the vanishing gradient
problem. However, a disadvantage of (4) is that we lose the flexibility to
take larger update steps when the gradient information would suggest it,
requiring more steps than the usual gradient boosting update (3.1). In
addition, (4) is unable to improve the coefficients in the convergence phase
of the algorithm, when much smaller updates than $\epsilon$ are required. We
fix both problems with further adjustments to (4), namely with a semi-constant
stagewise version similar to the gradient clipping in Bengio _et al._ (1994).
Moreover, after $\rho\cdot 100$% (e.g., 80%) of iterations, we no longer
impose a minimum length of the update step, allowing the algorithm to converge
extremely fast.

Gradient (or partial derivative) clipping is a technique often used when
training deep learning models to prevent gradients from becoming too large or
too small. It involves rescaling the gradients so that their (Euclidean-)
length is within a certain range. This is usually done by setting a maximum
(minimum) threshold for the gradients. If the length of a gradient is larger
(smaller) than this value, it is scaled so that it has the same length as the
threshold value.

##### Semi-constant step length

More precisely, we substitute the constant $\epsilon$ in (4) with the semi-
constant version $\epsilon_{j^{*}_{k}k}$, which depends on the partial
derivative of the normalized log-likelihood
$\ell(\boldsymbol{\beta};\mathbf{y},\mathbf{X})/n$ with respect to the
$j^{*}_{k}\text{-th}$ variable in the $k\text{-th}$ distributional parameter,
denoted by $\beta_{j^{*}_{k}k}$222Instead of the sum, we use the mean to
compute the gradients to avoid dependence on the sample size.. With this, we
arrive at the updating rule

$\displaystyle\partial\ell_{j^{*}_{k}k}$
$\displaystyle:=\frac{1}{n}\frac{\partial\ell(\boldsymbol{\beta}^{[t-1]};\mathbf{y},\mathbf{X})}{\partial{\beta_{j^{*}_{k}k}}},$
$\displaystyle\epsilon_{j^{*}_{k}k}$
$\displaystyle=\begin{cases}\nu\cdot\epsilon&\text{if
}|\partial\ell_{j^{*}_{k}k}|<\nu\cdot\epsilon\text{ and }t<\rho\cdot T\\\
|\partial\ell_{j^{*}_{k}k}|&\text{if
}|\partial\ell_{j^{*}_{k}k}|<\nu\cdot\epsilon\text{ and }t\geq\rho\cdot T\\\
|\partial\ell_{j^{*}_{k}k}|&\text{if
}\nu\cdot\epsilon\leq|\partial\ell_{j^{*}_{k}k}|\leq\epsilon\\\
\epsilon&\text{else},\end{cases}$ (SC-SDR)
$\displaystyle\boldsymbol{\beta}_{k}^{[t]}$
$\displaystyle=\boldsymbol{\beta}_{k}^{[t-1]}+\epsilon_{j^{*}_{k}k}\cdot\mathrm{sign}\left((\mathbf{X}_{k})_{\cdot
j^{*}_{k}}^{\top}\mathbf{g}_{k}\right)\cdot\mathbf{e}_{j^{*}_{k}}$
$\displaystyle j_{k}^{*}$
$\displaystyle=\underset{{j=1,\dots,J_{k}}}{\mathrm{argmax}}\,\left|({\mathbf{X}_{k}})_{\cdot
j}^{\top}\mathbf{g}_{k}\right|,$

where $\nu$ (e.g., 0.1) is a shrinkage parameter used to define the range of
the update. Essentially, our update stays between
$\nu\cdot\epsilon_{j^{*}_{k}k}\cdot\mathrm{sign}\left(({\mathbf{X}_{k}})_{\cdot
j^{*}_{k}}^{\top}\mathbf{g}_{k}\right)$ and
$\epsilon_{j^{*}_{k}k}\cdot\mathrm{sign}\left(({\mathbf{X}_{k}})_{\cdot
j^{*}_{k}}^{\top}\mathbf{g}_{k}\right)$ in the first $\rho\cdot 100$% of the
iterations. After that, we omit the minimum threshold and the update stays
between $0$ and
$\epsilon_{j^{*}_{k}k}\cdot\mathrm{sign}\left(({\mathbf{X}_{k}})_{\cdot
j^{*}_{k}}^{\top}\mathbf{g}_{k}\right)$. As mentioned earlier, this makes it
possible for the algorithm to converge properly. The tag of 4.1 stands for
semi-constant updating for stagewise distributional regression.

##### Illustration

We demonstrate the effectiveness of 4.1 by using the example on chronic
undernutrition of preschool children taken from Zhang _et al._ (2022). In this
manuscript, the authors develop a semi-analytical adaptive step-length (SAASL)
boosting algorithm that aims to address the problem of unequal updating in
distributional regression for the normal location-scale model in the context
of non-cyclical gradient boosting. The updating step length therein is 10% of
the optimal updating step length with regard to the risk reduction of the
respective parameters. We also compare our approach to the non-cyclical
gradient boosting algorithm as implemented in the R package gamboostLSS
(Hofner _et al._ , 2021). This application illustrates modeling chronic
undernutrition of preschool children aged 0-35 months—denoted by the variable
stunting—for the same subset used by Zhang _et al._ (2022) of the data taken
from the 1998-99 Indian Demographic and Health Survey (International Institute
for Population Sciences and ORC Macro(2000), IIPS). Similar to Zhang _et al._
(2022) a normal model $\texttt{NO}(\mu=\eta_{\mu},\log(\sigma)=\eta_{\sigma})$
is used, where predictors $\eta_{k}$, $k\in\\{\mu,\sigma\\}$, depend on a
quadratic polynomial of each of the variables cbmi (BMI of the child), cAge
(age of the child), mBMI (BMI of the mother), and mAge (age of the mother) and
are specified as follows:

$\displaystyle\eta_{k}$ $\displaystyle=$
$\displaystyle\beta_{0k}+\beta_{1k}\cdot\texttt{cBMI}+\beta_{2k}\cdot\texttt{caAge}+\beta_{3k}\cdot\texttt{mBMI}+\beta_{4k}\cdot\texttt{mAge}+$
$\displaystyle\beta_{5k}\cdot\texttt{cBMI2}+\beta_{6k}\cdot\texttt{cAge2}+\beta_{7k}\cdot\texttt{mBMI2}+\beta_{8k}\cdot\texttt{mAge2}.$

All variables are standardized and the added 2 at the end of the variable
names denotes the squared effects. In Figure 2, the corresponding coefficient
paths are shown. We observe that both the non-cyclical gradient boosting and
SAASL-boosting algorithm fail to converge within $5000$ iterations, while the
4.1 algorithm converges after $\approx 1300$ iterations. This is due to the
small updating steps that depend on the small gradients in the first two
methods and the clipped updating steps in the latter
($\nu=0.1,\epsilon=0.01,\rho=0.8$). The results demonstrate that stagewise
boosting distributional regression with the updating rule 4.1 can improve the
convergence speed and accuracy compared to other methods.

Why is this important? In a cross-validation scenario, consider the task of
identifying relevant variables. If the boosting algorithm fails to converge,
it may result in a high number of false negatives. While the user might
attempt to address this by increasing the number of iterations, our experience
suggests that doing so can render estimation infeasible. Additionally, users
may not always be aware of this potential issue.

Figure 2:  Comparison of coefficient paths obtained from different boosting
algorithms to illustrate the problem of chronic malnutrition of preschool
children in India. The three methods shown are gradient boosting, SAASL
boosting, and non-cyclical stagewise boosting with updating rule 4.1.

###  4.2 Best-Subset Updating

In the non-cyclical updating scheme, only one distributional parameter is
updated per iteration. In contrast, the cyclical updating scheme updates each
parameter in turn at every iteration, requiring the selection of optimum
stopping iterations for each parameter through a grid search. While the latter
approach is often computationally infeasible, e.g., for stability selection or
big data problems, it still offers some advantages in terms of variable
selection. This is because parameter-wise stopping iterations, or updates, can
be generated for parameters even when encountering issues like vanishing small
gradients, provided that a sufficient number of iterations are specified. To
leverage the advantages of both approaches, we propose a novel updating
variant called _best-subset updating_. This approach is particularly useful
for estimating models with complex distributions, such as ZANBI, where certain
distributional parameters may be dominated by others using non-cyclical
updating, making it impossible to select a balanced set of variables. Best-
subset updating addresses this issue by considering all possible combinations
of distributional parameters for updating, resulting in superior estimation of
all distributional parameters. This approach mirrors the computationally
demanding grid search employed in cyclic boosting algorithms, but instead of
tracking a grid of parameter iterations, only one iteration path is required.
This reduces the early stopping search to a single dimension, which makes the
best-subset updating very useful for big data problems. For example, this can
be seen in the simulation in Section 5, where we compare the true positive
rate of the three parameter distributions in Figure 7.

The algorithm works as follows: For every non-empty subset
$S\subset\\{1,\dots,K\\}$, we combine the partial derivatives corresponding to
$S$ into a gradient vector
$\nabla\mathbf{L}_{S}=\left(\partial\ell_{j^{*}_{s}s}\right)_{s\in S}$ which
will be used to define the updating step. These partial derivatives represent
a weighting of the different variables with regard to their relative
importance. To avoid an exploding or a vanishing gradient, we rescale the
gradient vector $\nabla\mathbf{L}_{S}$ if its Euclidean length exceeds a
certain value $\epsilon$. Additionally, if any individual partial derivative
$\partial\ell_{j^{*}_{s}s}$ fails to overcome a minimum threshold
$\nu\cdot\epsilon$ (e.g., $\nu=0.1$ and $\epsilon=0.01$) in the first
$\rho\cdot 100\%$ of the iterations (e.g., $\rho=0.8$), we rescale them to
this threshold value. The scaling factor controlling the maximum threshold is
defined as

$\displaystyle\epsilon_{S}$ $\displaystyle=\begin{cases}1&\text{if
}\|\nabla\mathbf{L}_{S}\|_{2}<\epsilon\\\
\frac{\epsilon}{\|\nabla\mathbf{L}_{S}\|_{2}}&\text{else}.\end{cases}$

Then, clipping from below for the individual derivatives ensures that the
updating steps do not vanish in the first $\rho\cdot 100$% of the iterations
and the final updating step is

$\displaystyle\epsilon_{Ss}$
$\displaystyle=\begin{cases}\nu\cdot\epsilon&\text{if
}\epsilon_{S}\cdot|\partial\ell_{j^{*}_{s}s}|<\nu\cdot\epsilon\text{ and
}t<\rho\cdot T\\\ \epsilon_{S}\cdot|\partial\ell_{j^{*}_{s}s}|&\text{if
}\epsilon_{S}\cdot|\partial\ell_{j^{*}_{s}s}|<\nu\cdot\epsilon\text{ and
}t\geq\rho\cdot T\\\ \epsilon_{S}\cdot|\partial\ell_{j^{*}_{s}s}|&\text{if
}\nu\cdot\epsilon\leq\epsilon_{S}\cdot|\partial\ell_{j^{*}_{s}s}|\leq\epsilon\\\
\epsilon&\text{else},\end{cases}$ (SC-BS-SDR)
$\displaystyle\boldsymbol{\beta}_{s}^{[t]}$
$\displaystyle=\boldsymbol{\beta}_{s}^{[t-1]}+\epsilon_{Ss}\cdot\mathrm{sign}\left(({\mathbf{X}_{s}})_{\cdot
j^{*}_{s}}^{\top}\mathbf{g}_{s}\right)\cdot\mathbf{e}_{j^{*}_{s}}$
$\displaystyle j_{s}^{*}$
$\displaystyle=\underset{{j=1,\dots,J_{s}}}{\mathrm{argmax}}\,\left|({\mathbf{X}_{s}})_{\cdot
j}^{\top}\mathbf{g}_{s}\right|,$

for $s\in S$. Similar to 4.1, after $\rho\cdot 100$% of the iterations, the
minimal length of the updating step $\epsilon_{Ss}$ is no longer imposed. The
subset $S^{*}$ which maximizes the log-likelihood $\ell(\cdot)$ is choosen for
an update in iteration $t$. The tag of SC-BS-SDR stands for semi-constant,
best subset updating for stagewise distributional regression.

To prevent overfitting, early stopping methods need to be applied to both 4.1
and SC-BS-SDR, as both algorithms gradually build up their coefficients. In
the following, we introduce a novel variable selection method, correlation
filtering, that takes advantage of the correlation of the data with the
gradient vectors in each iteration. This variable selection algorithm can be
used in combination with both 4.1 and SC-BS-SDR.

###  4.3 Early Stopping

Optimal stopping iterations in boosting settings are typically determined
through costly cross-validation, which can be impractical in many scenarios.
To address this, we propose using the Bayesian information criterion (BIC):

$\mathrm{BIC}^{[t]}=-2\cdot\ell(\boldsymbol{\beta}^{[t]};\mathbf{y},\mathbf{X})+\log(n)\cdot\mathrm{df}^{[t]},$

where $\mathrm{df}^{[t]}$ is the number of non-zero parameters at iteration
$t$. This approach usually leads to parsimonious models, since counting the
number of non-zero coefficients is more restrictive compared to estimating the
degrees of freedom using the hat matrix (which is very computationally
intensive). Zou _et al._ (2007) investigate the usefulness of this approach
for LASSO. We adopt this approach for stagewise boosting distributional
regression.

###  4.4 Correlation Filtering

Instead of computing the inner product of the variables with the gradient
vector ${(\mathbf{X}_{k}})_{\cdot j_{k}}^{\top}\mathbf{g}_{k}$ for each
distributional parameter used to determine the candidate variables to be
updated, one can determine the variable to be updated using the correlation
between the gradient vector and the variables if all variables are
standardized. The latter variant has the advantage that it can be regularized
by adding a threshold value $\kappa$ (e.g. $\kappa=0.15$) to filter out
variables with an absolut correlation value
$\left|c_{j_{k}k}\right|\leq\kappa$. Only the remaining variables with
$\left|c_{j_{k}k}\right|>\kappa$ are considered for a possible update. In the
case that no variable remains in a distributional parameter with a
sufficiently large correlation $c_{j_{k}k}$ with $\mathbf{g}_{k}$, then no
candidate variable is provided. In the extreme case where no candidate
variable is provided across all distributional parameters, no update is
carried out. Typically, in the early iterations of the updating, the
correlation between the covariates and $\mathbf{g}_{k}$ are large and with
continuing updates the correlated variables get absorbed into the model,
leading to a decline of the correlation values until they all no longer
overcome the minimum requiremend $\kappa$ for an update. At this stage the
updating stops and an implicit early stopping is found. The correlation
filtering (CF) in combination with SC-BS-SDR is presented in Algorithm SBDR in
Section 4.5.

Please note that CF only serves as a variable selection step. After selecting
the variables, their coefficients are boosted in a refitting step until
convergence. This process resembles the deselection algorithm presented by
Strömer _et al._ (2022), where a complete model, incorporating all variables,
is grown using boosting333utilizing non- cyclical gradient boosting in the
distributional regression setting. Subsequently, insignificant variables with
minimal impact on risk reduction are deselected, and the model is re-estimated
using only the selected variables (algorithm VarDes in the following).

To exemplify the usefulness of CF, a comparison between VarDes and CF with
best-subset selection is made in Figure 3.

Figure 3:  Comparison of variable selection steps of VarDes (top) vs. CF
(bottom). ZANBI distribution with 1000 observations and 2000 noise variables
for parameters $\mu$, $\sigma$ and $\nu$. The true non-noise variables for
$\mu$ are $x_{1},x_{3},x_{5},x_{6}$, for $\sigma$ are $x_{2},x_{4},x_{5}$ and
for $\nu$ are $x_{3},x_{4},x_{5}$. All variables are correlated according to
the correlation setting in the simulation Section 5 ($\rho=0.7$). The
vertical-gray-dotted line represents the early stopping induced by 10-fold
cross validation for the VarDes and BIC by CF. The displayed variable names
are the final set of selected variables for the respective predictors.

This case shows a high $p$ low $n$ setting ($p>n$, we simulate data according
to the data generating process presented in the simulation in Table 2). This
dataset has 1000 observations and 2000 correlated noise variables for each
predictor. While the VarDes algorithm effectively identifies true positive
variables for parameter $\nu$, it fails to select any variables for parameter
$\mu$, and only selects one false variable for $\sigma$. The absence of true
variable selections for parameter $\mu$ and $\sigma$ is attributed to small
gradients concerning this predictors, which make all variables appear
insignificant for the gradient boosting algorithm. Conversely, CF identifies
all true variables and demonstrates fewer falsely selected variables. Overall,
our correlation filtering selects 10 out of 10 true effects and 3 false
effects. Variable deselection identified 3 out of 10 true effects and selects
50 false effects.

####  4.4.1 Choice of the Correlation Threshold - Hypothesis Testing Framework

As the sample size and number of variables vary, the appropriateness of
setting $\kappa$ to an arbitrary value becomes questionable. A smaller value
might be necessary as the sample size increases, while a larger value may be
warranted with a greater number of variables. Therefore, $\kappa$ should be
adjusted accordingly to better reflect the nuances of each situation.
Therefore we propose to set $\kappa$ as the critical value in the following
hypothesis testing framework. Assume we are dealing with the variables
$(\mathbf{X}_{k})_{\cdot 1},\dots,(\mathbf{X}_{k})_{\cdot J_{k}}$ for a given
parameter $\theta_{k}$. Furthermore assume all variables and the gradient
vector $\mathbf{g}_{k}$ are standardized.

The sample version of the Pearson correlation coefficient is given by:

$\displaystyle c_{j_{k}k}=\frac{1}{n-1}({\mathbf{X}_{k}})_{\cdot
j_{k}}^{\top}\mathbf{g}_{k}=\frac{\sqrt{n}}{n-1}\frac{1}{\sqrt{n}}({\mathbf{X}_{k}})_{\cdot
j_{k}}^{\top}\mathbf{g}_{k}$

and let $C_{j_{k}k}$ denote the population version of the Pearson correlation
coefficient.
The hypothesis pair:

$\displaystyle H_{0}\text{: }R_{j^{*}_{k}k}=$
$\displaystyle\mathrm{max}\left\\{\lvert C_{1k}\rvert,\dots,\lvert
C_{J_{k}k}\rvert\right\\}=0,\text{ where we assume independence of
}({\mathbf{X}_{k}})_{\cdot j_{k}}\text{ and }\mathbf{g}_{k}.$ $\displaystyle
H_{1}\text{: }R_{j^{*}_{k}k}=$ $\displaystyle\mathrm{max}\left\\{\lvert
C_{1k}\rvert,\dots,\lvert C_{J_{k}k}\rvert\right\\}\neq 0.$

The standardization and independence under $H_{0}$ of
${(\mathbf{X}_{k}})_{\cdot j_{k}}$ and $\mathbf{g}_{k}$ implies the random
variable corresponding to the inner product ${(\mathbf{X}_{k}})_{\cdot
j_{k}}^{\top}\mathbf{g}_{k}$ is a sum of $n$ products, where each product has
mean zero and variance one, and thus
$\frac{1}{\sqrt{n}}({\mathbf{X}_{k}})_{\cdot
j_{k}}^{\top}\mathbf{g}_{k}\stackrel{{\scriptstyle
asy.}}{{\rightarrow}}\mathcal{N}\left(\mu=0,\sigma^{2}=1\right)$ according to
the central limit theorem. This implies $c_{j_{k}k}$ follows asymptotically a
normal distribution
$\mathcal{N}\left(\mu=0,\sigma^{2}=\frac{n}{(n-1)^{2}}\right)$.

We take a threshold value which corresponds to a significance level $\alpha$,
i.e., a threshold value satisfying

$\displaystyle 1-\alpha=\mathbb{P}\left(\mathrm{max}\left\\{\lvert
C_{1k}\rvert,\dots\lvert C_{J_{k}k}\rvert)\right\\}\leq\kappa|H_{0}\right).$

The independence stated in $H_{0}$ implies the pairwise independence of
$C_{1k},\dots,C_{J_{k}k}$, thus

$1-\alpha=\mathbb{P}\left(-\kappa\leq C_{1k}\leq\kappa|H_{0}\right)^{J_{k}}.$

As the asymptotic distribution for $C_{1k}$ is a normal distribution, we can
easily compute the corresponding quantiles and thus the threshold value
$\kappa$ with help of the inverse of the cumulative distribution function of
the standard normal distribution $\Phi^{-1}(p)$:

$\kappa=\Phi^{-1}\left(\frac{1+(1-\alpha)^{\frac{1}{J_{k}}}}{2}\right)\cdot\frac{\sqrt{n}}{n-1}.$

The independence assumption might appear stringent but at the later stages of
the updating, when a lot of the dependence is absorbed into the model, the
remaining dependence should be negligible and before that, the correlation is
typically high. We use $\alpha=0.05$ as the default setting. As the
simulations show, this approximation works quite well for the correlation
filtering.

####  4.4.2 Choice of the Correlation Threshold - Threshold Descent Algorithm

An alternative approach to determining the value of $\kappa$ for the CF method
involves employing a grid search framework. Threshold descent works by
starting with a high $\kappa$ value and by boosting the model until the
coefficients no longer overcome the initial threshold. This gives us the first
set of coefficients, which are then refitted in a later step. Then, the value
$\kappa$ is lowered and continuing the boosting, additional variables may
enter the model. The second set of coefficients is saved and, once again, the
threshold is lowered. This process is continued until the threshold
$\kappa=0$. These are the selection steps. For the final model, the different
sets of the selection process are then refitted until convergence. The final
model is then determined by the BIC. The threshold descent algorithm is
illustrated in Figure 4.

Figure 4:  Illustration of threshold descent on simulated NO data. Displayed
are the coefficients corresponding to parameter $\mu$. Starting with a high CF
value $\kappa=0.19$, the algorithm boosts for 200 iterations thereby selecting
a few variables. Every 200 iterations $\kappa$ gets lowered by 0.02, thereby
relaxing the correlation filtering threshold. This in turn allows for more
variables to enter the model the lower the threshold descends.

As the observation size and the number of variables grow ($n$ and $p$), the
computational burden as well as the need for fast algorithms increases.
Therefore we further propose a batched version of the CF, 4.1 and SC-BS-SDR,
similar to stochastic gradient descend.

###  4.5 Efficient Estimation for Large-Scale Data

In this subsection, we adopt a stochastic approximation approach for faster
memory efficient computation of inner products and updating steps, following
the methodology proposed by Umlauf _et al._ (2023). Specifically, we compute
the inner products
$(\mathbf{X}{k})_{[\mathbf{i}_{t}]j^{*}_{k}}^{\top}\mathbf{g}_{[\mathbf{i}_{t}]k}$
using only a randomly selected batch $\mathbf{i}_{t}\subset\\{1,\dots,n\\}$ of
size $bs$ at each iteration. These inner products are used to identify the
candidate variables for updating for each predictor.

For the updating step of 4.1 we replace the normalized derivatives with their
batchwise counterpart based only on the data $\mathbf{i}_{t}$:

$\displaystyle\partial\hat{\ell}_{j_{k}^{*}k}$
$\displaystyle=\frac{1}{bs}\frac{\partial\ell(\boldsymbol{\beta}^{[t-1]};\mathbf{y}_{[\mathbf{i}_{t}]},\mathbf{X}_{[\mathbf{i}_{t}]})}{\partial{\beta_{j_{k}^{*}k}}},$
$\displaystyle\hat{\epsilon}_{j^{*}_{k}k}$
$\displaystyle=\begin{cases}\nu\cdot\epsilon&\text{if
}|\partial\hat{\ell}_{j_{k}^{*}k}|<\nu\cdot\epsilon\text{ and }t<\rho\cdot
T\\\ |\partial\hat{\ell}_{j_{k}^{*}k}|&\text{if
}|\partial\hat{\ell}_{j_{k}^{*}k}|<\nu\cdot\epsilon\text{ and }t\geq\rho\cdot
T\\\ |\partial\hat{\ell}_{j_{k}^{*}k}|&\text{if
}\nu\cdot\epsilon\leq|\partial\hat{\ell}_{j_{k}^{*}k}|\leq\epsilon\\\
\epsilon&\text{else},\end{cases}$ (SC-BW-SDR )
$\displaystyle\boldsymbol{\beta}_{k}^{[t]}$
$\displaystyle=\boldsymbol{\beta}_{k}^{[t-1]}+\hat{\epsilon}_{j_{k}^{*}k}\cdot\mathrm{sign}\left((\mathbf{X}_{k})_{[\mathbf{i}_{t}]j_{k}^{*}}^{\top}\mathbf{g}_{[\mathbf{i}_{t}]k}\right)\cdot\mathbf{e}_{j_{k}^{*}}$
$\displaystyle j_{k}^{*}$
$\displaystyle=\underset{{j=1,\dots,J_{k}}}{\mathrm{argmax}}\,\left|({\mathbf{X}_{k}})_{[\mathbf{i}_{t}]j}^{\top}\mathbf{g}_{k}\right|.$

The batchwise variant of SC-BS-SDR is analogous to 4.5. Both batchwise
variants perform their preliminary variable selection on a subset
$\mathbf{i}_{t}$ of the data. The decision which distribution parameter or
combination of parameters to update, if any, is made on the next batch
$\mathbf{i}_{t+1}$. This leads to additional stability as the second step is
performed on data that is quasi out-of-sample. An additional advantage is also
that by randomization the risk of getting stuck in local optima can be
substantially minimized when using more complex distributions. The exact
algorithm is shown in Algorithm SBDR. Specifying a threshold of zero yields
the correlation filter-free algorithm, and replacing all possible subsets with
only the singleton subsets in the best-subset selection yields the non-cyclic
update.

For batchwise correlation filtering, a threshold derived from the assumption
$C_{j}\stackrel{{\scriptstyle
a}}{{\sim}}\mathcal{N}\left(\mu=0,\sigma^{2}=\frac{bs}{(bs-1)^{2}}\right)$ is
used.

Algorithm SBDR : Stagewise Boosting Distributional Regression.

Set $\epsilon>0$, $T>1$ and $\kappa,\rho,\nu\geq 0$, batch indices
$B=(\mathbf{i}_{1},\dots,\mathbf{i}_{T})$ of size $bs=|\mathbf{i}_{1}|$.

Initialize $\boldsymbol{\beta}^{[1]}$. Intercept-coefficients are initialized
at corresponding maximum likelihood estimates, all other coefficients at zero.
Set $\boldsymbol{\beta}^{[2]}=\boldsymbol{\beta}^{[1]}$.

for t from 2 to $T$ do

Set current batch $\mathbf{i}=\mathbf{i}_{t}$ and next batch
$\tilde{\mathbf{i}}=\mathbf{i}_{t+1}$ for $t<T$, if $t=T$, then
$\tilde{\mathbf{i}}=\mathbf{i}_{1}$.

for k from 1 to K do

$\boldsymbol{\eta}_{[\mathbf{i}]k}\leftarrow\mathbf{X}_{[\mathbf{i}]k}\boldsymbol{\beta}_{k}^{[t-1]}$
and
$\boldsymbol{\eta}_{[\tilde{\mathbf{i}}]k}\leftarrow\mathbf{X}_{[\tilde{\mathbf{i}}]k}\boldsymbol{\beta}_{k}^{[t-1]}$.

end for

Compute
$\texttt{logLik}_{[\mathbf{i}]}\leftarrow\ell(\mathbf{y}_{[\mathbf{i}]},\boldsymbol{\eta}_{[\mathbf{i}]1},\dots,\boldsymbol{\eta}_{[\mathbf{i}]K})$
and
$\texttt{logLik}_{[\tilde{\mathbf{i}}]}^{\texttt{old}}\leftarrow\ell(\mathbf{y}_{[\tilde{\mathbf{i}}]},\boldsymbol{\eta}_{[\tilde{\mathbf{i}}]1},\dots,\boldsymbol{\eta}_{[\tilde{\mathbf{i}}]K})$.

for k from 1 to K do

Update intercept $\beta_{0k}^{[t]}$ with clipped gradient update with
$\frac{1}{bs}\frac{\partial\ell(\mathbf{y}_{[\mathbf{i}]},\boldsymbol{\eta}_{[\mathbf{i}]1},\dots,\boldsymbol{\eta}_{[\mathbf{i}]K})}{\partial{\beta_{0k}}}$.

end for

for k from 1 to K do

$\mathbf{g}_{k}\leftarrow\left(\frac{\partial\log
d_{y}(y_{i},\boldsymbol{\eta}_{i1},\dots,\boldsymbol{\eta}_{iK})}{\partial{\eta_{k}}}\right)_{i\in[\mathbf{i}]}$

$j_{k}^{*}\leftarrow$ Determine non-intercept column of
$\mathbf{X}_{[\mathbf{i}]k}$ with highest correlation with $\mathbf{g}_{k}$ in
absolut value.

$c_{j^{*}_{k}k}\leftarrow$ Corresponding correlation in absolute value.

if $c_{j^{*}_{k}k}>\kappa$ then

$\partial\hat{\ell}_{j_{k}^{*}k}\leftarrow\frac{1}{bs}\frac{\partial\ell(\mathbf{y}_{[\mathbf{i}]},\boldsymbol{\eta}_{[\mathbf{i}]1},\dots,\boldsymbol{\eta}_{[\mathbf{i}]K})}{\partial{\beta_{j_{k}^{*}k}}}$

else

$\partial\hat{\ell}_{j_{k}^{*}k}\leftarrow 0$

end if

end for

Best-subset selection:

$\boldsymbol{\beta}_{1}^{\mathrm{best}},\dots,\boldsymbol{\beta}_{K}^{\mathrm{best}}\leftarrow\boldsymbol{\beta}_{1}^{[t-1]},\dots,\boldsymbol{\beta}_{K}^{[t-1]}$

$\boldsymbol{\beta}_{1}^{\mathrm{temp}},\dots,\boldsymbol{\beta}_{K}^{\mathrm{temp}}\leftarrow\boldsymbol{\beta}_{1}^{[t]},\dots,\boldsymbol{\beta}_{K}^{[t]}$

for non-empty subset $S\subset\\{1,\dots,K\\}$ do

$\nabla\hat{\mathbf{L}}_{S}\leftarrow\left(\partial\hat{\ell}_{j^{*}_{s}s}\right)^{\top}_{s\in
S}$

if $\left\|\nabla\hat{\mathbf{L}}_{S}\right\|_{2}>\epsilon$ then

Clip from above

for s in S do

$\partial\hat{\ell}_{j^{*}_{s}s}\leftarrow\frac{\partial\hat{\ell}_{j^{*}_{s}s}}{\left\|\nabla\hat{\mathbf{L}}_{S}\right\|_{2}}$

end for

end if

if $t<\rho\cdot T$ then

Clip from below

for s in S do

$\partial\hat{\ell}_{j^{*}_{s}s}=\mathrm{sign}(\partial\hat{\ell}_{j^{*}_{s}s})\cdot\mathrm{max}(\nu\cdot\epsilon,|\partial\hat{\ell}_{j^{*}_{s}s}|)$

end for

end if

for s in S do

Update coefficients

$\boldsymbol{\beta}_{s}^{[t]}\leftarrow\boldsymbol{\beta}_{s}^{[t]}+\partial\hat{\ell}_{j^{*}_{s}s}\cdot\mathbf{e}_{j_{s}^{*}}$

end for

for k in K do

$\boldsymbol{\eta}_{[\tilde{\mathbf{i}}]k}\leftarrow\mathbf{x}_{[\tilde{\mathbf{i}}]k}\boldsymbol{\beta}_{k}^{[t]}$

end for

$\texttt{logLik}_{[\tilde{\mathbf{i}}]}^{\texttt{new}}\leftarrow\ell(\mathbf{y}_{[\tilde{\mathbf{i}}]},\boldsymbol{\eta}_{[\tilde{\mathbf{i}}]1},\dots,\boldsymbol{\eta}_{[\tilde{\mathbf{i}}]K})$

if
$\texttt{logLik}_{[\tilde{\mathbf{i}}]}^{\texttt{new}}>\texttt{logLik}_{[\tilde{\mathbf{i}}]}^{\texttt{old}}$
then

If oos-likelihood improves, keep new coefficients.

$\boldsymbol{\beta}_{1}^{\mathrm{best}},\dots,\boldsymbol{\beta}_{K}^{\mathrm{best}}\leftarrow\boldsymbol{\beta}_{1}^{[t]},\dots,\boldsymbol{\beta}_{K}^{[t]}$

$\texttt{logLik}_{[\tilde{\mathbf{i}}]}^{\texttt{old}}\leftarrow\texttt{logLik}_{[\tilde{\mathbf{i}}]}^{\texttt{new}}$

end if

Reset coefficients

$\boldsymbol{\beta}_{1}^{[t]},\dots,\boldsymbol{\beta}_{K}^{[t]}\leftarrow\boldsymbol{\beta}_{1}^{\mathrm{temp}},\dots,\boldsymbol{\beta}_{K}^{\mathrm{temp}}$

end for

Keep best coefficients

$\boldsymbol{\beta}_{1}^{[t]},\dots,\boldsymbol{\beta}_{K}^{[t]}\leftarrow\boldsymbol{\beta}_{1}^{\mathrm{best}},\dots,\boldsymbol{\beta}_{K}^{\mathrm{best}}$

$\boldsymbol{\beta}_{1}^{[t+1]},\dots,\boldsymbol{\beta}_{K}^{[t+1]}\leftarrow\boldsymbol{\beta}_{1}^{\mathrm{best}},\dots,\boldsymbol{\beta}_{K}^{\mathrm{best}}$

end for

return best of
$\left(\boldsymbol{\beta}_{1}^{[t]},\dots,\boldsymbol{\beta}_{K}^{[t]}\right)_{t=1,\dots,T}$
according to the BIC.

###  4.6 Software Implementation

Our stagewise distributional regression methods are implemented in the R (R
Core Team, 2020) package stagewise (https://github.com/MattWet/stagewise). The
package supports all distribution families of the R package gamlss.dist
(Stasinopoulos and Rigby, 2022). Examples with code on how to fit the models
using the proposed Algorithm SBDR are provided in the help pages of the
package. In addition, the scripts for the following simulation study and the
lightning application are provided.

##  5 Simulation

We present a comprehensive evaluation of our proposed methods compared to
other benchmarks. Our focus is on assessing the performance using several
metrics, including the root mean squared error of the additive predictors
(RMSE), a (continuous) ranked probability score ((C)RPS; Gneiting and Raftery,
2007), the number of true positives (TP) and false positives (FP) in the
variable selection process, and the elapsed computation time. To accomplish
this, we employ three distinct distributions (NO, GA and ZANBI), each with
varying numbers of additional noise variables ($30$ and $100$), different
pairwise correlations ($0$ and $0.7$) between the variables, and various
sample sizes ($500$ up to $1000000$) to facilitate a thorough comparison of
the methods.

###  5.1 Methods

Please note that all of the boosting-type benchmark methods used in our
simulation employ non-cyclical updates (see Section 3.2). Specifically, we
consider the following benchmark methods:

  * •

True Model (TM): Here the true model is known and thus no non-informative
variables are incorporated in the model. The model is estimated with the
standard backfitting method implemented in the R package bamlss (Umlauf _et
al._ , 2021, 2022).

  * •

Gradient Boosting (GB): We use the non-cyclical gradient boosting version in
Thomas _et al._ (2018). The optimal stopping iteration (mstop) is selected by
ten-fold cross-validation. The non-cyclical gradient boosting algorithm is
implemented in the R package gamboostLSS (Hofner _et al._ , 2021).

  * •

Stability Selection (StabSel): Our third benchmark method is the non-cyclical
stability selection method (Meinshausen and Bühlmann, 2010; Thomas _et al._ ,
2018), which is also implemented in the R package bamlss. This method is based
on calculcating the variable selection frequencies on 100 resampled datasets
for the NO and GA. To make the ZANBI computational feasable, we incorporate
only 10 resampled datasets, which still takes around 250 minutes in a
challenging ZANBI setting with 10000 observations for a single replication.
For comparision a similar NO setting with 100 resampled datasets takes around
80 minutes and an increase from 10 to 100 resampled datasets increases the
computation time around 5 to 10-fold.

  * •

Variable Deselection (VarDes): First a full model via Gradient Boosting (GB
above) is estimated and then variables wich contribute less than 1% to the
risk reduction up until the mstop are deselected. After that the model with
the selected variables is refitted with mstop iterations (Strömer _et al._ ,
2022).

  * •

Semi-Analytical Adaptive Step-Length (SAASL): This method is only used for the
normal location-scale settings as it is not implemented for other
distributions (Zhang _et al._ , 2022). Early stopping is determined by ten-
fold cross-validation.

As we propose several adaptations of the general stagewise regression
algorithm – noncyclic stagewise distributional regression, best-subset
stagewise distributional regression, a new variable selection method,
correlation filtering, and a batched update expansion of them – we investigate
the performance measures described at the beginning of this section for all
possible combinations. Each method is using the BIC to determine the optimal
stopping iteration _mstop_. This is the first step, the variable selection
step. The degrees of freedom for the BIC penalty are chosen strictly as the
number of non-zero coefficients. In the batchwise variants, we use a moving
average over the computed BIC scores of the updates. In the second step, we
recompute the model using only the non-zero coefficients until convergence and
without correlation filtering. Specifically, we evaluate the performance of
the following variants of the proposed algorithm:

  * •

Standard: Non-cyclical updating with BIC. We use Algorithm SBDR with
$\kappa=0$ and the restriction that only subsets of the distributional
parameters with only a single element are allowed in the best-subset search.
This turns this algorithm effectively into the non-cyclical updating version
where the correlation filtering is turned off. Furthermore the batchsize is
set equal to the number of observations, i.e., the updating steps are selected
with full batchsize.

  * •

BS: Best-subset updating with BIC. We use Algorithm SBDR with $\kappa=0$ and
the batchsize set equal to the number of observations, i.e., the updating
steps are selected with full batchsize.

  * •

CF: Non-cyclical updating and correlation filtering with BIC. We use Algorithm
SBDR with the restriction that only subsets of the distributional parameters
with only a single element are allowed in the best-subset search. Furthermore
the batchsize is set equal to the number of observations, i.e., the updating
steps are selected with full batchsize.

  * •

BS+CF: Best-subset Updating and Correlation Filtering with BIC. We use
Algorithm SBDR and the batchsize is set equal to the number of observations,
i.e., the updating steps are selected with full batchsize.

  * •

BW: Batchwise non-cyclical updating with BIC. We use Algorithm SBDR with
$\kappa=0$ and the restriction that only subsets of the distributional
parameters with only a single element are allowed in the best-subset search.

  * •

BW+BS: Batchwise best-subset updating with BIC, i.e., Algorithm SBDR with
$\kappa=0$.

  * •

BW+CF: Batchwise non-cyclical updating and correlation filtering with BIC. We
use Algorithm SBDR with the restriction that only subsets of the
distributional parameters with only a single element are allowed in the best-
subset search.

  * •

BW+BS+CF: Batchwise best-subset updating and correlation filtering with BIC.
We use Algorithm SBDR.

  * •

ThresDesc: Best-subset updating with threshold descent for correlation
threshold selection. Final model is selected with BIC.

  * •

ThresDesc+BW: Batchwise best-subset updating with threshold descent for
correlation threshold selection. Final model is selected with BIC.

The different stagewise methods are summarized in Table 1.

Name | Update | Correlation Filtering | Batchwise
---|---|---|---
Standard | non-cyclical |  |
BS | best-subset |  |
CF | non-cyclical | ✓ |
BS + CF | best-subset | ✓ |
BW | non-cyclical |  | ✓
BW + BS | best-subset |  | ✓
BW + CF | non-cyclical | ✓ | ✓
BW + BS + CF | best-subset | ✓ | ✓
ThresDesc | best-subset | ✓* with descending |
correlation threshold |
ThresDesc + BW | best-subset | ✓* with descending | ✓
correlation threshold
Table 1:  Different methods used in the simulation.

###  5.2 Simulation Design

We simulate data from the normal distribution (NO), the Gamma distribution
(GA) and the zero-adjusted negative binomial type 1 distribution (ZANBI). All
distributions used in this paper are implemented in the R package gamlss.dist
(Stasinopoulos and Rigby, 2022). In the simulation study, we let all
distributional parameters depend on covariates which are drawn from a uniform
distribution $x_{1},\dots,x_{6}\sim\mathcal{U}(-1,1)$. The complete
specifications for each distribution are shown in Table 2.

Distribution | Parameters
---|---
$\texttt{NO}(\mu,\sigma)$ | $\mu=x_{1}+2x_{2}+0.5x_{3}-x_{4}$
$\log\left(\sigma\right)=0.5x_{3}+0.25x_{4}-0.25x_{5}-0.5x_{6}$
$\texttt{GA}(\mu,\sigma)$ | $\log\left(\mu\right)=x_{1}+2x_{3}+0.5x_{5}-x_{6}$
$\log\left(\sigma\right)=0.5x_{3}+0.75x_{4}-0.3x_{5}-0.5x_{6}$
$\texttt{ZANBI}(\mu,\sigma,\nu)$ | $\log\left(\mu\right)=0.5+0.5x_{1}-x_{3}+0.75x_{5}+0.75x_{6}$
$\log\left(\sigma\right)=-1+x_{2}-1.25x_{4}+x_{5}$
$\log\left(\frac{\nu}{1-\nu}\right)=-0.5+x_{3}-x_{4}-x_{5}$
Table 2:  Distributions and specifications for the linear predictors of the
respective parameters used in the simulation study.

Here, the mean of the GA is equal to $\mu$. The variance of the GA
distribution is given by $\mu^{2}\sigma^{2}$. For the ZANBI, $\mu$ is the mean
of the NBI part, $\mu+\sigma\mu^{2}$ is the variance of the NBI part and $\nu$
corresponds to the probability of the zero event, i.e.,
$\mathbb{P}\left(Y=0|\mu,\sigma,\nu\right)=\nu$

##### Further settings

  * •

To assess the performance of our methods across both small and large data
settings, we simulate data for varying numbers of observations including 500,
1000, 5000, 10000, 100000 and 1000000. For the largest sample sizes of 100000
and 1000000 we limit the evaluation to the true model and batchwise methods
only, as the computational burden of the other methods was deemed to be too
high.

  * •

In addition, the influence of an additional number of noise variables (denoted
with nnoise in the following) is considered with 30 and 100 noise variables.
Each predictor is modeled including all available covariates. Accordingly, for
each predictor 2+nnoise or 3+nnoise non-relevant covariates are included.

  * •

Two cases of correlation between the covariates are considered. Firstly,
covariates are uncorrelated ($\rho=0$), and secondly, random pairs of
covariates are correlated with $\rho=0.7$. The second case is generated as
follows: First, we set up a data-matrix with uniformly distributed entries
with nnobs rows and $l$ = 6+nnoise columns. Second, we multiply the data-
matrix with the cholesky factor $\mathbf{L}^{\top}$ from the cholesky
decomposition $\mathbf{L}\mathbf{L}^{\top}=\mathbf{\Sigma}$ of the covariance
matrix

$\mathbf{\Sigma}=\begin{pmatrix}1&\rho&\rho^{2}&\cdots&\rho^{l-1}\\\
\rho&1&\rho&\cdots&\rho^{l-2}\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
\rho^{l-1}&\rho^{l-2}&\rho^{l-3}&\cdots&1\end{pmatrix}$

to introduce the correlation to the data. At this point the correlation
between neighboring columns is 0.7 and the correlation between columns $i$ and
$j$, where $i,j=1,\dots,\text{6+\code{nnoise}}$, is $0.7^{\lvert i-j\rvert}$
and thus, the correlation propagates from column to column. Third, we apply a
random permutation to the columns. The columns of the final matrix correspond
to the covariates $x_{1},\dots,x_{6+\text{\code{nnoise}}}$. The final
permutation is a part of the random data generation process and therefore
varies across different replications of the simulation. For instance, two
variables may exhibit high correlation in one replication, while showing
almost zero correlation in another replication.

  * •

The used threshold values ($\kappa$) depend on the number of observations or
the batchsize used and is derived from the significance level $\alpha=0.05$.
Small sample sizes bear the risk of a small power for the correlation tests,
i.e., not rejecting $H_{0}$ although $H_{0}$ is wrong. To circumvent this we
enforce a maximum threshold $\kappa=0.175$. On the other extreme, when the
sample size is high, everything tends to be significant, thus we also enforce
a minimum threshold $\kappa=0.075$.

  * •

As mentioned before, to evaluate the performance of the methods, the (C)RPS,
and the predictor bias, i.e., the RMSE of the estimated additive predictors
and the true additive predictors
$\text{RMSE}_{k}=\sqrt{\frac{1}{n}\sum_{i}(\hat{\eta}_{ki}-\eta_{ki})^{2}}$,
for $k=\mu,\sigma,\nu$, are calculated based on an out-of-sample validation
data-set with $10000$ observations. Furthermore, the number of correctly
selected covariates and the number of falsely selected covariates in each
predictor (true and false positives, TP and FP) are counted. In addition, the
computational time is tracked for all methods.

  * •

Each combination of the simulation design is replicated 100 times.

Figure 5 illustrates the new methods, except the threshold descent method,
presented using one specific simulation setting for the NO distribution.

Figure 5:  Evolution of coefficient paths of variable selection step of new
methods illustrated with the setting NO, $n=5000$, $\rho=0.7$,
$\texttt{nnoise}=100$. Dashed vertical lines depict early stopping through
information criteria. df is the number of selected variables for the
respective distributional parameter.

It can be seen that all variants using correlation filtering have remarkably
low false positive rates. However, this approach can lead to excessive
shrinkage of the coefficients, as can be observed, for example, in the method
CF coefficient-paths for parameter $\sigma$. Therefore, as mentioned earlier,
the model is re-estimated in a second step without correlation filtering and
the optimal value of mstop is selected using the BIC.

###  5.3 Simulation Results

We focus on presenting the simulation results of the settings of the NO and
ZANBI with $\texttt{nnoise}=100$ and $\rho=0.7$ in this section. All other
settings can be found in Appendix Section 8.3. Figure 6 and Figure 7 show the
average results of 100 replications and we highlight these in the following:

####  5.3.1 False Positives and True Positives

For both NO and ZANBI, the comparison of false positive rates across various
methods indicates that correlation filtering and ThresDesc variants
consistently outperform the models Standard, BS, BW, BW+BS, GB, and SAASL in
terms of false positive rate. In smaller settings ($n\leq 1000$), the
batchwise variants (BW+CF and BW+BS+CF) tend to perform slightly worse than
the full batch variants (CF and BS+CF), though this difference diminishes for
settings with $n>1000$. Regarding true positives, CF and BS+CF perform
slightly worse in selecting all true variables in the normal distribution
setting with $n=500$ compared to other methods. However, this difference
disappears with larger sample sizes. In the ZANBI setting, best-subset methods
combined with correlation filtering (BS+CF and BW+BS+CF) achieve a marginally
better true positive rate than their non-cyclical counterparts (CF and BW+CF),
which in turn achive better true positive rates than GB and SAASL in the ZANBI
settings. The latter two methods face significant issues in the challenging
ZANBI scenario, exhibiting very low true positive rates and consequently high
predictor bias and (C)RPS values. The results are very promising, showing a
strong ability for the correlation filtering and ThresDesc variants to perform
variable selection.

####  5.3.2 (C)RPS and Predictor Bias

Regarding (C)RPS in the NO settings, all models converge to the true model
(TM) as the number of observations increases. However, due to the high false
positive rate, GB and SAASL models converge slightly slower than other
methods. The CF, BS+CF, StabSel, and ThresDesc variants all exhibit low bias
and a good (C)RPS in the NO settings. In the ZANBI settings, VarDes, GB, and
StabSel have difficulty accurately estimating some linear predictors due to
the vanishing gradient problem for both small and large $n$, resulting in high
(C)RPS values for VarDes and GB. On the other hand our proposed methods, which
all use the semi-constant step length, are able to converge to the true linear
predictors and the (C)RPS with increasing sample size.

####  5.3.3 Computation Time

The elapsed time to compute a single model is high for the methods GB, VarDes
and StabSel, as the respective models have to be computed multiple times on
different folds. This makes the use of these models in high data settings
infeasable. This is opposite to our methods which are comparable fast on
moderate sample sizes and our batchwise methods which can be computed very
efficiently even in high data settings. For example the most demanding
setting, ZANBI with $n=1000000$, takes around 145 minutes for the ThresDesc+BW
to compute. The GB with only $n=10000$ already takes around 150 minutes
despite the 100 times smaller dataset.

####  5.3.4 Summary

Correlation filtering and threshold descent are highly competitive tools for
variable selection, performing on par with the strong benchmark competitor
StabSel in the NO settings. In the highly complex ZANBI settings, all
competitors encounter some difficulties, but our proposed correlation
filtering and threshold descent methods, which use semi-constant step lengths,
are able to accurately model all linear predictors. Furthermore our methods do
not require any costly cross-validation or other subsampling approaches, which
make them computational very appealing. The batchwise variants enhance the
methods even further, making them highly scalable and fast for large datasets.
This makes it possible to handle datasets with millions of observations, where
all competitors quickly reach their computational limits.

Figure 6:  Summary of all results of the simulation study for all NO settings
with $\texttt{nnoise}=100$ and $\rho=0.7$. The average results of 100
replications vs. number of observations denoted on the abscissa is shown for
the respective metrics. Figure 7:  Summary of all results of the simulation
study for all ZANBI settings with $\texttt{nnoise}=100$ and $\rho=0.7$. The
average results of 100 replications vs. number of observations denoted on the
abscissa is shown for the respective metrics.

##  6 Application: Lightning Forecast in Austria

Lightning is a natural phenomenon that arises from imbalances in atmospheric
electrical charge during thunderstorms. The resulting electrical discharge can
cause significant damage to property and pose a danger to people and animals
(Yair, 2018). In addition to its immediate effects, lightning also has
important implications for the global climate, as it is a major source of
atmospheric nitrogen oxides (Schumann and Huntrieser, 2007), a potent
greenhouse gas that can affect climate change (Figure SPM.2 in Masson-Delmotte
_et al._ , 2021). However, the relationship between lightning and climate
change remains a subject of scientific debate (Murray, 2018), in part because
lightning processes are not yet fully understood and cannot be resolved by
weather forecast models used to predict the state of the atmosphere. As a
result, many studies of lightning rely on proxies such as cloud top height
(Price and Rind, 1992), ice flux in the mid atmosphere (Finney _et al._ ,
2014), or wind shear (Taszarek _et al._ , 2021), that are based on fairly
simple formulation and consider only certain aspects of the physical processes
involved.

To address these limitations, scholars have proposed machine learning (ML)
approaches that incorporate multiple physical processes in their analysis of
lightning (e.g., Ukkonen and Mäkelä, 2019; Simon _et al._ , 2023). These
approaches are capable of processing large amounts of data and identifying the
most relevant variables from a pool of inputs, but typically focus on
describing the occurrence of lightning via binary classification rather than
on the number of lightning counts, which is a crucial variable for
investigating flash rates (Cecil _et al._ , 2014). By using ML to gain a more
nuanced understanding of lightning processes, researchers may be better
equipped to elucidate the relationship between lightning and climate change
and to develop more accurate models for predicting lightning behavior in the
future.

The stagewise boosting algorithm we propose is ideally suited for modeling
flash rates. Firstly, it enables variable selection even with very large
datasets. Additionally, its numerical stability is crucial, particularly when
dealing with complex distributions such as the count data distribution
employed in this application.

##### Data

We analyze lightning counts using high-resolution data from the Austrian
Lightning Detection and Information System (ALDIS, Schulz _et al._ , 2005) and
reanalysis data from ERA5, the fifth generation of atmospheric reanalyses from
the European Centre for Medium-Range Weather Forecasts (ECMWF, Copernicus
Climate Change Service, 2017; Hersbach and et al., 2020). ERA5 provides hourly
globally complete and consistent pseudo-observations of the atmosphere with a
horizontal resolution of approximately $32~{}\mathrm{km}\times
32~{}\mathrm{km}$, spanning from 1950 to the present. Our study not only
provides a comprehensive description of lightning but also supports a complete
reanalysis to study climate trends in lightning (c.f. Simon _et al._ , 2023).
This is particularly important as homogeneous lightning observations from
ALDIS are available only for the period from 2010 to 2019. In total, our
dataset includes more than 9.1 million observations, which is particularly
challenging for the estimation of complex distributional regression models. We
use the data corresponding to 2010 up to 2018 ($\approx 8.2$ million
observations) as training data and the year 2019 as validation data.

##### Model Specification

Considering that only approximately $2.65\%$ of our observations exhibit
positive counts, we employ a zero-adjusted negative binomial model (ZANBI)
distribution to model lightning counts. We include $84$ physical variables
from ERA5 in our analysis, which we initially scale using a smooth kernel
density estimate of the empirical distribution function to ensure that all
variables fall within the interval $[0,1]$. Following this initial
transformation, we augment each variable with eight different transformations,
resulting in each linear predictor having a pool of $84\cdot 8=672$ variables
to select from:

$\begin{aligned} &x\mapsto x\\\ &x\mapsto x^{3}\\\ &x\mapsto\exp{(x)}\\\
&x\mapsto\mathrm{logm}(x)=\log{(1-x+0.01)}\end{aligned}\qquad\begin{aligned}
&x\mapsto x^{2}\\\ &x\mapsto\sqrt{x}\\\
&x\mapsto\mathrm{logp}(x)=\log{(x+0.01)}\\\
&x\mapsto\mathrm{logitc}(x)=\log{\left(\frac{0.999\cdot x+0.001\cdot
0.5}{1-(0.999\cdot x+0.001\cdot 0.5)}\right)}\end{aligned}$

The resulting variables are then standardized as the last variable preparation
step. To better capture the characteristics of rare positive lightning events,
we aim to improve our model by subsampling the zero count data during the
variable selection and refitting steps. We achieve this by specifying batches
consisting of $10000$ random samples each from the lightning
$\texttt{counts}=0$ and $\texttt{counts}>0$ observations. This gives us a
batch with 20000 observations for every iteration step. By carrying out the
batchwise correlation filtering with these subsampled batches, we can adjust
for the rarity of positive lightning events in our dataset. In contrast, the
approach taken by Simon _et al._ , 2019 is considerably more complex. They use
a two-stage hurdle model, consisting of a binomial model for lightning
occurrence (yes/no) and a zero-truncated negative binomial model for positive
counts. They also use stability selection with boosting to select the most
important variables for their model, which is computationally very demanding.

To account for the subsampling, we apply an intercept adjustment in the
logistic regression part of the model afterwards. This correction, derived
from King and Zeng, 2001, ensures that the corrected subsampled distribution
is consistent with the original distribution. In the Appendix Section 8.1 we
provide a proof that the parameter estimates do not require modification,
except for the logistic regression intercept of the ZANBI model
$\boldsymbol{\beta}_{0\nu}$. The subsampled parameter
$\boldsymbol{\beta}_{0\nu}^{\text{sub}}$ gets adjusted to

$\boldsymbol{\beta}_{0\nu}=\boldsymbol{\beta}_{0\nu}^{\text{sub}}-\log\left(\frac{1-\tau_{0}}{\tau_{0}}\cdot\frac{t_{0}}{1-t_{0}}\right),$

where $\tau_{0}$ is the proportion of zeros in the full dataset and $t_{0}$ is
the proportion of zeros in the subsampled dataset (i.e., in each batch). We
have $t_{0}=0.5$ and $\tau_{0}\approx 0.0265$.

We include all 672 variables as predictors and we follow a two-step procedure.
The first step involves variable selection using correlation filtering and
best-subset selection. The second step involves refitting with best-subset
selection until convergence using the variables selected in step one.

##### Results

Once variables are selected, the best-subset stagewise boosting algorithm is
applied to refit the model until convergence. The selection process is
depicted in Figure 9. It shows a balanced selection of variables for every
linear predictor. The variable importance score of the selected variables for
four years is presented in Figure 10 and for all years in the appendix Section
8.2. The algorithm selected various transformations of 22 variables into the
final model, with cswc2040, the mass of specific snow water content between
the $-20^{\circ}$C and $-40^{\circ}$C isotherms, being the most important one
according to the variable importance score (explained in caption of Figure
10). Other important variables are mcc, the medium cloud cover, cp, the
convective potential and cape, the convective available potential energy. Note
that the variable importance appears to differ for a number of variables in
the validation data, such as cape, in comparison to the scores observed in the
data used for model estimation. This also illustrates the difficulty of
developing a forecasting model for lightning counts, given the complex
physical phenomena that drive the data.

Results in Figure 8 based on the out-of-sample data show firstly, a QQ-plot
based on randomized quantile residuals (Dunn and Smyth, 1996), indicating a
quite good fit for all but some extreme values corresponding to observations
with very high counts. Secondly, the PIT histogram also indicates a overall
well calibrated model fit. The bottom row of Figure 8 shows the BIC curve of
the selection model with the minimum at iteration 1768 and a log-likelihood
plot of the selection model and the refitting model. At the transition from
the selection to the refitting model, the log-likelihood plot experience a
rapid jump, before it reaches a stable level for almost 30000 iterations.

Figure 11 exemplary shows a forecast based on the estimated model for Austria
on the afternoon (1 pm to 4 pm, local time) on July 27, 2019. Whilst the first
row shows the observed number of lightning strikes, rows 2-4 show the
probability of exceeding 0,9 and 19 strikes. It can be seen that the model is
able to properly depict areas where high lightning counts occur.

Figure 8:  Diagnostic plots are shown for the ZANBI lightning model. The top
row shows a quantil-quantil plot (top left) and a probability integral
transform plot (top right) for out-of-sample observations from the year 2019.
The BIC of the selection model is shown on the bottom left. The bottom right
shows the in-sample log-likelihood path of the selection model and the
refitting model, where the vertical red lines denote the early stopping in the
selection step induced by BIC.  Figure 9:  ZANBI Variable selection process.
Shown are the coefficient paths of the selected variables. Figure 10:  ZANBI
lightning model. Variable importance scores (VIS) - VIS are calculated using a
permutation-based approach (similar to the permutation-based VIM of Wei _et
al._ (2015)). For each variable, the log-likelihood of the model is compared
to the log-likelihood of the model where the values of that variable in the
year-specific subset of the data are permuted. The resulting scores represent
the change in log-likelihood attributable to each variable. Scores are then
normalized to the range [0, 1], with 0 indicating minimum impact and 1
indicating maximum impact on the model’s predictive performance for the
respective year. The model is trained using data from 2010 to 2018, with data
from 2019 serving as the true validation dataset.

##### Summary

Our algorithm successfully identifies relevant variables for the prediction of
lightning events in Austria. By using the stagewise boosting approach, a
complex distributional regression model could be estimated with this large
data set of over 8.2 million observations. The entire estimation process takes
about 11 1/2 hours (selection step: 4h 25min; refitting step: 7h 12min). This
processing time is considerably faster compared to other methods, for example
if stability selection would be used, and underlines the efficiency of our
algorithm. Furthermore, the numerical stability of the algorithm is crucial
for dealing with the zero-adjusted negative binomial model for flash counts,
as it ensures accurate and reliable variable selection.

Figure 11:  ZANBI lightning model. Compared is a (out of sample) forecast for
the date 2019-07-27 with the observed counts. The time evolution from 1 pm up
to 4 pm is depicted along the horizontal axis in the panel plot. The top row
shows the observed number of lightning counts, the second row shows the
probability for at least one count, the third row shows the probability for at
least 10 counts and the fourth row shows the probability for at least 20
counts.

##  7 Summary

Our contributions include the concept of stagewise boosting distributional
regression, which offers clear advantages over standard gradient boosting for
the estimation of complex models. This is achieved by using a semi-constant
step length to avoid the vanishing gradient problem. Furthermore, we introduce
a novel updating scheme, best subset updating, which iteratively updates the
best subsets of the possible distribution parameters. Additionally, we
introduce a new variable selection method, correlation filtering, which proves
to be very effective in high-dimensional data settings without the need for
costly cross-validation or other subsampling approaches, making it
computationally very attractive. Our proposed algorithm can be further
optimized applying a batchwise variant for big data problems. This variant
combines the principles of stochastic gradient descent and component- wise
boosting to compute update steps on randomly selected subsets of the data.
This approach enables the estimation of models with very large data sets and
reduces the risk of getting stuck in local optima.

In an extensive simulation study, we demonstrate that our methods, evaluated
with metrics such as false positive rates and CRPS, are at least equivalent to
other benchmark methods, but are much faster, because no complex methods for,
e.g., finding the optimal stop iteration are necessary, and also no complex
resampling. Furthermore, we show that stagewise boosting has clear advantages
in more demanding situations, e.g., when estimating a zero-adjusted negative
binomial distribution. The latter distribution was used in a challenging
lightning forecasting application, which incorporates over 9.1 million
observations with 672 possible variables to choose from for each of the three
distributional parameters of the ZANBI distribution.

One aspect not considered in this work is the inclusion of regression splines
instead of purely linear effects. In the future, we plan to integrate these
into the algorithm to enable the estimation of spatial effects, for example.

## Acknowledgements

This project was partially funded by the Austrian Science Fund (FWF) grant
number $33941$. We are grateful for data support by Gerhard Diendorfer and
Wolfgang Schulz from OVE-ALDIS. The computational results presented here have
been achieved (in part) using the LEO HPC infrastructure of the University of
Innsbruck.

## References

  * Bengio _et al._ (1994) Bengio Y, Simard P, Frasconi P (1994).  “On the Difficulty of Training Recurrent Neural Networks.”  _Technical report, Department d’Informatique et de Recherche Opérationnelle, Université de Montréal_.
  * Cecil _et al._ (2014) Cecil DJ, Buechler DE, Blakeslee RJ (2014).  “Gridded Lightning Climatology from TRMM-LIS and OTD: Dataset Description.”  _Atmospheric Research_ , 135, 404–414.  10.1016/j.atmosres.2012.06.028.
  * Copernicus Climate Change Service (2017) Copernicus Climate Change Service (2017).  “ERA5: Fifth Generation of ECMWF Atmospheric Reanalyses of the Global Climate.”  Copernicus Climate Change Service Climate Date Store (CDS).  Date of access: June 2019, https://cds.climate.copernicus.eu/cdsapp#!/home.
  * Dunn and Smyth (1996) Dunn PK, Smyth GK (1996).  “Randomized Quantile Residuals.”  _Journal of Computational and Graphical Statistics_ , 5(3), 236–244.  10.2307/1390802.
  * Finney _et al._ (2014) Finney DL, Doherty RM, Wild O, Huntrieser H, Pumphrey HC, Blyth AM (2014).  “Using Cloud Ice Flux to Parametrise Large-Scale Lightning.”  _Atmospheric Chemistry and Physics_ , 14(23), 12665–12682.  10.5194/acp-14-12665-2014.
  * Gamerman (1997) Gamerman D (1997).  “Sampling from the Posterior Distribution in Generalized Linear Mixed Models.”  _Statistics and Computing_ , 7(1), 57–68.  10.1023/a:1018509429360.
  * Gneiting and Raftery (2007) Gneiting T, Raftery AE (2007).  “Strictly Proper Scoring Rules, Prediction, and Estimation.”  _Journal of the American Statistical Association_ , 102(477), 359–378.  10.1198/016214506000001437.
  * Groll _et al._ (2019) Groll A, Hambuckers J, Kneib T, Umlauf N (2019).  “LASSO-Type Penalization in the Framework of Generalized Additive Models for Location, Scale and Shape.”  _Computational Statistics & Data Analysis_, 140, 59–74.  10.1016/j.csda.2019.06.005.
  * Hastie _et al._ (2009) Hastie TJ, Tibshirani RJ, Friedman J (2009).  _The Elements of Statistical Learning_.  2nd edition. Springer-Verlag, New York.  10.1007/978-0-387-84858-7.
  * Hersbach and et al. (2020) Hersbach H, et al (2020).  “The ERA5 Global Reanalysis.”  _Quarterly Journal of the Royal Meteorological Society_ , 146(730), 1999–2049.  10.1002/qj.3803.
  * Hofner _et al._ (2021) Hofner B, Mayr A, Fenske N, Schmid M (2021).  _gamboostLSS: Boosting Methods for GAMLSS Models_.  R package version 2.0-7, URL https://CRAN.R-project.org/package=gamboostLSS.
  * Hofner _et al._ (2016) Hofner B, Mayr A, Schmid M (2016).  “gamboostLSS: An R Package for Model Building and Variable Selection in the GAMLSS Framework.”  _Journal of Statistical Software_ , 74(1), 1–31.  10.18637/jss.v074.i01.
  * International Institute for Population Sciences and ORC Macro(2000) (IIPS) International Institute for Population Sciences (IIPS) and ORC Macro (2000).  “National Family Health Survey (NFHS-2), 1998–99 [Datasets]. IAKR42.DTA.”
  * King and Zeng (2001) King G, Zeng L (2001).  “Logistic Regression in Rare Events Data.”  _Political Analysis_ , 9, 137–163.
  * Klein _et al._ (2015a) Klein N, Kneib T, Klasen S, Lang S (2015a).  “Bayesian Structured Additive Distributional Regression for Multivariate Responses.”  _Journal of the Royal Statistical Society C_ , 64, 569–591.  10.1111/rssc.12090.
  * Klein _et al._ (2015b) Klein N, Kneib T, Lang S (2015b).  “Bayesian Generalized Additive Models for Location, Scale and Shape for Zero-Inflated and Overdispersed Count Data.”  _Journal of the American Statistical Association_ , 110(509), 405–419.  10.1080/01621459.2014.912955.
  * Masson-Delmotte _et al._ (2021) Masson-Delmotte V, Zhai P, Pirani A, Connors SL, Péan C, Berger S, Caud N, Chen Y, Goldfarb L, Gomis M, _et al._ (2021).  “Climate Change 2021: The Physical Science Basis.”  _Contribution of working group I to the sixth assessment report of the intergovernmental panel on climate change_ , pp. 3–32.
  * Meinshausen and Bühlmann (2010) Meinshausen N, Bühlmann P (2010).  “Stability Selection.”  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 72(4), 417–473.  https://doi.org/10.1111/j.1467-9868.2010.00740.x.
  * Murray (2018) Murray LT (2018).  “An Uncertain Future for Lightning.”  _Nature Climate Change_ , 8(3), 191–192.  10.1038/s41558-018-0094-0.
  * Price and Rind (1992) Price C, Rind D (1992).  “A Simple Lightning Parameterization for Calculating Global Lightning Distributions.”  _Journal of Geophysical Research: Atmospheres_ , 97(D9), 9919–9933.  10.1029/92JD00719.
  * R Core Team (2020) R Core Team (2020).  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria.  URL https://www.R-project.org/.
  * Rigby and Stasinopoulos (2005) Rigby RA, Stasinopoulos DM (2005).  “Generalized Additive Models for Location, Scale and Shape.”  _Journal of the Royal Statistical Society C_ , 54(3), 507–554.  10.1111/j.1467-9876.2005.00510.x.
  * Schulz _et al._ (2005) Schulz W, Cummins K, Diendorfer G, Dorninger M (2005).  “Cloud-to-Ground Lightning in Austria: A 10-Year Study Using Data from a Lightning Location System.”  _Journal of Geophysical Research: Atmospheres_ , 110(D9).  10.1029/2004JD005332.
  * Schumann and Huntrieser (2007) Schumann U, Huntrieser H (2007).  “The Global Lightning-Induced Nitrogen Oxides Source.”  _Atmospheric Chemistry and Physics_ , 7(14), 3823–3907.  10.5194/acp-7-3823-2007.
  * Simon _et al._ (2018) Simon T, Fabsic P, Mayr GJ, Umlauf N, Zeileis A (2018).  “Probabilistic Forecasting of Thunderstorms in the Eastern Alps.”  _Monthly Weather Review_ , 146, 2999–3009.  10.1175/MWR-D-17-0366.1.
  * Simon _et al._ (2023) Simon T, Mayr GJ, Morgenstern D, Umlauf N, Zeileis A (2023).  “Amplification of Annual and Diurnal Cycles of Alpine Lightning.”  _Climate Dynamics_.  10.1007/s00382-023-06786-8.
  * Simon _et al._ (2019) Simon T, Mayr GJ, Umlauf N, Zeileis A (2019).  “NWP-Based Lightning Prediction Using Flexible Count Data Regression.”  _Advances in Statistical Climatology, Meteorology and Oceanography_ , 5(1), 1–16.  10.5194/ascmo-5-1-2019.
  * Stasinopoulos and Rigby (2022) Stasinopoulos DM, Rigby RA (2022).  _gamlss.dist: Distributions for Generalized Additive Models for Location, Scale and Shape_.  R package version 6.1-1, URL https://CRAN.R-project.org/package=gamlss.dist.
  * Strömer _et al._ (2022) Strömer A, Staerk C, Klein N, Weinhold L, Titze S, Mayr A (2022).  “Deselection of Base-Learners for Statistical Boosting—with an Application to Distributional Regression.”  _Statistical Methods in Medical Research_ , 31(2), 207–224.  10.1177/09622802211051088.
  * Taszarek _et al._ (2021) Taszarek M, Allen JT, Brooks HE, Pilguj N, Czernecki B (2021).  “Differing Trends in United States and European Severe Thunderstorm Environments in a Warming Climate.”  _Bulletin of the American Meteorological Society_ , 102(2), 296–322.  10.1175/BAMS-D-20-0004.1.
  * Thomas _et al._ (2017) Thomas J, Hepp T, Mayr A, Bischl B (2017).  “Probing for Sparse and Fast Variable Selection with Model-Based Boosting.”  _Computational and Mathematical Methods in Medicine_ , 2017, 1421409.  ISSN 1748-670X.  10.1155/2017/1421409.
  * Thomas _et al._ (2018) Thomas J, Mayr A, Bischl B, Schmid M, Smith A, Hofner B (2018).  “Gradient Boosting for Distributional Regression: Faster Tuning and Improved Variable Selection via Noncyclical Updates.”  _Statistics and Computing_ , 28, 673–687.  /10.1007/s11222-017-9754-6.
  * Tibshirani (2015) Tibshirani RJ (2015).  “A General Framework for Fast Stagewise Algorithms.”  _Journal of Machine Learning Research_ , 16(78), 2543–2588.  URL http://jmlr.org/papers/v16/tibshirani15a.html.
  * Ukkonen and Mäkelä (2019) Ukkonen P, Mäkelä A (2019).  “Evaluation of Machine Learning Classifiers for Predicting Deep Convection.”  _Journal of Advances in Modeling Earth Systems (JAMES)_ , 11(6), 1784–1802.  10.1029/2018MS001561.
  * Umlauf _et al._ (2021) Umlauf N, Klein N, Simon T, Zeileis A (2021).  “bamlss: A Lego Toolbox for Flexible Bayesian Regression (and Beyond).”  _Journal of Statistical Software_ , 100(4), 1–53.  10.18637/jss.v100.i04.
  * Umlauf _et al._ (2018) Umlauf N, Klein N, Zeileis A (2018).  “BAMLSS: Bayesian Additive Models for Location, Scale, and Shape (and Beyond).”  _Journal of Computational and Graphical Statistics_ , 27(3), 612–627.  10.1080/10618600.2017.1407325.
  * Umlauf _et al._ (2022) Umlauf N, Klein N, Zeileis A, Köhler M (2022).  _bamlss: Bayesian Additive Models for Location Scale and Shape (and Beyond)_.  R package version 1.2-4, URL http://CRAN.R-project.org/package=bamlss.
  * Umlauf and Kneib (2018) Umlauf N, Kneib T (2018).  “A Primer on Bayesian Distributional Regression.”  _Statistical Modelling_ , 18(3-4), 219–247.  10.1177/1471082X18759140.
  * Umlauf _et al._ (2023) Umlauf N, Seiler J, Wetscher M, Simon T, Lang S, Klein N (2023).  “Scalable Estimation for Structured Additive Distributional Regression.”  10.48550/ARXIV.2301.05593.
  * Wei _et al._ (2015) Wei P, Lu Z, Song J (2015).  “Variable importance analysis: A comprehensive review.”  _Reliability Engineering & System Safety_, 142, 399–432.  ISSN 0951-8320.  https://doi.org/10.1016/j.ress.2015.05.018.
  * Yair (2018) Yair Y (2018).  “Lightning hazards to human societies in a changing climate.”  _Environmental Research Letters_ , 13(12), 123002.  10.1088/1748-9326/aaea86.
  * Zhang _et al._ (2022) Zhang B, Hepp T, Greven S, Bergherr E (2022).  “Adaptive Step-Length Selection in Gradient Boosting for Gaussian Location Scale Models.”  _Computational Statistics_ , 37, 2295–2332.  /10.1007/s00180-022-01199-3.
  * Zou _et al._ (2007) Zou H, Hastie T, Tibshirani R (2007).  “On the "Degrees of Freedom" of the Lasso.”  _The Annals of Statistics_ , 35(5), 2173 – 2192.  10.1214/009053607000000127.

##  8 Appendix

###  8.1 Proof of Consistency of Subsampling Correction and Intercept
Adjustment

King and Zeng, 2001 prove the consistency of the subsampling correction.
Following their notation, let us consider random variables $X,Y$ with density
$\mathbb{P}(X,Y)$ representing the full sample and random variables $x,y$ with
density $\mathbb{P}(x,y)$ representing the subsampled case, with all positives
and only a random selection of zeros from $X,Y$. Furthermore let $d$ and $D$
be random samples of size $n$ drawn from $\mathbb{P}(X,Y)$ and
$\mathbb{P}(x,y)$ respectively. In this setting, they showed, as
$n\rightarrow\infty$:

$\mathbb{P}(y|x,d)\cdot\underbrace{\frac{\mathbb{P}(Y|D)}{\mathbb{P}(y|d)}\cdot\frac{\mathbb{P}(x|d)}{\mathbb{P}(X|D)}}_{\text{subsampling
correction}}\stackrel{{\scriptstyle d}}{{\rightarrow}}\mathbb{P}(Y|X)$

and similar

$\mathbb{P}(Y|X,D)\stackrel{{\scriptstyle d}}{{\rightarrow}}\mathbb{P}(Y|X),$

Thus, the corrected subsampled distribution is consistend for the distribution
of interests. For the intercept adjustment in a zero adjusted count data model
the general probability mass function reads:

$\displaystyle\mathbb{P}{(y=j|x,d)}$
$\displaystyle=\begin{cases}\nu(y=0|x,d)&\text{if }j=0\\\
(1-\nu(y=0|x,d))\cdot f(y=j|x,d)&\text{else},\end{cases}$

where $\nu(y=0|x,d)=(1+\mathrm{exp}(-x_{k}\beta))^{-1}$ is the logistic
regression function and $f(\cdot)$ is a zero truncated count distribution. We
abbreviate the correction factor
$\frac{\mathbb{P}(Y=j|D)}{\mathbb{P}(y=j|d)}=\frac{\tau_{j}}{t_{j}}$ and
rewrite
$\frac{\mathbb{P}(x|d)}{\mathbb{P}(X|D)}=(\sum_{j=0}^{\infty}\mathbb{P}(y=j|x,d)\cdot\frac{\tau_{j}}{t_{j}})^{-1}$.
Furthermore, assume we only subsample the zeros, then
$\frac{\tau_{j}}{t_{j}}=\frac{1-\tau_{0}}{1-t_{0}}$ for $j>0$. The corrected
probabilities are then:

$\displaystyle\mathbb{P}(Y=0|X,D)=\mathbb{P}(y=0|x,d)\cdot\underbrace{\frac{\tau_{0}}{t_{0}}\cdot(\sum_{j=0}^{\infty}\mathbb{P}(y=j|x,d)\cdot\frac{\tau_{j}}{t_{j}})^{-1}}_{\text{subsampling
correction}}$ $\displaystyle=$
$\displaystyle\frac{\mathbb{P}(y=0|x,d)\cdot\frac{\tau_{0}}{t_{0}}}{\mathbb{P}(y=0|x,d)\cdot\frac{\tau_{0}}{t_{0}}+(1-\mathbb{P}(y=0|x,d))\cdot\frac{1-\tau_{0}}{1-t_{0}}}$
$\displaystyle=$
$\displaystyle\left(1+\frac{t_{0}}{\tau_{0}}\cdot\frac{1-\tau_{0}}{1-t_{0}}\cdot\frac{1}{\mathbb{P}(y=0|x,d)}-1\right)^{-1}$
$\displaystyle=$
$\displaystyle\left(1+\mathrm{exp}\left(-x_{k}\beta+\log\left(\frac{t_{0}}{\tau_{0}}\cdot\frac{1-\tau_{0}}{1-t_{0}}\right)\right)\right)^{-1}$

and for $i>0$:

$\displaystyle\mathbb{P}(Y=i|X,D)=\mathbb{P}(y=i|x,d)\cdot\underbrace{\frac{1-\tau_{0}}{1-t_{0}}\cdot(\sum_{j=0}^{\infty}\mathbb{P}(y=j|x,d)\cdot\frac{\tau_{j}}{t_{j}})^{-1}}_{\text{subsampling
correction}}$ $\displaystyle=$
$\displaystyle\frac{\mathbb{P}(y=i|x,d)\cdot\frac{1-\tau_{0}}{1-t_{0}}}{\mathbb{P}(y=0|x,d)\cdot\frac{\tau_{0}}{t_{0}}+(1-\mathbb{P}(y=0|x,d))\cdot\frac{1-\tau_{0}}{1-t_{0}}}$
$\displaystyle=$
$\displaystyle\left(\frac{1-\mathbb{P}(y=0|x,d)}{\mathbb{P}(y=i|x,d)}+\frac{\tau_{0}}{1-\tau_{0}}\cdot\frac{1-t_{0}}{t_{0}}\cdot\frac{\mathbb{P}(y=0|x,d)}{\mathbb{P}(y=i|x,d)}\right)^{-1}$
$\displaystyle=$ $\displaystyle
f(y=i|x,d)\left(1+\frac{\tau_{0}}{1-\tau_{0}}\cdot\frac{1-t_{0}}{t_{0}}\cdot\frac{\mathbb{P}(y=0|x,d)}{(1-\mathbb{P}(y=0|x,d))}\right)^{-1}$
$\displaystyle=$
$\displaystyle\left(1-\left(1+\mathrm{exp}\left(-x_{k}\beta+\log\left(\frac{t_{0}}{\tau_{0}}\cdot\frac{1-\tau_{0}}{1-t_{0}}\right)\right)\right)^{-1}\right)\cdot
f(y=i|x,d).$

In both cases only a constant term was added to the linear predictor of the
logistic regression part and the estimates of the non intercept variables are
not changed. This is in effect an intercept adjustment.

###  8.2 Flash Model - VIS all years

Figure 12:  ZANBI lightning model. Variable importance scores (VIS) - VIS are
calculated using a permutation-based approach. For each variable, the log-
likelihood of the model is compared to the log-likelihood of the model where
the values of that variable in the year-specific subset of the data were
permuted. The resulting scores represent the change in log-likelihood
attributable to each variable. Scores are then normalized to the range [0, 1],
with 0 indicating minimum impact and 1 indicating maximum impact on the
model’s predictive performance for the respective year. The model was trained
using data from 2010 to 2018, with data from 2019 serving as the true
validation dataset.

###  8.3 Simulation Results

Figure 13:  Simulation study. Average (C)RPS rates of 100 replications vs.
number of observations for the different settings. Figure 14:  Simulation
study. Average (C)RPS rates of 100 replications vs. number of observations for
the different settings. Figure 15:  Simulation study. Average RMSE of the
$\mu$-predictor of 100 replications vs. number of observations for the
different settings. Figure 16:  Simulation study. Average RMSE of the
$\mu$-predictor of 100 replications vs. number of observations for the
different settings. Figure 17:  Simulation study. Average RMSE of the
$\sigma$-predictor of 100 replications vs. number of observations for the
different settings. Figure 18:  Simulation study. Average RMSE of the
$\sigma$-predictor of 100 replications vs. number of observations for the
different settings. Figure 19:  Simulation study. Average RMSE of the
$\nu$-predictor of 100 replications vs. number of observations for the
different settings. Figure 20:  Simulation study. Average true positive rate
of 100 replications vs. number of observations for the different settings.
Figure 21:  Simulation study. Average true positive rate of 100 replications
vs. number of observations for the different settings. Figure 22:  Simulation
study. Average false positive rates of $100$ replications vs. number of
observations for the different settings. Figure 23:  Simulation study. Average
false positive rates of $100$ replications vs. number of observations for the
different settings. Figure 24:  Simulation study. Average time needed to
compute a single model of $100$ replications vs. number of observations for
the different settings. Figure 25:  Simulation study. Average time needed to
compute a single model of $100$ replications vs. number of observations for
the different settings.