Automated Model Selection for Generalized Linear Models
Benjamin Schwendinger Florian Schwendinger Laura Vana-Gür












ict]Benjamin Schwendingercor1 benjaminschwe@gmail.com
[cor1]Corresponding author.
aau]Florian Schwendinger
cstat]Laura Vana-Gür
[ict]organization=Institute of Computer Technology, TU Wien,
            addressline=Gußhausstraße 27-29,
            city=Vienna,
            postcode=1040,
            country=Austria
[aau]organization=Department of Statistics, University of Klagenfurt,
            addressline=Universitätsstraße 65-57,
            city=Klagenfurt,
            postcode=9020,
            country=Austria
[cstat]organization=Institute of Statistics and Mathematical Methods in Economics, TU Wien,
            addressline=Wiedner Hauptstraße 7,
            city=Vienna,
            postcode=1040,
            country=Austria


In this paper, we show how mixed-integer conic optimization can be used to combine feature subset selection with holistic generalized linear models to fully automate the model selection process.
Concretely, we directly optimize for the Akaike and Bayesian information criteria while imposing constraints designed to deal with multicollinearity in the feature selection task.
Specifically, we propose a novel pairwise correlation constraint that combines the sign coherence constraint with ideas from classical statistical models like Ridge regression and the OSCAR model.







Feature Subset Selection Holistic Generalized Linear Models Mixed-Integer Conic Optimization Best Subset Selection Computational Statistics Generalized Linear Models













§ INTRODUCTION

Model selection is an important but typically time-consuming task.
Several authors suggest using mixed-integer optimization to automate model selection.
The classical best subset selection (BSS) problem <cit.> identifies the best k features out of p possible features according to some goodness-of-fit (GOF) measure.
For linear regression, the BSS problem, in the least-squares sense is non-convex and known to be NP-hard <cit.>.
In practice, the BSS problem is commonly approached using methods such as stepwise regression or exhaustive enumeration <cit.>.
<cit.> highlight the significant advancements in mixed-integer optimization solvers since the 1970s when BSS problems were first investigated.
As a result, the BSS problem, which was once perceived as intractable, can now be solved by mixed-integer linear optimization solvers and used to find the solution of real-world linear regression problems.
In the case of logistic regression, the BSS problem becomes a mixed-integer non-linear optimization problem.
To address this, <cit.> explore the utilization of several general-purpose non-linear mixed-integer solvers.
Another approach to tackle the BSS problem in linear and logistic regression is presented by <cit.>, who employ a hybrid approach combining cyclic coordinate descent with a local combinatorial search.

An extension to the BSS problem is the integration of an information criterion (IC) such as the Akaike information criterion (AIC) or the Bayesian information criterion (BIC) into the optimization problem's objective. 
This enables the selection of not only the best features for a fixed k but the overall optimal features according to the chosen goodness-of-fit measure.
This approach is commonly referred to in the optimization literature as feature subset selection (FSS).
Based on the division of subset selection algorithms into classes of filter, wrapper and embedded methods as done by <cit.>, FSS classifies as an embedded method. 
FSS has been applied for various model classes including generalized linear models (GLMs).
<cit.> apply FSS to linear regression using Mallows' C_p statistic as a GOF measure.
By using intelligent branch-and-bound strategies, <cit.> demonstrate the ability to solve the FSS problem for linear regression with 1,000s of observations and 100s of features within seconds. 
<cit.> propose a linear approximation approach to solve the FSS problem in logistic regression models.
They utilize a tangent line approximation of the logistic loss function, resulting in a mixed-integer linear problem. 
Likewise, <cit.> suggest using a piecewise-linear approximation, solvable through mixed-integer linear programming, to address the FSS problem in Poisson models.
Many heuristic approaches for solving the FSS problem for GLMs are based on known metaheuristics.
These include simulated annealing <cit.>, genetic algorithms (), ant colony optimization <cit.> and particle swarm optimization <cit.>.

In this paper, we aim to provide a unified framework for performing automated model selection in GLMs by solving the FSS problem using mixed-integer conic programming and to extend existing approaches to FSS to effectively handle two common challenges in GLM optimization: separation and multicollinearity.
Furthermore, our paper is the first to propose utilizing conic optimization specifically for the FSS problem in Poisson regression models.

The proposed framework builds on the class of holistic generalized linear models (HGLMs) introduced by <cit.> in  package holiglm for R.
HGLMs extend linear, binomial and Poison GLMs by adding constraints designed to improve the model quality (e.g., restricting the number of variables entering the model, enforcing coherent signs of coefficients, etc.).
Instead of approximating the log-likelihood or its components, holiglm formulates the underlying optimization problems as conic optimization problems, providing a more reliable <cit.> and accurate solution approach.
The parameters of these constrained models are obtained by using 
(mixed-integer) conic optimization.

While the framework in <cit.> allows for the specification of an upper bound on the number of variables to be selected in the model, this bound must still be predetermined prior to solving the underlying optimization problem.
In particular, we modify the framework of HGLMs for the purpose of FSS and directly integrate the AIC and BIC into the objective function rather than using the likelihood function as an objective. 
This allows us to obtain an exact solution to the FSS problem without resorting to piecewise-linear approximation methods while treating FSS for linear, binomial and Poison GLMs in a unified way.

