Integer constraints for enhancing interpretability in linear regression

SORT 44 (1) January-June 2020, 67-98 DOI: 10.2436/20.8080.02.95 Integer constraints for enhancing interpretability in linear regression Emilio Carrizosa1, Alba V. Olivares-Nadal2and Pepa Ram´ırez-Cobo3 Abstract One of the main challenges researchers face is to identify the most relevant features in a prediction model. As a consequence, many regularized methods seeking sparsity have flourished. Although sparse, their solutions may not be interpretable in the presence of spurious coefficients and correlated features. In this paper we aim to enhance interpretability in linear regression in presence of multicollinearity by: (i) forcing the sign of the estimated coefficients to be consistent with the sign of the correlations between predictors, and (ii) avoiding spurious coefficients so that only significant features are represented in the model. This will be addressed by modelling constraints and adding them to an optimization problem expressing some estimation procedure such as ordinary least squares or the lasso. The so-obtained constrained regression models will become Mixed Integer Quadratic Problems. The numerical experiments carried out on real and simulated datasets show that tightening the search space of some standard linear regression models by adding the constraints modelling (i) and/or (ii) help to improve the sparsity and interpretability of the solutions with competitive predictive quality. MSC: 62J05, 90C11. Keywords: Linear regression, Multicollinearity, Sparsity, Cardinality constraint, Mixed Integer Non Linear Programming. 1 Introduction A plethora of real world data involve multiple features interacting between them. As a consequence, one of the most common research challenges is trying to predict a variable by making use of attributes that are deterministic or easier to access. A widely studied tool to achieve this is the linear regression model Y=β0+β β βX+a(1) 1Institute of Mathematics of the University of Seville (IMUS). 2The University of Chicago Booth School of Business, 5751 S. Woodlawn Ave., Chicago, Illinois 60637. Email: [email protected] 3Department of Statistics and Operational Research, Universidad de C´adiz. Received: January 2019 Accepted: April 2020 68 Integer constraints for enhancing interpretability in linear regression where Y= (y1,...,yK)′contains the Krealizations of the random variable to be predicted, X∈RK×Ncontains the observations of the attributes X1,...,XNthat influence on Y, and a∈RKdenotes the error term. In practice, the coefficients β β βneed to be estimated and thus, the user needs to select an estimation method, which is usually derived from solving an optimization problem of the form: min β β βf(β β β) s.t β β β∈B(2) where Bdenotes the feasible region. Since the data collection technologies are improving altogether with communication systems and computers’ memories and processors, the dimension of the data sets to be handled is increasing drastically. As a consequence, nowadays researchers aim to return an interpretable output which explains the main interactions between the features that conform the pile of data and the dependent variables. Usually, this is understood as a problem of choosing the most relevant features for prediction (Friedman, Hastie and Tibshirani, 2001; Cai, Tsay and Chen, 2009; Hastie, Tibshirani and Wainwright, 2015). Sparse methods will yield solutions β β βin (2) with a large number of zero coefficients, in which only the most significant features are associated with the non-zeroes (Tibshirani, 1996; Hastie et al., 2015). Although sparsity may be a desirable property for our solution, we should take into account that other characteristics need to be sought in order to obtain a more interpretable output. First, correlated variables can provide highly variable estimated coefficients that make it difficult to understand the impact of a feature on the predictive variable. Second, spurious coefficients complicate the judgement of whether a feature is truly relevant for prediction or not. We will explain these two issues with further detail in what follows and motivate why we aim to alleviate them in this paper while still seeking for a sparse solution. It is known that ordinary least squares (OLS) provides solutions that may be highly dense. A good representative of a possibly sparse estimation method in the form (2) is the lasso (Tibshirani, 1996), which adds a ℓ1-norm penalization term to the OLS objective: f(β β β) = ky−β0−β β βXk2 2+λkβ β βk1(3) B=RN. The lasso encompasses OLS when the penalty parameter λequals zero, but when λincreases the solution becomes more sparse. The lasso is computationally feasible and, under certain technical conditions on the data matrix X, it enjoys good statistical properties, see Friedman et al. (2001); B¨uhlmann and van de Geer (2011). However the lasso also presents certain shortcomings well-documented in the literature (for a brief review, see Bertsimas et al. (2016) and the references therein). In particular, it is known that estimation through OLS or the lasso may be quite unstable in the presence of strong Emilio Carrizosa, Alba V. Olivares-Nadal and Pepa Ram´ırez-Cobo 69 collinearity on the data (Silvey, 1969; Sengupta and Bhimasankaram, 1997; Hesterberg et al., 2008). On one hand, the presence of correlated variables may yield a high variability in the estimated coefficients, complicating thus the interpretation of the results (Farrar and Glauber, 1967; Watson and Teelucksingh, 2002; Montgomery, Peck and Vining, 2012). On the other hand, a consequence of collinearity that leads to problems for interpreting the effect of the regressors is that two variables that are highly positively (negatively) correlated may have associated estimated coefficients with different (same) signs. This problem can be illustrated through the following numerical example in Hesterberg et al. (2008). Consider the diabetes database (Efron and Hastie, 2003), which consists of the measures of 10 variables (age, sex, body mass index, average blood pressure and six different blood serums) on 442 patients. The top panel of Figure 1 depicts the path of solutions of the lasso for this database; that is to say, the estimates of the coefficients β β βobtained are depicted against the different values of the penalty λ. As noted by Hesterberg et al. (2008), features tc and ldl (bottom left panel), have a correlation of 0.89. However, their estimated coefficients take opposite signs, which is in contradiction with their dependence degree. Similarly, the coefficients for variables hdl and tch (bottom right panel), which show a correlation of −0.73, have the same sign when estimated by the OLS (the case λ=0 in Figure 1). Hence it seems that the coefficients of highly correlated variables may take values that compensate each other. Finally, an additional inconsistency is that the sign of the estimated coefficient of hdl (squared blue line, left bottom panel) varies depending on the level of sparsity required. The negative effects of collinearity have been differently addressed in the literature. On one hand, some authors (Chatterjee and Hadi, 2015; Montgomery et al., 2012) suggest to remove variables that are highly correlated or unimportant, often carrying out significancy tests to determine if a variable can be discarded. However, the results of these tests may be misleading in the presence of strong collinearity (Watson and Teelucksingh, 2002). In this line, the recent paper by Bertsimas and King (2015) proposes to tighten the estimation procedure (2) by adding constraints that explicitly forbid the coefficients of variables with a high pairwise correlation to be simultaneously non-zero. Nonetheless, as it will be seen in Section 3.3.1, these approaches may be detrimental if highly correlated features own a strong predictive power. On the other hand, some authors encourage highly correlated predictors to be altogether in the model. The graph-guided fused lasso (GFlasso hereafter, proposed in Kim and Xing, 2009) encourages two highly correlated variables to have similar estimated coefficients by adding the penalization γ|βi−sign(ρi j)βj|to the lasso objective function, where ρi j denotes the correlation between Xiand Xj. The SRIG method (Sparse Regression Incorporating Graphical Structure Among Predictors, introduced in Yu and Liu, 2016), determines the value of the coefficient βjnot only by feature Xjbut also by all features Xisuch that ρi j is large in absolute value. Under certain technical conditions, the SRIG is endowed with nice properties that, for instance, ensure the recovery of the original model. However, as it will be shown in Section 3.2, these conditions may not be fulfilled in real databases, yielding outputs that may not improve the performance of the 70 Integer constraints for enhancing interpretability in linear regression −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 0.6 Lasso λ β 243.35 123.93 86.49 35.61 24.29 18.87 5.47 1.5 1.39 0 −0.4 −0.2 0.0 0.2 0.4 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 Lasso λ β 889.32 452.9 316.07 130.13 88.78 68.97 19.98 5.48 5.09 0 tc ldl −0.4 −0.2 0.0 0.2 0.4 Lasso λ β −0.4 −0.2 0.0 0.2 0.4 Lasso λ β 889.32 452.9 316.07 130.13 88.78 68.97 19.98 5.48 5.09 0 hdl tch Figure 1:Top: path of solutions of the lasso for the diabetes database (the size of the coefficients β β βare depicted against the values of the penalty λ). Bottom: problematic paths. current methodologies. Our approach is less restrictive than the previous one, since it neither encourages the removal of variables nor the presence of groups of correlated features. Instead, we propose a constraint (called sign coherence constraint) that aims to elude the signs’ inconsistency phenomenon related to collinearity shown in Figure 1, avoiding that the coefficients of highly correlated variables compensate each other. This constraint provides more flexible models, as for some cases two highly correlated variables may appear altogether in the output, and for some other cases one feature of a highly correlated pair may be removed. This will be illustrated in Figure 5, Section 3.2. Our approach is not the first one to restrict the sign of the estimated coefficients in linear regression. For instance, Meinshausen (2013) makes use of non-negative least squares to recover the real sparsity pattern in high-dimensional data under certain conditions. Also, the LARS algorithm (Efron, Hastie, Johnstone and Tibshirani, 2004) emulates the lasso Emilio Carrizosa, Alba V. Olivares-Nadal and Pepa Ram´ırez-Cobo 71 solution by requiring the sign of the coefficients to match the sign of their correlation with the residuals. Another well-known example is the non-negative garotte (Breiman, 1995), which performs subset selection while forcing the signs of the coefficients to match the signs of the OLS estimates. To conclude, there are further remedies to alleviate collinearity issues, consisting of harvesting more observations (Sengupta and Bhimasankaram, 1997; Montgomery et al., 2012) or applying methods that decorrelate the data (Cao, Guo and Bouman, 2010; Massy, 1965). Nevertheless, the later approaches imply the transformation of the variables and thus complicate the interpretation of the final models with respect to the original features. More recently, optimization approaches bounding the Variance Inflation Factor (VIF) and condition number of the correlation matrix have also been proposed (Tamura et al., 2019; Jou, Huang and Cho, 2014; Tamura et al., 2017). On top of the unreliable interpretation of coefficients in presence of high correlations, the lasso also suffers a drawback, mitigated in this paper: the presence of spurious coefficients. For λ>0 the ℓ1-norm performs a shrinkage of the coefficients in the lasso solution that eventually attains sparsity as a side effect. However, as will be shown in the numerical section, the solutions of the lasso may be still dense for large datasets due to these spurious coefficients. In this paper we avoid this negative effect of shrinkage by defining a novel constraint (called significance constraint) that forces the estimated coefficients to be either zero or larger than a fixed value (to be tuned). In summary, in this paper we model two novel constraints which will tighten the search space for β β βin Problem (2). As a result, the interpretability of the solutions is improved since (i) the signs of the coefficients are coherent with the sign of the correlations between highly or moderately correlated predictors, and (ii) the shrinkage is combatted while avoiding spurious coefficients, which may lead to the annihilation of some coefficients, thus increasing the sparsity. As will be shown in the numerical experiments, such better interpretability is obtained without damaging the predictive power of the model. When discerning the suitability of these constraints for a particular database, the user should realize that constraints modelling (i) become inactive if no highly correlated predictors are found, while constraints expressing (ii) do if all the variables have nonspurious estimated coefficients. Hence, the user might want to analyse the correlations before adding the sign coherence constraints. However, we do recommend adding the significance constraint if a regularized method is used for estimation. The resulting optimization problems will belong to the class of Mixed Integer Quadratic Programs (MIQP), which have recently proven very suitable in different statistics problem as linear regression (Tamura et al., 2019; Bertsimas and King, 2015), time series (Carrizosa, Olivares-Nadal and Ram´ırez-Cobo, 2016), classification (Carrizosa, Nogales-G´omez and Morales, 2016; Carrizosa, Nogales-G´omez and Morales, 2017), or dimensionality reduction (Carrizosa and Guerrero, 2014). Indeed, Bertsimas and King (2015); Bertsimas et al. (2016) use MIQP theory to solve (in tractable way) the best 72 Integer constraints for enhancing interpretability in linear regression subset selection problem (Miller, 2002): f(β β β) = ky−β0−β β βXk2 2(4) B=β β β∈RN:kβ β βk0≤VT, where the ℓ0-norm is the cardinality function kβ β βk0=#(j:βj6=0)and kβ β βk0≤VT,(5) denotes the cardinality constraint which leads to attain the desired level of sparsity given by the value VT. In this work, the two novel constraints (sign coherence and significance constraints) will be combined with the cardinality constraint (5), so that sparsity is also achieved in addition to a better interpretability. The paper is structured as follows. In the next section we model the new constraints to be added to Problem (2) in order to enhance interpretability through mathematical programming. The numerical experiments are carried out in Section 3, where the estimation methods under comparison and the design of experiments are also discussed. The last section is devoted to concluding remarks and extensions. 2 Mathematical model formulation In this paper, it is our aim to enhance the interpretability of the outputs by replacing any estimation procedure in the form (2) by: min β β βf(β β β) s.t β β β∈B∩S(6) where Swill gather the proposed constraints. Tightening an estimation procedure by adding constraints, i.e., solving (6) instead of (2), has already been considered in the literature in order to improve the performance of linear regression estimation methods like (2), see for example Bertsimas and King (2015). In this section we model the tightening set Sby defining constraints that can be added to a classic (possibly sparse) linear regression estimation method (2), in order to enhance the interpretability of the outcome as well as improving its sparsity. As commented in the previous section, the first novel constraint, called sign coherence constraint, imposes coherence between the signs of the estimated coefficients and the signs of large pairwise correlations between predictors. The second novel constraint, the so-called significance constraint, allows only for truly significant features to be considered in the model. The idea is that the user should feel free to add any of these constraints, when compatible, to her selected estimation method given by (2), yielding (6). Emilio Carrizosa, Alba V. Olivares-Nadal and Pepa Ram´ırez-Cobo 73 2.1 The sign coherence constraint The presence of correlated variables in the data is demonstrated to lead to undesired consequences, such as a high variability on the estimated coefficients and the sign inconsistencies explained in the introduction; see e.g. Bartholomew et al. (2008). As it was commented, the traditional procedure to avoid these undesired behaviour consists of removing highly correlated variables. In particular, Bertsimas and King (2015) forbids two highly correlated variables to be simultaneously non-zero. Specifically, the following pairwise correlation constraints are modelled γi+γj≤1∀(i,j)∈Ωη,(7) where Ωη={(i,j):|ρi j| ≥ η}is the set of pairs of features considered to be highly correlated, and γi,γjare defined as γj=(1 if βj6=0 0 if βj=0. (8) From now on, the constraint defined by Bertsimas and King (2015) and stated as (7) will be called correlation constraint. In contrast to Bertsimas and King (2015), we propose here a less restrictive approach that allows two highly correlated variables to be in the model at the same time, but forbids misleading interpretations and misrepresentative coefficients. Our aim is to avoid sign inconsistencies while allowing the model to include two correlated variables if they contribute to improve or maintain the prediction quality. Therefore, we propose to model constraints that avoid the compensation of coefficients for correlated variables. Under the light of the example illustrated in Figure 1, these are the requirements we aim to gather when modelling the sign coherence constraint: 1. The coefficients of two features that are moderately or highly positively correlated must have the same sign. 2. The coefficients of two features that are moderately or highly negatively correlated must have opposite signs. In order to model these constraints, we introduce the following binary variables: ν+ j=(1 if βj>0 0 otherwise ν− j=(1 if βj<0 0 otherwise 74 Integer constraints for enhancing interpretability in linear regression Therefore, the previous requirements 1-2 can be easily formulated as constraints as: ν+ i+ν− j≤1∀(i,j)∈Ω+ α(9) ν− i+ν+ j≤1∀(i,j)∈Ω+ α(10) ν+ i+ν+ j≤1∀(i,j)∈Ω− α(11) ν− i+ν− j≤1∀(i,j)∈Ω− α(12) where Ω+ αand Ω− αare the sets of pairs of features that are moderately or highly correlated, expressed as Ω+ α={(i,j):ρi j ≥α}and Ω− α={(i,j):ρi j ≤ −α}. That is to say, constraints (9)-(10) mean that, if two variables i,jare highly positively correlated (i.e. (i,j)∈Ω+ α), then we do not allow one of the coefficients to be positive and the other negative. Similarly, constraints (9)-(10) imply that, if two variables i,jare highly negatively correlated (i.e. (i,j)∈Ω− α), we forbid their coefficients to be both positive or both negative. Note that variables ν+ j,ν− jare linked with γj, defined in Equation (8), as follows: γj=ν+ j+ν− j, and thus the cardinality constraint (5) can be also written as: N ∑ j=1 (ν+ j+ν− j)≤VT. In order to illustrate the impact of these constraints we compare the path of solutions depicted in Figure 1 for the lasso applied to the diabetes dataset, against the lasso tightened with the sign coherence constraints, as depicted in Figure 2. As it can be observed, the use of constraints (9)-(12) to tighten the feasible region of the lasso might avoid the sign of a coefficient to vary depending on the level of sparsity required, easing the interpretation of the impact of the predictors over the response variable. 2.2 The significance constraint In this section we formulate a novel constraint that helps combatting the negative effects of shrinkage of the lasso while discarding spurious coefficients. The idea is to allow only for significant variables to be represented in the model and to improve the sparsity of the output. Intuitively, large coefficients are identified with the significance of a feature once the data are normalized. Following this reasoning we propose to establish a threshold of significance that a feature must be able to exceed to be allowed in the model. We model the significance constraint as follows: |βj| ∈ {0} ∪[ǫ,+∞)j=1,...,N(13) Emilio Carrizosa, Alba V. Olivares-Nadal and Pepa Ram´ırez-Cobo 75 −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β −0.4 −0.2 0.0 0.2 0.4 Lasso+sign coherence λ β 8.9e+02 4.5e+02 3.2e+02 1.3e+02 8.9e+01 6.9e+01 2.0e+01 5.5e+00 5.1e+00 0.0e+00 Figure 2:Path of solutions of the tightened lasso by the sign coherence constraint for the diabetes database. where ǫis called here the significance threshold, to be fixed by the user or to be tuned (see Section 3.1). In (13), a forbidden region in (0,ǫ)is defined with the aim to avoid shrinkage and to forbid spurious coefficients in the solution. This constraint was already used by Carrizosa et al. (2016) to discover potential causalities in multivariate time series. Note that using the binary variables ν+ jand ν− jdefined in Section 2.1, constraints (13) can be expressed in a more manageable form via two sets of linear constraints βj≥ǫν+ j−ν− jM∀j=1,...,N βj≤ −ǫν− j+ν+ jM∀j=1,...,N(14) where Mis a large constant. This big M, often appearing when modelling problems with integer variables, is large enough so it does not exclude reasonable values of the parameters βj(see, e.g. Camm, Raturi and Tsubakitani, 1990). In order to clarify the effect of the significance constraint (13), consider the heat map given by Figure 3. The left panel depicts the values of the estimated coefficients for the lasso, in the first column, and the lasso with the significance constraint (taking ǫ=0.3), in the second column, for the golf2009 database (Winner, 2016). The right panel represents the values of the estimated coefficients for the OLS and its counterpart tightened with the significance constraint (taking ǫ=0.05)for the compact database (Torgo, 2016). Such datasets will be described with further details in Section 3. The colour represents the sign of the coefficients βj(blue for negative, red for positive) and the intensity is related to the magnitude of such coefficients. In Figure 3 it can be observed that adding significance constraints establishes a clearer cut between zero and non-zero coefficients. Also note 82 Integer constraints for enhancing interpretability in linear regression set S1). In this heatmap, features X6and X10, with correlation 0.91, do not appear simultaneously in the same shuffle when adding the correlation constraints. Nonetheless, their coefficients are non-zero simultaneously in various shuffles when considering sign coherence constraints. β1 β2 β3 β4 β5 β6 β7 β8 β9 β10 β11 OLS ∩S1 OLS ∩S2 Lasso OLS Estimated coefficients for golf2009 database Figure 5:Heatmap representing the coefficients βiestimated by OLS, lasso and OLS tightened by adding correlation constraints (tightening set S2), or sign coherence constraints (tightening set S1) for golf2009 database. The above results were obtained for fixed values of ηand α, which were set to 0.6 and 0.8, respectively. These numerical examples followed the experimental design in Bertsimas and King (2015), who fixed the correlation threshold α. However, in Appendix B we explore the sensitivity of the tightening procedures to changes in ηand α. We conclude that, in general, the calibration of these parameters seem to yield less sparse solutions with a similar MSE. As a consequence the results shown in our numerical experiments disregard the calibration of the correlation thresholds. 3.2.2 The effect of the significance constraint Now we aim to study the impact of adding significance constraints (13) to the OLS and the lasso. This was briefly analysed in Section 2.2, where heat maps representing the estimated coefficients β β βwere represented in Figure 3. We observed that imposing a threshold ǫto the estimates may lead to more sparse solutions by avoiding spurious Emilio Carrizosa, Alba V. Olivares-Nadal and Pepa Ram´ırez-Cobo 83 coefficients and discarding unimportant variables. The first two rows of Table 3 display the results for the untightened OLS and its counterpart tightened with S3(significance constraint). Analogously, the last two rows display the results for the unrestricted and tightened lasso. In this table, we observe that tightening the feasible region of the OLS and the lasso by using the set S3always improves the sparsity of the output while usually attaining a competitive predictive quality. For instance, the significance constraints improves the MSE of the OLS and the lasso in a 4% and a 2.6% for yacht database, respectively, while reducing in 3.5 and 1.5 the number of non-zeroes. However, an exception is found in golf2009 dataset, where the novel constraint worsens the accuracy of the OLS and the lasso, although yielding 3 more zeroes in both cases. Table 3:Predictive quality (MSE) and sparsity degree (NZ) for the baseline methods and the approach tightened by the significance constraint. Cpu Yacht Whitewine Redwine Golf2008 Golf2009 Compact MSE NZ MSE NZ MSE NZ MSE NZ MSE NZ MSE NZ MSE NZ OLS 1.000 6 1.000 6 1.000 11 1.000 11 1.000 6 1.000 11 1.000 21 S30.978 4.5 0.960 2.5 0.999 8 1.000 7.5 1.007 5 1.205 8 1.002 15 lasso 1.000 5 0.959 2.5 1.000 10.5 0.997 9.5 1.000 4 0.798 9.5 1.000 20.5 S30.989 4 0.934 1 0.999 8 1.001 7 0.987 3 0.936 6.5 1.003 13 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement 2 4 6 8 10 −0.06 −0.04 −0.02 0.00 0.02 λ Improvement Median improvement over the MSE in whitewine database 0.05 0.06 0.08 0.1 0.125 0.15 0.175 0.2 0.25 0.3 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement 2 4 6 8 10 −0.10 −0.05 0.00 0.05 0.10 λ Improvement Median improvement over the MSE in compact database 0.05 0.06 0.08 0.1 0.125 0.15 0.175 0.2 0.25 0.3 Figure 6:Improvements over the MSE of the lasso when tightened via the significance constraint. Since it is not straightforward to choose a grid of thresholds ǫto calibrate from, we will try to gain some intuition by studying the improvement of the predictive quality when adding the significance constraints for each value of ǫconsidered. To this aim, Figure 6 shows the median improvement on the lasso MSE for each value of the penalty λand each threshold in the proposed grid for whitewine (left panel) and compact (right panel) database. 84 Integer constraints for enhancing interpretability in linear regression From Figure 6 we observe that the behaviour for the largest ǫ(0.3, 0.25) can be more effective specially for the larges values of λ. That is to say, a large ǫmay help combating the strong shrinkage of the lasso when highly sparse solutions are sought. Still large but more conservative values of the threshold (ǫ=0.2,0.175,0.15) seem to also improve the MSE of the lasso with large λ, also providing less extreme behaviours than the choices ǫ=0.3,0.25. Finally, the smallest values of the threshold (ǫ=0.06,0.05) may slightly improve the lasso with small λand the OLS. In conclusion, there is no straightforward a priori choice for the parameter ǫ, which should be calibrated. Nevertheless, if the user is seeking a highly sparse solution (i.e., the user is choosing a high penalty λ) it seems advisable to choose larger values for ǫin order to combat the shrinkage more effectively. On the other hand, when estimating via OLS (or lasso with small values of λ) the user might want to focus on a grid with a majority of small values of ǫ. 3.2.3 Global performance Finally, we show in Table 4 the results when our novel constraints are jointly considered in combination to the cardinality constraint, or equivalently, when the set S4(cardinality + sign coherence + significance constraints) is used to tighten the OLS or lasso approaches. For comparison reasons, the table also shows the results under the tightening set S5(cardinality + correlation constraints) proposed in Bertsimas and King (2015). Table 4 also displays the predictive quality and number of non-zero coefficients for the Enet, the SRIG method in Yu and Liu (2016) and the GFlasso method in Kim and Xing (2009). Table 4:Predictive quality (MSE) and sparsity degree (NZ) for the baseline estimation methods (OLS, lasso, SRIG, Enet and GFlasso), and OLS and lasso tightened with S4(cardinality + novel constraints) or S5(cardinality + correlation constraint), for the real datasets. Method Cpu Yacht Whitewine Redwine Golf2008 Golf2009 Compact MSE NZ MSE NZ MSE NZ MSE NZ MSE NZ MSE NZ MSE NZ OLS 1.000 6 1.000 6 1.000 11 1.000 11 1.000 6 1.000 11 1.000 21 S40.990 4 0.934 1 1.011 6 1.002 6 0.999 3 0.963 4 1.012 9 S51.002 4 0.934 1 1.008 8 1.000 6 1.006 3 1.024 5 1.016 11 lasso 1.000 5 0.959 2.5 1.000 10.5 0.997 9.5 1.000 4 0.798 9.5 1.000 20.5 S40.993 4 0.934 1 1.013 6 0.998 6 0.982 3 0.971 4.5 1.013 9 S51.049 4 0.934 1 1.008 7.5 0.995 6 0.988 3 1.045 4 1.015 11.5 SRIG 0.988 6 0.942 2 1.000 11 0.999 11 0.983 6 1.016 11 0.999 21 Enet 0.917 5.5 0.948 2 1.000 11 0.998 11 0.994 4 0.805 9.5 1.000 21 GFlasso 0.972 6 0.959 2.5 1.000 11 0.996 11 1.000 4 0.866 9.5 1.001 20 First, we analyse the performance of the baseline estimation procedures (OLS and lasso) against their tightened counterparts. We can conclude that reducing the search space of the coefficients β β βby intersecting with either S4or S5always improves the sparsity of the solutions. Moreover, the predictive quality is usually similar to that of the OLS and the lasso. In particular, for yacht database, both tightened approaches Emilio Carrizosa, Alba V. Olivares-Nadal and Pepa Ram´ırez-Cobo 85 improve the MSE of the OLS and the lasso estimates by a 6.6% and a 2.6%, respectively. However, these tightened procedures worsens the accuracy of the lasso in the golf2009 database. Second, we compare the performance of the two tightening sets. From the table it can be observed that they deliver a similar accuracy-sparsity trade-off for 3 out of the 7 real datasets (yacht, redwine, golf2008). For the remaining databases, the approaches attain different trade-offs between sparsity and predictive quality. Indeed, the novel set S4maintains the sparsity attained by the set S5og Bertsimas and King (2015) on the cpu database, while slightly improving the MSE in a 1.2% and a 5.3% for the tightened OLS and lasso, respectively. In contrast, the proposed tightening set provides more sparse solutions, with a similar predictive quality for whitewine and compact datasets. Finally, S4improves the MSE of the OLS and the lasso in a 6% and 7.1% in golf2009 database, but it provides a slightly more dense solution for the lasso. Note that the best accuracy for the lasso in this database was attained when adding exclusively sign coherence constraints (see Table 2), although the solution provided was more dense. Third, both the SRIG and GFlasso are clearly outperformed by the tightened approaches. On top of this, for the real datasets considered here there is no guarantee that the necessary assumptions to preserve the theoretical properties of the SRIG are fulfilled. In order to compare the performance of our approach against these methods under a more favourable scenario for the later, in the next section we replicate the simulation study of Yu and Liu (2016), hence assuring the non-violation of the conditions for the SRIG method. 3.3 Simulations In the previous section, the behaviour of the novel constraints for the case of real datasets with a small number of predictors was analysed. In this section we aim to examine the sensitivity of the tightened procedures under various settings. First, we will simulate data as in Yu and Liu (2016) to understand the behaviour of the methodology under different correlation structures and for different sizes of the training sample. Second, we aim to test the proposed methodology for larger datasets and for different correlations intensities. To do so we simulate data as in Bertsimas and King (2015). 3.3.1 Sensitivity to correlation structure and training sample size In this section we aim to test the proposed methodology for datasets simulated following the three examples and training sizes described in Yu and Liu (2016). Ten instances have been generated for each example and for training and testing sizes of 40, 80 and 120. For both Examples 1 and 2, βi=3 for i=1,..,15, and βi=0 for i=16,...,100. In Example 1, however, the predictors were generated as follows: Xi=Zj+0.4ai X,Zj∼N(0,1),ai X∼N(0,1) 86 Integer constraints for enhancing interpretability in linear regression for 5(j−1) + 1≤i≤5jand j=1,2,3. For i>15, Xi∼N(0,1). In Example 2, the k-th vector of observations X1 k,...,XN kwas generated following a multivariate normal distribution with zero mean and covariance matrix Σ= (σi j), where σi j =0.5|i−j|. Analogously, in Example 3 the vector of observations was also drawn from a multivariate normal distribution with zero mean but with covariance matrix Σ= (B+ψI)−1, where bi j =0 for i=j, and bi j =0.5δi j,δi j ∼Be(0.05)otherwise. Parameter ψis fixed so that the condition number of Σ−1equals N. The real coefficients are β β β=Σ−1Σxy, where Σxy is the cross-covariance vector whose elements equal 10 for the four predictors with the largest degrees and 0 otherwise. These three examples have different structures of correlation and, as it will be seen later, this may influence on the performance of the approaches making use of constraints taking into account pairwise correlations. As an illustration, Figure 8 of the Appendix C of the Supplementary Material displays the pairs of variables appearing in the correlation constraints of Bertsimas and King (2015) and our sign coherence constraints. As done in Bertsimas and King (2015), the grid of values of the parameter λto be tuned for the tightened MIQP with lasso objective function is logarithmically generated in the interval (0,λmax], where λmax is the penalty provided by the lars for which only one coefficient is non-zero. Analogously to the real datasets results, Table 5 shows the median MSE and number of non-zeroes (NZ). In particular, for each example (rows) and training sizes (columns), each row shows the results obtained for an estimation method, namely the OLS, the lasso, and their tightened counterparts, which take the form of Problem (6) with tightening sets S4(cardinality + sign coherence + significance constraints) and S5(cardinality + correlation constraints). The last rows for each example display the results for the methods dealing with correlated variables: the Enet, the SRIG, and the GFlasso. To make it easier to discuss these results, in Figure 7 we have represented the MSE against the number of nonzeroes for the three simulated examples of Table 5. The different estimation methods have been assigned different colours, the solid items representing the approaches making use of our proposed tightening set. The diversity of training samples have been represented by unalike symbols. In the top panel of Figure 7 we can observe that, in Example 1, the MIQPs (i.e., the approaches proposed in Bertsimas and King (2015) and in this paper) tend to attain more sparse solutions than the continuous optimization methods (i.e., the lasso, SRIG and GFlasso). As it can be observed in Figure 8 in the Appendix, the true generating model contains predictors that hold a high pairwise correlation and, therefore, are forbidden to appear simultaneously in the outputs yielded by the tightening set of Bertsimas and King (2015). As a consequence, the approaches making use of this set S5are not able to recover the original graph, hence delivering more sparse solutions but with worse predictive quality than our tightening set. In constrast, the sign coherence constraints flexibility allows the simultaneous presence of all the variables in the generating model. As a consequence, S4yields outputs with similar performance to that of the SRIG, that is the best approach amongst those based on continuous optimization models. Emilio Carrizosa, Alba V. Olivares-Nadal and Pepa Ram´ırez-Cobo 87 ll 0 10 20 30 40 0.0 0.2 0.4 0.6 0.8 1.0 1.2 NZ MSE 0 10 20 30 40 0.0 0.2 0.4 0.6 0.8 1.0 1.2 NZ MSE 0 10 20 30 40 0.0 0.2 0.4 0.6 0.8 1.0 1.2 NZ MSE OLS BS OLS S4 OLS S5 Lasso Lasso BS Lasso S4 Lasso S5 SRIG Enet GFlasso l Ntrain=40 Ntrain=80 Ntrain=120 Example 1 ll 0 20 40 60 80 0.0 0.2 0.4 0.6 0.8 1.0 1.2 NZ MSE 0 20 40 60 80 0.0 0.2 0.4 0.6 0.8 1.0 1.2 NZ MSE 0 20 40 60 80 0.0 0.2 0.4 0.6 0.8 1.0 1.2 NZ MSE OLS BS OLS S4 OLS S5 Lasso Lasso BS Lasso S4 Lasso S5 SRIG Enet GFlasso l Ntrain=40 Ntrain=80 Ntrain=120 Example 2 l l 0 10 20 30 40 50 0.0 0.1 0.2 0.3 0.4 0.5 0.6 NZ MSE 0 10 20 30 40 50 0.0 0.1 0.2 0.3 0.4 0.5 0.6 NZ MSE 0 10 20 30 40 50 0.0 0.1 0.2 0.3 0.4 0.5 0.6 NZ MSE OLS BS OLS S4 OLS S5 Lasso Lasso BS Lasso S4 Lasso S5 SRIG Enet GFlasso l Ntrain=40 Ntrain=80 Ntrain=120 Example 3 Figure 7:Median MSE and NZ for each method in the simulated datasets of Yu and Liu (2016). In the bottom panels of Figure 7 we observe that the MIQP approaches are clustered together. For both Examples 2 and 3 it is evident that the methods based on integer optimization improve the sparsity of the outputs with no damage to the predictive quality. More specifically, in Example 3 all the methods based on continuous optimization conform a unique cluster with higher density and similar or slightly worse MSE than the MIQP approaches. Nonetheless, in Example 2 the SRIG and the lasso attain more sparse solutions than the rest of continuous optimization approaches, yet still yielding more non-zeroes than the tightened procedures. In this case, all the approaches attain a similar accuracy. 88 Integer constraints for enhancing interpretability in linear regression Table 5:Predictive quality (MSE) and sparsity degree (NZ) for the baseline estimation methods (OLS, lasso, SRIG, Enet and GFlasso), and OLS and lasso tightened with S4(cardinality + novel constraints) or S5(cardinality + correlation constraint), for the simulated datasets with N =100. Ntrain=40 Ntrain=80 Ntrain=120 MSE NZ MSE NZ MSE NZ OLS 1.000 100.0 1.000 100.0 1.000 100.0 S40.689 7.5 0.240 7.5 0.081 12.0 S51.138 3.5 0.767 3.0 0.579 3.0 lasso 0.358 24.0 0.131 27.5 0.043 31.0 S40.306 14.0 0.154 17.0 0.062 20.0 S51.231 4.5 0.731 4.0 0.564 3.5 SRIG 0.082 15.0 0.050 15.0 0.082 15.0 Enet 0.351 27.0 0.111 32.0 0.351 27.0 Example 1 GFlasso 0.188 26.0 0.064 31.5 0.188 26.0 OLS 1.000 100.0 1.000 100.0 1.000 100.0 S41.380 5.0 0.197 6.5 0.133 9.0 S51.560 4.0 0.206 6.0 0.142 9.0 lasso 0.756 27.5 0.096 30.0 0.080 31.5 S41.160 7.0 0.115 12.0 0.105 15.5 S50.841 18.5 0.115 15.5 0.102 15.5 SRIG 0.913 9.5 0.149 28.5 0.111 28.5 Enet 0.862 53.0 0.157 53.0 0.125 50.0 Example 2 GFlasso 0.638 38.5 0.053 80.0 0.038 53.5 OLS 1.000 100.0 1.000 100.0 1.000 100.0 S40.412 5.0 0.101 13.5 0.106 15.0 S50.424 2.0 0.110 8.5 0.101 11.5 lasso 0.494 18.5 0.135 44.0 0.151 47.5 S40.391 7.5 0.083 18.0 0.097 14.5 S50.325 7.5 0.078 13.5 0.097 16.5 SRIG 0.537 34.5 0.110 43.5 0.185 51.0 Enet 0.525 28.0 0.148 45.5 0.165 49.5 Example 3 GFlasso 0.494 23.5 0.135 44.0 0.151 51.0 3.3.2 Scalability and sensitivity to correlation intensity In this section we aim to test the proposed methodology for larger datasets with diverse correlation intensities. As the overall design of experiments, the synthetic generation of the data is similar to that of Bertsimas and King (2015). The k-th vector of observations X1 k,...,XN kwas generated following a multivariate normal distribution with zero mean and covariance matrix Σ= (σi j), where σi j =ρ|i−j|. In particular, we chose ρ=−0.9 and ρ=−0.5 so as to test the performance of the constraints under highly and moderate correlations. The regression model is taken in small dimension, but with quite a number of irrelevant covariates. More precisely, the number of features was set to 500, only 10 of Emilio Carrizosa, Alba V. Olivares-Nadal and Pepa Ram´ırez-Cobo 89 which corresponding to explanatory variables, the remaining 490 being noise. For each value of ρ, 10 instances were generated as follows. The βiwere uniformly generated in the interval (−2,2)for isuch that imodulus N/10 =0. The response was generated following (1), with β0=0 and the error terms i.i.d. following a normal distribution with zero mean and variance as in Bertsimas and King (2015). As done in Section 3.3.1, the sequence of λhas been logarithmically generated. Table 6 shows the median MSE and number of non-zeroes (NZ) of the solutions for the OLS, the lasso and their tightened counterparts, as well as the SRIG, Enet and GFlasso. Table 6:Predictive quality (MSE) and sparsity degree (NZ) for the baseline estimation methods (OLS, lasso, SRIG, Enet and GFlasso), and OLS and lasso tightened with S4(cardinality + novel constraints) or S5(cardinality + correlation constraint), for the simulated datasets with N =500. r=0 5–.r=–0 9. OLS 1.000 500 1.000 500 S40.745 6 0.732 6 S50.849 6 1.474 4 lasso 0.535 80.5 0.533 183 S40.526 10.5 0.528 11 S50.501 13 0.807 19 SRIG 0.630 34 0.792 102 Enet 0.535 50 0.534 51.5 GFlasso 0.537 70.5 0.533 54.5 MSE NZMSE NZ Note that, for the simulated data with moderately correlated features (ρ=−0.5), the correlation constraints (7) and the sign coherence constraints (9)-(12) are inactive, since the highest pairwise correlation in absolute value is roughly 0.5. In this case, both sets of tightening constraints help to considerably improve the sparsity of the baseline estimation procedures. Indeed, the density of the OLS is drastically reduced by increasing in 494 the zeroes of the output, while its predictive quality is also substantially improved. However, the novel tightening set S4provides an accuracy 12.2% better than the attained with the set S5proposed in Bertsimas and King (2015). On the other hand, the benchmark sparse regression method, the lasso, attains a median of 80.5 nonzeroes, while its tightening counterparts produce much more sparse solutions with better MSEs. Although the predictive quality of the outputs of the tightened procedure with S5 is slightly better than that obtained with S4, this comes at the price of yielding more dense solutions. More generally, the methods based on MIQP solvers clearly outperform the approaches relying on continuous optimization techniques: the former manage to considerably reduce the sparsity of the later with a slightly better predictive quality. Amongst the later approaches, the GFlasso is outperformed by the Enet which, with a similar accuracy, provides a much more sparse solution. Nonetheless, the most sparse of these methods is the SRIG, although yielding the worst predictive quality. 90 Integer constraints for enhancing interpretability in linear regression Table 7:Predictive quality (MSE) and sparsity degree (NZ) for the baseline estimation methods (OLS, lasso, SRIG, Enet and GFlasso), and OLS and lasso tightened with S4(cardinality + novel constraints) or S5(cardinality + correlation constraint), for the simulated datasets with ρ=−0.9and N =50,500,1000. N=50 N500=N1000= MSE NZ MSE NZ MSE NZ OLS 1.000 50 1.000 500 1.000 1000 S40.662 8 0.732 6 1.035 5 S50.634 8 1.474 4 1.310 5.5 lasso 0.671 24.5 0.533 183 0.555 56 S40.637 8.5 0.528 11 0.904 8 S50.595 19 0.807 19 0.406 10 SRIG 0.127 42 0.792 102 0.835 103.5 Enet 0.108 26 0.534 51.5 0.538 49.5 GFlasso 0.693 41 0.533 54.5 0.615 359 The simulated instances with highly correlated features, ρ=−0.9, show a similar behaviour. Nevertheless, the lasso provides significantly more dense outputs in this case, with a median of 183 non-zeroes. The tightening set S5proposed in Bertsimas and King (2015) worsens its predictive quality although improving substantially its sparsity. In contrast, the novel tightening set S4attains around eight more zeroes than the later while maintaining the lasso’s accuracy. Regarding the OLS, both tightened estimation methods reduce drastically the density of the solutions, although the S4obtains around two more non-zero coefficients than the tightening set S5. However, the price paid for two extra non-zeroes is shown in the accuracy of the solution: while adding cardinality, sign coherence and significance constraints improves the MSE of the OLS in a 26.8%, adding cardinality and correlation constraints instead worsens the accuracy in 47.4%. Analysing the results for the continuous optimization based methods, we observe that the best performance is yielded by the Enet and GFlasso which, with a similar MSE, significantly enhance the sparsity of the classic lasso. Nonetheless, when tightened with S4, the later attains outputs with comparable predictive quality and many more zero coefficients (around 40 more). In order to analyse the scalability of our methodology we have also simulated data as in Bertsimas and King (2015) with 50 and 1000 variables, and a maximum correlation of 0.9. For the later, the MIQP to be solved would have 3001 variables and more than 2000 constraints. As the size of the problem is considerably larger, we have allowed for a time limit of 40 seconds in this case, which was also the time limit considered for high dimensional data where N>K. The results are collected in Table 7. As it can be observed, MIQP approaches attain solutions that are considerably more sparse than their continuous counterparts while still delivering a good predictive quality. In particular, SRIG delivers the most dense outputs from the sparse continuous methods, while the GFlasso is outperformed both in terms of accuracy and sparsity. Regarding the MIQP, Emilio Carrizosa, Alba V. Olivares-Nadal and Pepa Ram´ırez-Cobo 91 they attain different accuracy-sparsity trade-offs. Although the proposed constraints attain better accuracy when combined with the OLS objective, this is not true when the regularization penalty λis positive. Nonetheless, in this case the tightening set S4yields more sparse outputs than S5. Summarizing, tightening the search space of the OLS and the lasso provides considerably more sparse solutions, although the sign coherence and significance constraints yield a better accuracy-sparsity trade-off than the correlation constraint. Indeed, the later can substantially worsen the predictive quality of the baseline methods in order to reduce the density of the outputs, while the former entails a more competitive MSE. 4 Concluding remarks The aim of this paper is to enhance the interpretability in a regression model without worsening its predictive quality. We assume we have a baseline regression estimation procedure based on solving an optimization problem (e.g. OLS or lasso), and then the underlying optimization problems are modified by adding new constraints to those defining the search space. These constraints avoid misleading estimators that may be obtained in the presence of highly correlated variables and detect the most important features for the prediction. In order to assess the impact of adding the two novel constraints over various estimation procedures, in our numerical experiments we consider the OLS and the lasso. The search space of β β βis reduced in these methods by using the tightening set composed by the sign coherence constraints (9)-(12) and/or the significance constraints (14), possibly in combination with the cardinality constraint (5). The first constraint forces the sign of the coefficients to be coherent with the sign of large and moderately large pairwise correlations between features, while the second avoids spurious coefficients and combats the shrinkage of regularized regression. We compare the performance of our tightening set, including all the proposed constraints, with the recent tightening set by Bertsimas and King (2015), which also defines a MIQP and includes the cardinality constraint (5) and the correlation constraint (7) which explicitly forbids two highly correlated variables to be simultaneously in the regression model. These methods are compared against other approaches also dealing with correlated variables but based on continuous optimization techniques: the Enet (Zou and Hastie, 2005), the SRIG (Yu and Liu, 2016), and the GFlasso (Kim and Xing, 2009). The results show that the novel constraints yield tractable optimization problems, solvable in short time by standard solvers, and may enhance the interpretability while often improving or maintaining the predictive quality and level of sparsity. More specifically, the MIQP approaches attain a different trade off between sparsity and predictive quality than the methods based on continuous optimization, usually yielding more sparse solution with similar or better MSE. Amongst the former methods, the novel constraints tend to improve the predictive quality of the outputs obtained with the tightening set proposed in Bertsimas and King (2015). 98 Integer constraints for enhancing interpretability in linear regression DAMPL code #DATA PARAMETERS param N; #Number of predictors param K; #Number of observations #SETS OF INDICES set Nvar:=1..N; set Nobs:=1..K; #THE MATRIX OF DATA param X {i in Nobs,j in 1..(N+1) };#Kx(N+1) matrix of data #Includes response variable in position N+1 #PARAMETERS OF THE METHODS param lambda; #Lasso penalty param NZ; #Upper-bound on the total number of non-zeroes param eps; #Significance threshold param M default 100; #Upper bound for the coefficients beta #SETS OF HIGHLY/MODERATELY CORRELATED VARIABLES set conjcorrpos dimen 2; #Positively correlated features set conjcorrneg dimen 2; #Negatively correlated features #VARIABLES var c ; #Intercept var beta {j in Nvar };#Slopes var nupos {j in Nvar }, binary; var nuneg {j in Nvar }, binary; var v {j in Nvar }>=0 ; #Auxiliar variables to express the absolute value #OBJECTIVE FUNCTION minimize fun: (1/N)*sum{i in Nobs }(X [i,p+1] -c-sum{j in Nvar }(beta[j]*X[i,j]))ˆ2 +(1/N)*lambda*sum{j in Nvar }v[j]; #CONSTRAINTS #SPARSITY CONSTRAINT subject to sparsity: sum{j in Nvar }(nuneg[j]+nupos[j])<=NZ; #SIGNIFICANCE CONSTRAINTS subject to significancepos {j in Nvar }: beta[j]>=eps*nupos[j]-nuneg[j]*M; subject to significanceneg {j in Nvar }: beta[j]<=-eps*nuneg[j]+nupos[j]*M; #SIGN COHERENCE CONSTRAINTS subject to coherencepos1 {(j,r) in conjcorrneg }: nupos[j]+nupos[r]<=1; subject to coherencepos2 {(j,r) in conjcorrneg }: nuneg[j]+nuneg[r]<=1; subject to coherenceneg1 {(j,r) in conjcorrpos }: nupos[j]+nuneg[r]<=1; subject to coherenceneg2 {(j,r) in conjcorrpos }: nuneg[j]+nupos[r]<=1; #AUXILIAR CONSTRAINTS subject to abs1 {j in Nvar }: v[j]>=beta[j]; subject to abs2 {j in Nvar }: v[j]>=-beta[j]; subject to sumnusj {j in Nvar }: nuneg[j]+nupos[j]<=1;