Assumptions

We assume that our training set is given by a data matrix 𝒟^train ∈ ℝ^{n × (p + 1)} where n is the number of observations and p is the number of predictors, i.e., our data matrix can be written as 𝒟^train = (X^train, Y^train) for the regressor matrix X^train ∈ ℝ^n × p and the response vector Y^train ∈ ℝⁿ.

The relation between the response and the regressors can be described by some map f : ℝ^p ↦ ℝ in the sense that y = f(x) + ϵ for some error term ϵ with mean zero and variance σ² ∈ ]0, ∞[.

The goal is to approximate f by a sparse and stable model where ‘’stable’’ refers to the property that the set of selected predictors does only change insignificantly on similar data.

Boosting families

Our package uses the family objects from the package mboost to define new loss functions. We refer to mboost (Hothorn et al. (2017)) for more details. For our applications, we allow for non-differentiable loss functions like ranking losses, i.e., losses that do not provide a gradient, leaving the ngradient argument of the family objects empty. We provide some specific losses for ranking problems. As a simple example, we demonstrate the family Rank():

y<-c(-3, 10.3,-8, 12, 14,-0.5, 29,-1.1,-5.7, 119)
yhat<-c(0.02, 0.6, 0.1, 0.47, 0.82, 0.04, 0.77, 0.09, 0.01, 0.79)
Rank()@risk(y,yhat)
#> [1] 0.1777778

The Rank family implements the hard ranking loss (S. Clémençon, Lugosi, and Vayatis (2008)) which compares the ordering of the two inserted vectors. More precisely, it sums up all misrankings, i.e., all pairs (i, j) where Y_i is higher/lower than Y_j but where Ŷ_i is lower/higher than Ŷ_j. At the end, this sum is standardized such that a hard ranking loss of 0 indicates that there are no misrankings while a hard ranking loss of, say, 0.5, indicates that in the half of all such possible pairs a misranking occurs. We provide further families for the so-called localized and weak ranking loss (Stéphan Clémençon and Vayatis (2007)), including a normalized variant for the latter.

For more details about ranking problems, see S. Clémençon, Lugosi, and Vayatis (2008) or Werner (2019a) as well as references therein.

SingBoost

The main function of this package is the singboost function which is based on glmboost (Hothorn et al. (2017)). In this version of the package, we restricted ourselves to linear regression base learners.

The usage of singboost is very similar to that of glmboost. The main difference is that SingBoost includes so-called ‘’singular iterations’’. In contrast to the gradient descent step executed in standard Boosting algorithms where the baselearner that fits the current negative gradient vector best, the singular iterations correspond to a secant step, i.e., the least-squares baselearner that fits the current residuals, evaluated in the (possibly complicated) loss function, best is taken. Therefore, singboost allows for loss functions that do not provide a (unique) gradient, so even a family object that has an empty ngradient argument can be inserted as the family argument of the singboost function.

The motivation behind these singular iterations comes from a measure-theoretical perspective in the sense that the selection frequencies can be regarded as a ‘’column measure’’ (cf. Werner and Ruckdeschel (2019)). Clearly, different loss functions lead to different column measures and potentially to different sets of selected variables, i.e., the columns that are relevant for our target loss function, maybe a ranking loss, and which are not selected by some L−Boosting for another loss function L, may be identified with a ‘’singular part’’ of column measures (cf. Werner and Ruckdeschel (2019)). We always use L₂−Boosting (Bühlmann and Yu (2003), Bühlmann (2006)) as reference Boosting model. Nevertheless, we assume that L₂−Boosting already selects relevant variables which would make a Boosting algorithm only containing singular iterations and therefore being computationally quite expensive unnecessary. Therefore, SingBoost alternates between singular iterations and standard L₂−Boosting iterations where the frequency of the singular iterations can be set by the parameter M, i.e., every M−th iteration is a singular iteration. For the sake of a wide applicability, we do not exclude standard (i.e., differentiable) loss functions. Note that our implementation contains an additional LS argument indicating whether singular iterations are gradient- or secant-based, requiring to be set to TRUE if a non-differentiable loss function is used.

The further input arguments are the number of iterations m_iter (miter), the learning rate κ (kap) and the data set 𝒟 (D). It is important to note that the inserted data set must not contain an intercept column since this column will be added automatically internally in singboost. If one does not want an intercept coefficient, one has to center the data beforehand. Furthermore, the last column must contain the response variable, so singboost always treats the last column as the response column and all other columns as regressor columns. Since a formula argument cannot be inserted, one has to generate the model matrix first if basic functions, interaction terms or categorical variables are included. For the particular case of localized ranking losses, the argument best controls how large the proportion of the best instances is.

Consider the following simple example:

glmres<-glmboost(Sepal.Length~.,iris)
glmres
#> 
#>   Generalized Linear Models Fitted via Gradient Boosting
#> 
#> Call:
#> glmboost.formula(formula = Sepal.Length ~ ., data = iris)
#> 
#> 
#>   Squared Error (Regression) 
#> 
#> Loss function: (y - f)^2 
#>  
#> 
#> Number of boosting iterations: mstop = 100 
#> Step size:  0.1 
#> Offset:  5.843333 
#> 
#> Coefficients: 
#>      (Intercept)      Sepal.Width     Petal.Length Speciesvirginica 
#>      -3.36560001       0.53639026       0.46090783      -0.01924627 
#> attr(,"offset")
#> [1] 5.843333
attributes(varimp(glmres))$self
#> [1] 0.00 0.38 0.57 0.00 0.00 0.05
attributes(varimp(glmres))$var
#> [1] (Intercept)       Sepal.Width       Petal.Length      Petal.Width      
#> [5] Speciesversicolor Speciesvirginica 
#> 6 Levels: (Intercept) < Petal.Width < ... < Petal.Length
firis<-as.formula(Sepal.Length~.)
Xiris<-model.matrix(firis,iris)
Diris<-data.frame(Xiris[,-1],iris$Sepal.Length)
colnames(Diris)[6]<-"Y"
coef(glmboost(Xiris,iris$Sepal.Length))
#>      (Intercept)      Sepal.Width     Petal.Length Speciesvirginica 
#>      -3.36560001       0.53639026       0.46090783      -0.01924627 
#> attr(,"offset")
#> [1] 5.843333
singboost(Diris)
#> $`Selected variables`
#> [1] "Petal.Length>=Sepal.Width>=Speciesvirginica"
#> 
#> $Coefficients
#> [1]  2.47773332  0.53639026  0.46090783  0.00000000  0.00000000 -0.01924627
#> 
#> $Freqs
#> [1] 0.00 0.38 0.57 0.00 0.00 0.05
#> 
#> $VarCoef
#>        Intercept      Sepal.Width     Petal.Length Speciesvirginica 
#>       2.47773332       0.53639026       0.46090783      -0.01924627
singboost(Diris,LS=TRUE)
#> $`Selected variables`
#> [1] "Petal.Length>=Sepal.Width>=Speciesvirginica"
#> 
#> $Coefficients
#> [1]  2.47773332  0.53639026  0.46090783  0.00000000  0.00000000 -0.01924627
#> 
#> $Freqs
#> [1] 0.00 0.38 0.57 0.00 0.00 0.05
#> 
#> $VarCoef
#>        Intercept      Sepal.Width     Petal.Length Speciesvirginica 
#>       2.47773332       0.53639026       0.46090783      -0.01924627

The application to the iris data set should demonstrate that L₂−Boosting is a special instance of our SingBoost algorithm (Werner and Ruckdeschel (2019)). More precisely, using the squared loss function, our SingBoost algorithm results in exactly the same model and coefficients, provided that the hyperparameters are identical. We did not specify them explicitly in the example, but the number of iterations in singboost is m_iter = 100 per default and the learning rate is κ = 0.1, mimicking the default setting for these parameters for L₂−Boosting in mboost (Hothorn et al. (2017)). Note that we did not explicitly insert a Boosting family. Per default, singboost uses Gaussian() which corresponds to standard L₂−Boosting. The LS argument indicates that we indeed perform so-called singular iterations. Clearly, for the squared loss function, both results are the same since L₂−Boosting with least-squares baselearners implicitly also takes the best baselearner in the sense that it fits the current residuals best, evaluated in the squared loss (and which is implemented however much more sophisticatedly in mboost by using correlation updates), see e.g. Zhao and Yu (2004), Zhao and Yu (2007). For clarity, we manually renamed the response variable Sepal.Length as Y in advance in the previous example (which is not necessary). singboost reports the selected variables, the coefficients and the selection frequencies, i.e., the relative number of iterations where a particular predictor has been selected. Note that singboost does not report the offset and the intercept separately.

singboost of course also allows for categorical variables, interaction terms etc.. See the following example:

glmres2<-glmboost(Sepal.Length~Petal.Length+Sepal.Width:Species,iris)
finter<-as.formula(Sepal.Length~Petal.Length+Sepal.Width:Species-1)
Xinter<-model.matrix(finter,iris)
Dinter<-data.frame(Xinter,iris$Sepal.Length)
singboost(Dinter)
#> $`Selected variables`
#> [1] "Petal.Length>=Sepal.Width.Speciessetosa>=Sepal.Width.Speciesvirginica"
#> 
#> $Coefficients
#> [1] 4.00831289 0.45662497 0.08896750 0.00000000 0.01751541
#> 
#> $Freqs
#> [1] 0.00 0.60 0.36 0.00 0.04
#> 
#> $VarCoef
#>                    Intercept                 Petal.Length 
#>                   4.00831289                   0.45662497 
#>    Sepal.Width.Speciessetosa Sepal.Width.Speciesvirginica 
#>                   0.08896750                   0.01751541
coef(glmres2)
#>                  (Intercept)                 Petal.Length 
#>                  -1.83502045                   0.45662497 
#>    Sepal.Width:Speciessetosa Sepal.Width:Speciesvirginica 
#>                   0.08896750                   0.01751541 
#> attr(,"offset")
#> [1] 5.843333

The example above shows how to deal with interaction terms. In fact, instead of inserting a formula argument as input for singboost, build the model matrix (without the intercept column) and enter this matrix as input argument. The same procedure is valid for categorical variables or basis functions. Note that we did not yet implement a group structure, i.e., we always treat all columns separately and do not perform block-wise updates as suggested in Gertheiss et al. (2011).

We now demonstrate how to apply SingBoost for quantile regression and hard ranking regression:

glmquant<-glmboost(Sepal.Length~.,iris,family=QuantReg(tau=0.75))
coef(glmquant)
#>  (Intercept)  Sepal.Width Petal.Length 
#>   -3.1274211    0.5191849    0.4756341 
#> attr(,"offset")
#> 50% 
#> 5.8
attributes(varimp(glmquant))$self
#> [1] 0.24 0.27 0.49 0.00 0.00 0.00
singboost(Diris,singfamily=QuantReg(tau=0.75),LS=TRUE)
#> $`Selected variables`
#> [1] "Petal.Length>=Sepal.Width>=Speciesvirginica"
#> 
#> $Coefficients
#> [1]  2.47836568  0.53612658  0.46095525  0.00000000  0.00000000 -0.01925953
#> 
#> $Freqs
#> [1] 0.00 0.38 0.57 0.00 0.00 0.05
#> 
#> $VarCoef
#>        Intercept      Sepal.Width     Petal.Length Speciesvirginica 
#>       2.47836568       0.53612658       0.46095525      -0.01925953
singboost(Diris,singfamily=QuantReg(tau=0.75),LS=TRUE,M=2)
#> $`Selected variables`
#> [1] "Petal.Length>=Sepal.Width>=Speciesvirginica"
#> 
#> $Coefficients
#> [1]  2.44900980  0.54702733  0.45925162  0.00000000  0.00000000 -0.01196688
#> 
#> $Freqs
#> [1] 0.00 0.42 0.55 0.00 0.00 0.03
#> 
#> $VarCoef
#>        Intercept      Sepal.Width     Petal.Length Speciesvirginica 
#>       2.44900980       0.54702733       0.45925162      -0.01196688
singboost(Diris,singfamily=Rank(),LS=TRUE)
#> $`Selected variables`
#> [1] "Petal.Length>=Sepal.Width>=Speciesvirginica>=Speciesversicolor"
#> 
#> $Coefficients
#> [1]  2.485989684  0.534163798  0.460436981  0.000000000  0.000344599
#> [6] -0.018630537
#> 
#> $Freqs
#> [1] 0.00 0.38 0.56 0.00 0.01 0.05
#> 
#> $VarCoef
#>         Intercept       Sepal.Width      Petal.Length Speciesversicolor 
#>       2.485989684       0.534163798       0.460436981       0.000344599 
#>  Speciesvirginica 
#>      -0.018630537

The last simulation has the argument M=2 which indicates that we alternate between singular and L₂−Boosting iterations where the second simulation uses the default M=10.

More details concerning the implementation of the SingBoost algorithm can be found in Werner (2019b).

Coefficient paths

If coefficient paths have to be drawn, use the path.singboost function followed by the singboost.plot function instead of the singboost function since only path.singboost also saves the intercept and coefficients paths. In the singboost.plot function, M and m_iter have to be inserted again:

singpath<-path.singboost(Diris)
singboost.plot(singpath,10,100)

Column Measure Boosting (CMB)

CMB intends to aggregate SingBoost models that have been fitted on subsamples of the training data. The main motivation behind this idea is that SingBoost is expected to detect relevant variables for the respective (complicated) loss function that L₂−Boosting would not select. From a measure-theoretical perspective, these columns may be identified with a ‘’singular part’’ (cf. Werner and Ruckdeschel (2019)). In order to propose the opportunity to study this singular part for specific loss functions, we provide the CMB procedure that aggregates SingBoost models in order to get a more ‘’stable’’ singular part. Note that the CMB procedure does not replace a Stability Selection. Furthermore, if one is only interested in getting a suitable stable model, we recommend to only use a Stability Selection without the CMB procedure by setting Bsing=1 and nsing=ncmb when using the CMB3S function (see below) for the sake of computational costs.

The main idea is to draw B^sing (Bsing) subsamples with n_sing (nsing) rows from the training set and to compute the singboost model on each. Then, the models are aggregated where different aggregation procedures can be used, i.e., either a simple average (weights1) or a weighted average based on the reciprocal of the out-of-sample loss (weights2) computed on the complement of the current subsample.

Let us however provide an example for CMB:

set.seed(19931023)
cmb1<-CMB(Diris,nsing=100,Bsing=50,alpha=0.8,singfam=Rank(),
evalfam=Rank(),sing=TRUE,M=10,m_iter=100,
kap=0.1,LS=TRUE,wagg='weights1',robagg=FALSE,lower=0)
cmb1
#> $`Column measure`
#> [1] 0.000 1.000 1.000 0.325 0.525 0.800
#> 
#> $`Selected variables`
#> [1] "Sepal.Width>=Petal.Length>=Speciesvirginica>=Speciesversicolor>=Petal.Width"
#> 
#> $`Variable names`
#> [1] "Intercept"         "Sepal.Width"       "Petal.Length"     
#> [4] "Petal.Width"       "Speciesversicolor" "Speciesvirginica" 
#> [7] "Y"                
#> 
#> $`Row measure`
#>   [1] 0.725 0.650 0.600 0.650 0.750 0.650 0.775 0.700 0.600 0.500 0.725 0.725
#>  [13] 0.675 0.650 0.675 0.700 0.725 0.600 0.575 0.600 0.625 0.750 0.650 0.675
#>  [25] 0.650 0.650 0.750 0.625 0.625 0.650 0.675 0.700 0.700 0.725 0.600 0.600
#>  [37] 0.550 0.675 0.625 0.725 0.700 0.675 0.625 0.600 0.650 0.650 0.650 0.675
#>  [49] 0.825 0.525 0.525 0.750 0.625 0.600 0.625 0.750 0.600 0.700 0.700 0.625
#>  [61] 0.575 0.675 0.625 0.475 0.775 0.600 0.725 0.700 0.675 0.550 0.575 0.700
#>  [73] 0.600 0.700 0.775 0.625 0.700 0.700 0.525 0.700 0.800 0.700 0.650 0.625
#>  [85] 0.650 0.825 0.675 0.725 0.675 0.675 0.800 0.775 0.675 0.650 0.700 0.725
#>  [97] 0.400 0.700 0.775 0.775 0.750 0.625 0.700 0.650 0.650 0.675 0.725 0.725
#> [109] 0.700 0.700 0.625 0.650 0.675 0.600 0.650 0.650 0.525 0.650 0.650 0.650
#> [121] 0.625 0.700 0.600 0.625 0.575 0.675 0.750 0.700 0.700 0.575 0.700 0.650
#> [133] 0.725 0.725 0.700 0.725 0.775 0.700 0.725 0.650 0.725 0.700 0.750 0.525
#> [145] 0.650 0.625 0.725 0.650 0.600 0.725

In this example, we apply SingBoost with the hard ranking loss (inserted through singfam=Rank()) to B^sing = 50 subsamples with n_sing = 100 rows from the iris data set each. After computing the SingBoost models, their out-of-sample performance is evaluated based on the complement of the subsample used for training (and thus always contains 50 observations) using again the hard ranking loss (inserted through evalfam=Rank()). The parameter alpha=0.8 indicates that we only use the best 80% (i.e., the best 40) of the SingBoost models for aggregation and wagg=‘weights1’ leads to a simple average of the selection frequencies. Finally, the aggregated selection frequencies (which form an empirical aggregated ‘’column measure’‘) are reported as well as the selected variables. The last part of the output, the empirical aggregated’‘row measure’‘, reports the’‘selection frequencies’’ for the rows. This is not trivial since these frequencies are not just the relative numbers of occurences of the respective rows in the training sets due to the partitioning but also contain information about the ‘’quality’’ of the instances since the relative numbers of occurences in the training sets corresponding to the best models are computed.