As mentioned above, a primary focus of the paper is addressing the issue of separation and multicollinearity in the automated model selection process.
For tackling strong multicollinearity, we propose a novel pairwise correlation constraint that combines ideas from the sign coherence constraint <cit.> with a simultaneous restriction of equal-magnitude coefficients.
The idea of restricting coefficients to have the same magnitude can also be found in other statistical models, such as Ridge regression <cit.>, where the coefficients shrink equally or the OSCAR (Octagonal Shrinkage and Clustering Algorithm for Regression) model <cit.>, where exact equality of clustered coefficients is required.
Separation on the other hand is characterized by extreme overlap or distinct separation in the data and can result in unbounded optimization problems that common solvers often fail to detect.
To ensure reliable results, it is therefore vital to verify the existence of solutions.
For the binomial family with logit, probit, log and complementary log-log link, this can be done with a linear program <cit.> which we employ to avoid wrongly reporting solutions for problems whose solution is actually not determinable.

Through extensive simulation studies, we demonstrate the feasibility and practicality of our approach.
The results clearly indicate that our proposed constraint outperforms existing methods, offering a more accurate and efficient model selection process for GLMs, effectively addressing the challenges posed by multicollinearity.







The remainder of this paper is structured as follows: 
Section <ref> introduces best subset selection, feature subset selection and holistic generalized linear models.
Section <ref> explores common pitfalls that arise when estimating certain GLMs, such as failure to converge or the nonexistence of a solution and possible solutions to mitigate these problems.
In Section <ref>, we present our proposed optimization problem for automated model selection in GLMs. The importance of dealing with multicollinearity when aiming for automated feature selection is illustrated through a simulation study in Section <ref>.
Section <ref> concludes the paper.



§ FEATURE SUBSET SELECTION IN GLMS



In this section we set the stage for introducing the proposed modeling approach 
for performing FSS in GLMs. More specifically, we start with a brief introduction
to GLMs in Section <ref> and show how to formulate
the likelihood optimization problem such that it can be solved using
(mixed-integer) conic programming.
FSS is an extension to BSS, in the sense that the objective (i.e., log-likelihood) 
in the BSS problem is replaced by an information criterion. Before introducing
the FSS problem, we introduce the formulation of the BSS problem as a conic 
program in Section <ref>. 
Section <ref> introduces the information criteria we use to extend BSS to FSS.


 §.§ Generalized linear models and conic optimization

In this work we focus on solving the FSS problem for GLMs with 
linear, binomial and Poisson families. Generally, GLMs as introduced by <cit.>, are a classs of models with probability density functions that belong to the exponential dispersion model (EDM) family with probability density function:

    f(y; θ, ϕ) = exp((yθ - b(θ))/ϕ + c(y, ϕ)).

Here, b(·) and c(·) are well-defined functions that vary depending on the specific distribution.
In addition, in the presence of a (design) matrix of covariates X with p+1 columns
(including a column of ones), a GLM has a linear predictor η = Xβ, and a link function g that establishes the relationship between the linear combination of the p+1 covariates and the mean of response y_i: g(E(y_i))=η_i. 
Given g and b, θ is then a function of η and therefore of β. 
Given a sample of n independent and identically distributed response observations
y^⊤ = (y_1, …, y_n) and observed covariates, the estimation of the parameters is usually done by maximum likelihood,
and the maximum likelihood estimate (MLE) of (β, ϕ) are the values (β^*, ϕ^*) that maximize the (log)-likelihood function for the EDM family:

    logℒ(β; y) = ∑_i=1^n log f(y; θ, ϕ) = ∑_i=1^n y_i θ_i(β) - b(θ_i(β))/ϕ_i + c(y_i, ϕ_i).



Conic optimization provides a framework for expressing the maximization of the 
log-likelihood function of various GLMs as convex optimization problems. 
A conic optimization problem is designed to model convex problems by
optimizing a linear objective function over the intersection of an affine hyperplane and a
nonempty closed convex cone.
The log-likelihood maximization can be reformulated as a conic problem.
The reason for this lies in the fact that the log-likelihood of common GLMs includes functions that can be represented by convex cones, which in turn can be solved by modern conic optimization solvers.
The estimation by means of the conic programming is in 
turn feasible given that efficient optimization solvers exist which can provide exact solutions. 
More specifically, the MLE for the linear regression model (Gaussian family with identity link) is the solution of a convex optimization problem which uses the second-order cone ^n_soc := {(t,x) ∈^n  |  x ∈^n-1, t ∈ R, ||x||_2 ≤ t}.
Since both logistic regression and Poisson regression involve exponential and logarithmic terms in their log-likelihoods, the primal exponential cone 𝒦_expp := {(x,y,z) ∈^3  |  y > 0, ye^x/y≤ z }∪{(x,0,z) ∈^3  |  x ≤ 0, z ≥ 0} is utilized to represent them.


A comprehensive introduction to conic optimization can be found in <cit.>.
For the detailed explanation and derivation of conic formulations for various GLMs based on the family and link information, we refer to the appendix provided by <cit.>.



 §.§ Best subset selection

We introduce the classical best subset selection (BSS) problem before presenting its 
extension to FSS.
The classical best subset selection (BSS) problem is concerned with determining the best k features from a set of p potential features using some goodness-of-fit (GOF) metric.
However, the BSS problem is non-convex and NP-hard <cit.>. 
Surrogate models are often used to overcome this computational burden, such as those that incorporate an L_1 penalty or a combination of L_1 and L_2 penalties <cit.>.
Adding an L_1 penalty gives the least absolute shrinkage and selection operator (LASSO) <cit.>.
However, the performance of BSS and LASSO depends on the signal-to-noise ratio.
A comprehensive overview by <cit.> highlights situations for the linear regression case where BSS outperforms LASSO and vice versa.
<cit.> approximate the BSS problem by solving a sequence of weighted LASSO problems, with the weights determined progressively.

For generalized linear models where the log-likelihood is used as the goodness-of-fit measure, the BSS problem can be formulated as an optimization problem.
The goal is to minimize the negative log-likelihood (equivalent to maximizing the likelihood) with respect to the parameter vector β while constraining the number of selected features to be k.
This results in the following optimization problem:

    βminimize -logℒ(β; y)
          subject to  ∑_i=1^p 𝕀_{β_i ≠ 0}≤ k.

Here, k ∈1, …, p is a user-defined parameter that restricts the size of the subset.
We can formulate the whole problem as a mixed-integer convex optimization problem containing the constraint with the ℓ_0 pseudo-norm, which counts the number of non-zero entries in β.
To do so, we introduce p binary variables z_i that indicate whether the covariate β_i is selected for the model or not.
Note that β_0 denotes the intercept, which is always included in the model.
Now, the BSS problem becomes the following problem:

    [                            βminimize                         - logℒ(β; y);                           subject to -M z_i ≤β_i ≤ M z_i,    i = 1, …, p,;                                                          ∑_i=1^p z_i ≤ k,;                                                      β∈^p+1, z ∈{0, 1}^p. ]

Here, M ≥ ||β̂||_∞ is a constant that ensures that a coefficient β_i is zero if the corresponding binary variable z_i is zero.
In other words, z_i indicates whether β_i is included in the model. 
These types of constraints are often referred to as big-M constraints.
It is well documented that problems containing big-M constraints depend on a good choice of M.
If M is too small, the convergence to the same optimum as the original problem is not guaranteed. 
If M is too large, loss of accuracy and numerical instability may occur.
One should also be aware that choosing an arbitrarily large M results in an unnecessarily large feasible region for the LP relaxation <cit.>.
To address this challenge, <cit.> proposed a data-driven method to determine lower and upper bounds for β̂_i in linear regression.
In the following, we extend this approach to convex GLMs.
Let UB be an upper bound on the MLE of Problem <ref>. 
Then one can find lower and upper bounds by solving the following convex optimization problems:

    [         u_i^+ :=   βmaximize  β_i;            s. t. -logℒ(β; y) ≤UB, ]


    [        u_i^- :=  βminimize  β_i;           s. t. -logℒ(β; y) ≤UB ]

where u_i^+ is an upper bound and u_i^- is a lower bound to β̂_i.
Now, let M_i := max{|u_i^+|, |u_i^-|} and one can choose the big-M as M = imax  M_i.

This procedure involves solving several convex optimization problems and estimating an upper bound (UB) beforehand.
An alternative simpler approach, which we have also employed in our simulation studies, is to standardize the design matrix X and choose an M that works well for most settings.
We have found that a value of M = 100 works well for many data sets when using a standardized design matrix.


 §.§ Information criteria for feature subset selection

While BSS obtains the best subset of features for a fixed number of maximal active coefficients k (where a coefficient is said to be active if it is non-zero), 
information criteria such as the Akaike Information Criterion <cit.> or the Bayesian Information Criterion <cit.> are often used (typically in a second stage, see <cit.>)
to select the best model out of the candidate models with different number of active coefficients.
For linear regression, it is computationally advantageous not to optimize the AIC directly but to use Mallows' C_p statistic (CP) <cit.> instead. 
<cit.> show that for linear regression, the AIC and CP are equivalent in the sense that both reach their minimum objective value with the same set of active coefficients.

Instead of having a two stage procedure, in FSS we will replace the objective function in Equation <ref> by an IC. Although many other ICs exist, we will focus on the AIC and BIC in this paper, as they are the most commonly used in practice.
For linear regression, however, we use the computationally advantageous CP to replace the AIC and a modified CP (CP2) to replace the BIC.
The proof of the equivalence of BIC and CP2 can be found in <ref>.


The Akaike Information Criterion (AIC) is defined as follows:

    AIC = 2k - 2 log(ℒ)

where ℒ is the likelihood function and k represents the number of active covariates.
Similarly, the Bayesian Information Criterion (BIC) is defined as:

    BIC = k log(n) - 2 log(ℒ).

The Mallows' C_p statistic (CP) for linear regression is defined by:

    CP = 2k + 1/σ^2 ||y - X β||_2^2 - n.

Similarly, the modified CP, which is equivalent to the BIC, is defined as follows:

    CP2 = k log(n) + 1/σ^2 ||y - X β||_2^2 - n.

Here σ^2 is the residual variance.
The following observation, which is sometimes referred to as monotonicity of the GOF measure, applies to AIC, BIC, CP and CP2.
When we encounter two models with the same likelihood value, then the model with fewer selected coefficients is better.
We use the unbiased estimator of the variance σ̂^2 = ||y - X β||^2_2/n - p from the full regression model.









§ ISSUES IN GLM OPTIMIZATION

In this section, we investigate common data settings of GLMs that cause optimizers to recover wrong solutions.
Two main reasons are causing optimizers to recover an incorrect solution for GLMs.
Firstly, for GLMs with a binomial response, separation in the data can cause the maximum likelihood estimate to contain infinite components.
Secondly, cases where the data exhibits strong multicollinearity can lead to instability in the estimation and in inconsistent signs in the regression coeffients.

Before introducing our approach, we present in the following some proposed approaches
in the literature to address these issues.





 §.§ Separation

<cit.> show that for logistic regression and probit models, the finiteness and uniqueness of the MLE are connected to the overlap
of the data.
They identify three different data settings, complete separation, quasi-complete separation and overlap and show that for logistic regression overlap is necessary and sufficient for the finiteness of the MLE.
<cit.> translate these criteria into a linear problem, which can be checked to verify the existence of the MLE.

Although, in theory, the solvers should be able to detect the unboundedness of the problem, we found in our experiments that all the solvers we tried failed to detect the unboundedness for GLMs with a binomial response.
Therefore, it is important to check in the simulation that, indeed, a solution exists;
otherwise, we would mainly compare the default tolerance settings of the solvers.

<cit.> suggest using the following linear program (LP) to verify that the solution exists:

    [                            βmaximize    ∑_i ∈ I x^⊤_i β - ∑_i ∈ J x^⊤_i β;                           subject to   x^⊤_i β≤ 0   ∀ i ∈ J = {i|y_i = 0};                                      x^⊤_i β≥ 0   ∀ i ∈ I =  {i|y_i = 1}.;                                      ]

If the solution is a zero vector, this verifies that the data is overlapping and that the solution of the corresponding logistic regression model is finite and unique.
In case the solution of the LP is unbounded, the data is separated and the MLE does not exist.

In the preparation of this paper, we found more than two examples of peer-reviewed articles where the authors did not check their data for separation, which led them to report results based on unbounded optimization problems.
This likely occurred as a consequence of the solver incorrectly signaling convergence.




 §.§ Multicollinearity

In the presence of strong collinearity, the matrix X^⊤ X in linear regression becomes ill-conditioned, leading to unstable estimates of the coefficients (β̂).
This instability can result in inflated estimates and inconsistent signs of the coefficients.
Different approaches to solving this problem have been proposed in the optimization and statistics literature.
One approach, as suggested by <cit.>, is to limit collinearity by incorporating the at most one constraint: 

    z_i + z_j ≤ 1   ∀ (i, j) ∈ℋ𝒞 = {(i, j) :  ν≤ |ρ_ij|},

where the binary z variables are the ones introduced in Equation <ref>, 
ν is a predefined constant and ρ_ij denotes Pearson's empirical correlation coefficient between the i-th and j-th columns of the design matrix X. 
Hence, ℋ𝒞 represents the set of highly correlated features.
This constraint limits the pairwise correlation by ensuring that, at most one of the variables among a pair with a correlation exceeding ν is selected in the regression model.

<cit.> relax this constraint to the sign coherence constraint (see Equation <ref>)

    -M(1-u_ij) ≤β_i,  sign(ρ_ij) β_j ≤ M u_ij  ∀ (i, j) ∈𝒢 = {(i, j) :  τ≤ |ρ_ij| }

to force coefficients of covariates with large pairwise multicollinearity to have coherent signs.
Hence, highly positively correlated features must have coefficients with the same sign, while highly negatively correlated features must have coefficients with opposite signs. 
Again, M is a sufficiently large enough constant and u_ij is a binary variable introduced to enforce the sign coherence.
One can see that for u_ij=1 and positive ρ_ij it holds that 0 ≤β_i, β_j ≤ M.
Similarly, for u_ij=0 and negative ρ_ij we have that -M ≤β_i ≤ 0 and 0 ≤β_j ≤ M.
Clearly, the sign coherence constraint is less restrictive than the at most one constraint.

On the other hand, a different but related approach in statistics is to assume that strongly correlated features should have similar estimates.
This assumption is utilized in models like Ridge regression <cit.> and the OSCAR model <cit.>.
In Ridge regression, a L_2 penalty term is added to the objective function, encouraging similar estimates for strongly correlated features. 
The OSCAR model goes further by enforcing exactly the same coefficient for strongly correlated features, inducing a clustering behavior among the coefficients.




§ SUGGESTED MODEL



 §.§ FSS for the Poisson model

We present in this section the formulation of the FSS problem for the Poisson model
using the AIC, which, to the best of our knowledge, has not been proposed before in the literature. The corresponding AIC and BIC formulations for linear, logistic, and Poisson regression can be found in <ref>.

We formulate the following mixed-integer conic program for feature subset selection:

    [                                     β, δminimize 2( ∑_j=1^p z_j ) - 2 ( ∑_i=1^n y_i x_i^⊤β - δ_i);                                       subject to         (x_i^⊤β, 1, δ_i) ∈_expp,    i = 1, …, n,;                                                              -M z_i ≤β_i ≤ M z_i,    i = 1, …, p,;                                                                        β∈^p+1, z ∈{0, 1}^p, δ∈^n. ]

This problem can be solved by off-the-shelf mixed-integer conic optimization solvers like ECOS <cit.> or MOSEK <cit.>.
It is worth noting that <cit.> have highlighted the concave but non-linear nature of the log-likelihood function for Poisson regression and proposed a piecewise-linear approximation method. On the other hand, we are the first to suggest utilizing conic optimization to solve the feature subset selection problem specifically for Poisson regression.
Similar problem formulations can be established for all family-link combinations introduced in <cit.>.




 §.§ Combined constraint for multicollinearity

To further enhance this FSS model to handle multicollinearity, we propose the so-called
combined constraint where we integrate the sign coherence constraint and a equal magnitude constraint as a unified criterion. Based on the idea of similar estimates for highly correlated features, we believe that an equal magnitude constraint would be beneficial in much the same way that the sign coherence constraint extends the at most one constraint.
Assuming that the design matrix X has been standardized, we can define the equal magnitude constraint as follows:

    β_i = sign(ρ_ij) β_j   ∀ (i, j) ∈ℋ𝒞 = {(i, j) :  ν≤ |ρ_ij| }

here ρ_ij again refers to the Pearson's correlation coefficient, and ν is a predefined constant threshold. 
This constraint ensures that the coefficients of strongly correlated features have the same magnitude but possibly different signs based on the correlation direction.


This combination allows us to automatically control pairwise multicollinearity 
in addition to the feature subset selection process.
The sign coherence constraint is applied to variables exhibiting a moderate to strong pairwise correlation, while the equal magnitude constraint is exclusively imposed on strongly correlated pairs.
By employing the equal magnitude constraint instead of the previous at most one constraint, users gain more insightful information about the underlying data structure, instead of making a near-random selection.



 §.§ Final model

By incorporating FSS with our novel combined constraint, we formulate the following optimization problem:

    [                                              βminimize                                                     IC;                                             subject to                   -M z_i ≤β_i ≤ M z_i,    i = 1, …, p,;                                                        -M(1-u_ij) ≤β_i,  sign(ρ_ij) β_j ≤ M u_ij  ∀ (i, j) ∈𝒢;                                                                           β_i = sign(ρ_ij) β_j   ∀ (i, j) ∈ℋ𝒞;                                                                             β∈^p+1, z ∈{0, 1}^p, u ∈{0,1}^|𝒢| ]

where 𝒢 = {(i, j) :  τ≤ |ρ_ij| < ν} and ℋ𝒞 = {(i, j) :  ν≤ |ρ_ij| }, with 0 ≤τ < ν≤ 1.
Again, M represents a sufficiently large positive constant.
This model simultaneously enforces coherent coefficient signs for moderately correlated features and equal coefficient magnitudes for highly correlated features.
The thresholds for medium and high correlations are defined by τ and ν, respectively. 
Moreover, before estimating any binomial model, we employ the linear program in <cit.> to identify the problems whose solution is actually not determinable.

It is worth noting that the objective function of our final model, as shown in Equation (<ref>), incorporates an arbitrary information criterion.
Consequently, our model can accommodate multiple information criteria, such as AIC, BIC, CP, and CP2, which can be expressed using the respective formulas in Equation (<ref>, <ref>, <ref>, <ref>).
This not only allows for greater flexibility, but also opens to the door to more customized models.
The model can be further extended by using further holistic constraints and expert knowledge can be seamlessly integrated into the model.
Our entire approach should be viewed as an additional tool in the modern data scientist's tool belt.




§ SIMULATION STUDIES

In our simulation study, we present two key findings.
Firstly, our approach demonstrates the ability to recover the true predictors.
Specifically, we achieve selection accuracy comparable to that of exact methods, highlighting the quality of our solutions.
Secondly, our newly integrated constraint proves successful in estimating variables within a multicollinearity context.
This finding emphasizes the effectiveness of our approach in overcoming challenges posed by multicollinearity and recovering accurate estimates.

All computational experiments were conducted on a Dell XPS15 laptop with an Intel Core i7–8750H CPU @ 2.20GHzx12 processor and 32 GB of RAM.
We utilized three mixed-integer optimization solvers: Gurobi 9.1.2 <cit.>, MOSEK 10.0.34 <cit.>, and ECOS 2.0.5 <cit.>.
The R package ROI was employed for representing the optimization problems.

To obtain the exact reference solutions, we utilized the R packages lmSubsets <cit.> and bestglm <cit.>.
The bestglm package can also solve linear regression problems, but we exclude it from the comparison as lmSubsets exhibits significantly faster performance.
There exist many more R packages for subset selection, such as glmulti <cit.>, L0Learn <cit.> or abess <cit.>, but only the selected two ensure that the solutions are indeed globally optimal.
In order to ensure a fair comparison, all solvers were restricted to utilizing only a single core.


 §.§ Simulation without multicollinearity

In this simulation, we compare the performance of the FSS formulations suggested in Equation (<ref>, <ref>, <ref>, <ref>, <ref>, <ref>) with special purpose solvers for FSS and BSS.

We use true positives (TP) and true negatives (TN) to calculate the selection accuracy. 
The true positives are the number of features j for which both the estimated coefficient (β̂_̂ĵ) and the true coefficient (β_j^true) are non-zero:

    TP(β) = |j : β̂_̂ĵ≠ 0, β_j^true≠ 0|.

Analogously, the true negatives represent the number of features j for which both the estimated coefficient (β̂_̂ĵ) and the true coefficient (β_j^true) are zero:

    TN(β) = |j : β̂_̂ĵ = 0, β_j^true = 0|.

Once TP and TN are calculated, the selection accuracy (A) is determined as the sum of TP and TN divided by the total number of potential features (p):

    A(β) = TP + TN/p.

The runtime provides an indication of the computational efficiency of the solvers.
At the same time, the selection accuracy measures how well the solvers are able to correctly identify the relevant and irrelevant features.
For linear regression, we employ the specialized solver lmSelect from the lmSubsets package, which utilizes a branch-and-bound strategy tailored for this problem.
For logistic and Poisson regression, we use the dedicated solver bestglm.
Unlike lmSubsets, bestglm employs complete enumeration and is only suitable for regression problems with a moderate number of features.
To avoid excessively long runtimes, bestglm restricts the maximum number of features (for non-Gaussian families) to 15.

Following the setting of <cit.> the simulation study adopts the following design: 
the design matrix X is generated from a multivariate normal distribution, with x_i ∼𝒩(0, Σ) for i = 1, …, n with a mean of zero and the covariance matrix Σ is the identity matrix I_p. 
Each scenario consists of 5 different runs, with n = 1000 observations.
The number of features p varies, and for each scenario, the first ⌈p/2⌉ coefficients of β are set to 1, while the remaining coefficients are set to 0.

In the linear regression setting, we generate 175 datasets.
The number of features p varies among {20, 25, 30, 35, 40, 45, 50}, and different standard deviations σ∈{0.05, 0.10, 0.50, 1.00, 5.00} are employed.
The response variable y_i is generated according to y_i = x_i^⊤β + ϵ_i, where ϵ_i follows a normal distribution with mean 0 and variance σ^2.
To avoid long runtimes in the brute force setting, we limit the allowed time for each run to 4200 seconds.

In the logistic and Poisson regression settings, we vary the number of features within p ∈{5, 10, 15, 20, 25, 30, 35, 40}.
The inverse link function is used to calculate μ_i = g^-1(η_i), where η_i = x_i^⊤β.
Depending on the model, the response y_i is sampled from either y_i ∼Binomial(1, p = μ_i) or y_i ∼Poisson(λ = μ_i), based on the underlying distribution assumption.

Tables <ref>–<ref> summarize the results of the computation time and selection accuracy for the different GLMs.

For linear regression, Table <ref> shows that all tested methods achieve a high selection accuracy.
However, in some cases, lmSelect finds models with a slightly lower AIC or BIC, including more features.
Upon inspecting the data, it is discovered that the true coefficients of the additionally selected features are all 0.
This discrepancy can be attributed to the fact that the residual variance σ^2 is only estimated, and as a result, the minima of AIC and Mallows' C_p may not coincide perfectly.
As for the computation time, Table <ref> indicates that our FSS model is only slightly slower than lmSelect. 
Compared to the enumerating BSS approach of trying all different values for the sparsity parameter k, the FSS approach is faster and scales better with the number of features.

Regarding logistic regression, FSS can recover the actual coefficients with high accuracy, as can be seen in Table <ref>. 
In terms of scalability, our FSS approach scales better than the complete enumeration of bestglm. See Table <ref> for the full timings.

Our findings for Poisson regression are similar to those for logistic regression. 
The computation time and selection quality results are summarized in Table <ref> and Table <ref>.
Here we achieve excellent selection accuracies and FSS scales better than complete enumeration.
For Poisson regression, we found that ECOS encounters numerical problems for the simulation setting with more than 20 features (see also Table <ref>).



 §.§ Simulation with multicollinearity

In this section, we compare the performance of FSS with and without pairwise correlation constraints. 
We use the mean squared error (MSE) as a performance measure of how well a method performs at recovering the estimates.
We consider three types of pairwise correlation constraints: the at most one constraint, the sign coherence constraint and the combined constraint.
Moreover, to establish a baseline, we also benchmark a model selected based on the information criterion alone without additional constraints which we call no constraint model.
Following the recommendation of <cit.> we set the threshold for the at most one constraint to ν = 0.7.
If the correlation between two features exceeds 0.7, at most, one of the features can be selected in the feature subset.
The sign coherence constraint evaluates two thresholds, τ = 0.5 and τ = 0.7.
This constraint ensures that if two features have a correlation magnitude above the threshold, they must have a coherent sign in the feature subset.
Thus, highly positively correlated features are forced to have the same sign, while highly negatively correlated features must have opposite signs.
The combined constraint combines the ideas of the equal magnitude constraint and the sign coherence constraint by enforcing the equal magnitude for highly correlated feature pairs and ensuring sign coherence for moderately correlated feature pairs.
To distinguish between moderately and highly correlated feature pairs, we used the thresholds τ = 0.5 and ν = 0.7.
 

It is important to note that the simulation study design will impact the study's outcome, especially when comparing pairwise correlation constraints.
By using the simulation setting proposed by <cit.>, we aim to provide a common ground for comparing the constraints. 
This simulation setting is widely used in statistics to simulate collinearity among features.

In this setting, the design matrix X is constructed as follows:

    x_ij = √((1 - α^2)) z_ij + α_i(p+1),    where  i = 1, …, n  and  j = 1, …, p.

Here, n = 100 represents the number of observations, p = 3 denotes the number of features, and z_ij is drawn from a normal distribution with mean 0 and standard deviation σ. The design matrix X is standardized, resulting in X^⊤ X being in correlation form.
The parameter α takes values from the set {0.7, 0.8, 0.9, 0.95, 0.99}.
In this particular setup, ρ_ij = α^2 gives the correlation between any two features.
Consequently, the pairwise correlations are ρ_ij = {0.49, 0.64, 0.81, 0.90, 0.98}.
The β coefficients are generated by computing the eigenvectors of X^⊤ X and selecting the eigenvector associated with the largest eigenvalue.

    η_i = β_0 + β_1 x_1i + β_2 x_2i + ⋯ + β_p x_pi = x_i^⊤β

Subsequently, Equation (<ref>) is employed to calculate η_i where β_0 is set to zero and μ_i is given by μ_i = g^-1(η_i). 
Finally, the response y_i is sampled either from
y_i ∼Normal(μ_i, σ) or
y_i ∼Binomial(1, μ_i) or
y_i ∼Poisson(μ_i).

This simulation setup allows for the generation of data with controlled correlations among features, enabling the evaluation and comparison of the different pairwise correlation constraints in the context of feature subset selection.
Table <ref>, Table <ref> and Table <ref> summarize the results for the different constraints and correlation settings for 1,000 simulations and 100 observations.

In the linear regression setting, the AIC and BIC exhibited almost identical results, leading us to report only the AIC outcomes. 
Analysis of Table <ref> suggests that when the standard deviation σ is very small (specifically, σ = 0.01) and the α values range from small to medium (0.7, 0.8, and 0.9), both the no constraint and the sign coherence constraint slightly outperform the combined constraint.
The at most one constraint performs best for small correlations.
This observation can be attributed to the decreased likelihood of the constraint becoming active at lower α values, rendering its performance similar to the one of the no constraint. 
Overall, the simulation results indicate that the combined constraint performs exceptionally well in this simulation setting, effectively harnessing the advantages of the sign coherence constraint and the equal magnitude constraint.

Table <ref> presents the results for the logistic regression setting, where the generated μ_i values implicitly determine the standard deviation.
In this scenario, the combined constraint outperformed all other constraints.
While the MSE values for all constraints were similar under low correlation conditions, the combined constraint exhibited significantly lower MSE values compared to the other constraints under strong correlation.

Similarly, Table <ref> showcases the results for the Poisson regression setting.
In this case, all constraints performed similarly under low correlation conditions.
However, as the correlation increased, the combined constraint emerged as advantageous in this simulation setting.
Furthermore, the results indicate that the MSE values for the no constraint and the sign coherence constraint were almost identical, suggesting that inconsistent signs were rare in this simulation setup.



§ CONCLUSION

This paper proposes an automated model selection approach by combining FSS with a correlation constraint for specific types of generalized linear models (GLMs), including linear, logistic, and Poisson regression.
Our method utilizes conic optimization and information criteria such as AIC or BIC for feature subset selection (FSS).
Moreover, it also enables the integration of additional constraints such as limiting pairwise correlation.
We have shown that our approach achieves high selection accuracy and improves computational efficiency as the dimensionality of the problem increases, outperforming naive enumeration methods.
The key contribution of our work is the development of a new mixed-integer conic programming (MICP) problem for model selection under multicollinearity.
Our experiments have demonstrated that our proposed combined constraint surpasses the existing collinearity constraints in the literature regarding performance.
Future research directions could explore incorporating additional regularization terms or leveraging linear approximations to speed up computation times. 





§ ACKNOWLEDGEMENTS

This work was supported by the Austrian Science Fund (FWF) under grant number ZK 35.






§ GLMS

The PDF of the Normal distribution can be expressed as:

    f(y; μ, σ^2) = 1/√(2 πσ^2)exp(-(y-μ)^2/2 σ^2) = exp((y μ - μ^2/2)/σ^2 + (-y^2/2 σ^2 - 1/2log(2 πσ^2)))

where one can see that it is indeed part of the EDM family with θ = μ, ϕ = σ^2, b(θ) = θ^2/2
and c(y, ϕ) = -y^2/2 ϕ - 1/2log(2 πϕ).

The PDF of the binomial distribution can be written as:

    f(y; p) = n  y p^y (1-p)^n-y = exp(logn  y + y log(p/1-p)+ n log(1-p) ).

Together with the link function θ = logit(p) = log(p/1-p) we get

    f(y; p) = exp(y θ + n log(1/1 + exp(θ)) + logn  y).

Here we can see that the binomial distribution is part of the EDM family with
θ = log(p/1-p), ϕ = 1, b(θ) = -n log(1/1 + exp(θ)) and c(y, ϕ) = logn  y.

The PDF of the Poisson distribution can be expressed by:

    f(y; μ) = μ^y exp(-μ)/y! = exp(y log(μ) - μ - log(y!)).

Now using the link function θ = log(μ) which is equal to exp(θ) = μ, we get

    f(y; μ) = exp(y θ - exp(θ) - log(y!)).

We can see that the Poisson distribution is also part of the EDM family
with θ = log(μ), ϕ = 1, b(θ) = exp(θ) and c(y, ϕ) = - log(y!).



§ MALLOWS C_P FOR BIC

The log-likelihood for linear regression is given by:

    logℒ(β; y) = - 1/2 σ^2 (y_i - x_i^⊤β)^2 - 1/2log(2 πσ^2).

Given the definition of BIC (in Equation (<ref>)) and CP2 (in Equation (<ref>)) we have

    [                                                   BIC =                                     k log(n) - 2 log(ℒ);                                                       = k log(n) + ∑_i=1^n 1/σ^2 (y_i - x_i^⊤β)^2 + log(2 πσ^2);                                                       =         k log(n) + n log(2 πσ^2) + ||y - X β ||_2^2/σ^2;                                                       =                                CP2 + n(log(2 πσ^2) + 1) ]

This implies that they reach their respective minima at the same coefficient values.



§ MIXED-INTEGER CONIC OPTIMIZATION PROBLEMS FOR FSS



 §.§ FSS for AIC



  §.§.§ Gaussian


    [                                                           β, ζminimize                                              2 (∑_j=1^p z_j ) + 1/σ^2ζ;                                                             subject to ( ζ + 1, ζ - 1, 2( y_1 -   x_1^⊤β ), …, 2( y_n -   x_n^⊤β )) ∈_soc^n+2;                                                                                                          -M z_i ≤β_i ≤ M z_i,    i = 1, …, p,;                                                                                                                      β∈^p+1, z ∈{0, 1}^p, ζ∈. ]




  §.§.§ Binomial


    [                                      β,δ,γminimize 2 ( ∑_j=1^p z_j ) - 2 ( ∑_i=1^n y_i x_i^⊤β - δ_i );                                         subject to             (δ_i, 1, 1+γ_i)∈𝒦_expp,    i = 1, …, n;                                                              ( x_i^⊤β, 1, γ_i)∈𝒦_expp,    i = 1, …, n;                                                                  -M z_i ≤β_i ≤ M z_i,    i = 1, …, p,;                                                                     β∈^p+1, z ∈{0, 1}^p, δ∈^n,  γ∈^n. ]




  §.§.§ Poisson

The optimization problem is already given in Equation (<ref>).



 §.§ FSS for BIC




  §.§.§ Gaussian


    [                                                           β, ζminimize                                         log(n) (∑_j=1^p z_j ) + 1/σ^2ζ;                                                             subject to ( ζ + 1, ζ - 1, 2( y_1 -   x_1^⊤β ), …, 2( y_n -   x_n^⊤β )) ∈_soc^n+2;                                                                                                          -M z_i ≤β_i ≤ M z_i,    i = 1, …, p,;                                                                                                                      β∈^p+1, z ∈{0, 1}^p, ζ∈. ]




  §.§.§ Binomial


    [                                           β,δ,γminimize log(n) ( ∑_j=1^p z_j ) - 2 ( ∑_i=1^n y_i x_i^⊤β - δ_i );                                              subject to                  (δ_i, 1, 1+γ_i)∈𝒦_expp,    i = 1, …, n;                                                                        ( x_i^⊤β, 1, γ_i)∈𝒦_expp,    i = 1, …, n;                                                                            -M z_i ≤β_i ≤ M z_i,    i = 1, …, p,;                                                                               β∈^p+1, z ∈{0, 1}^p, δ∈^n,  γ∈^n. ]




  §.§.§ Poisson


    [                                            β, δminimize log(n) ( ∑_j=1^p z_j ) - 2 ( ∑_i=1^n y_i x_i^⊤β - δ_i );                                              subject to                (x_i^⊤β, 1, δ_i) ∈_expp,    i = 1, …, n,;                                                                            -M z_i ≤β_i ≤ M z_i,    i = 1, …, p,;                                                                                      β∈^p+1, z ∈{0, 1}^p, δ∈^n. ]





§ RESULT TABLES


 
 

 
 

 
 

 
 
 






elsarticle-num-names 








