Interval Coded Scoring: a toolbox for interpretable scoring systems

Linear Programming ICS (lpICS)

This ICS core approach is based on the soft-margin SVM formulation (Vapnik, 1995). It can be restated and solved efficiently as a linear programming problem. Normally, MATLAB’s solver is used, but IBM’s CPLEX is selected if available, since it drastically improves execution time. The lpICS formulation is presented in matrix notation in Eq. (1). The same symbols as in Algorithm 1 have been used, with a tilde to refer to its specific use in lpICS. (1)${\displaystyle \begin{array}{cccccc}\hfill & \underset{\tilde{\mathbf{w}},b,\epsilon }{min}\hfill & \hfill & \parallel \mathbf{D}\tilde{\mathbf{w}}{\parallel }_1+\gamma {\epsilon }^T\mathbf{1},\hfill & \hfill & \mathbf{D}\in {\mathbb{R}}^{N_{ds}\times N_s},\tilde{\mathbf{w}}\in {\mathbb{R}}^{N_s}\hfill \\ \hfill & s.t.\hfill & \hfill & \mathbf{Y}\left(\tilde{\mathbf{Z}}\tilde{\mathbf{w}}+b\right)\ge \mathbf{1}-\epsilon ,\hfill & \hfill & \mathbf{Y}\in {\mathbb{R}}^{N\times N}{,}^{\text{lp}}\tilde{\mathbf{Z}}\in {\mathbb{B}}^{N\times N_s}\hfill \\ \hfill & \hfill & \hfill & \epsilon \ge \mathbf{0},\hfill & \hfill & \epsilon \in {\mathbb{R}}^N\hfill \end{array}.}$ Similarly to the original SVM formulation, it tries to fit a model $\tilde{\mathbf{Z}}\tilde{\mathbf{w}}+b$ on the class labels y in the diagonal matrix Y. The data matrix $\tilde{\mathbf{Z}}$ contains the binary transformed data, yielding N observations of N_s intervals. Furthermore, the constraint enforces correct classification, allowing slack ε with respect to the classification boundary. In the usual SVM objective function, the total slack is balanced against the classification margin through a hyperparameter γ, but lpICS minimizes $\parallel \mathbf{D}\tilde{\mathbf{w}}{\parallel }_1$ instead. This penalty term corresponds to Total Variation Minimization (Rudin, Osher & Fatemi, 1992). The binary matrix D encodes N_ds differences between coefficients of adjacent intervals of the same variable x^p or variable interaction (x^p₁, x^p₂). In the latter case, both the horizontal and vertical differences in the unvectorized matrix shape are taken into account. Hence, the sparsity-inducing ℓ₁ norm is applied on the difference of the weights. As a result, lpICS favours piecewise constant models as depicted in Fig. 1 because they have a sparse difference vector, or, equivalently, identical coefficients for adjacent intervals.

Data transformation. As mentioned in Algorithm 1, the matrix $\tilde{\mathbf{Z}}$ is the result of a transformation of the original data X guided by thresholds τ. The lpICS transformation $\tilde{\phi }$ supports both main effects (only a single variable x^p is affected) and interactions (x^p₁, x^p₂). The transformation of a single observation x_i is given as: (2)${\displaystyle \begin{array}{cc}\hfill & {\tilde{\mathbf{z}}}_i=\left[{\tilde{\mathbf{z}}}_i^1,{\tilde{\mathbf{z}}}_i^2,\dots ,{\tilde{\mathbf{z}}}_i^{N_d},{\tilde{\mathbf{z}}}_i^{1,2},\dots ,{\tilde{\mathbf{z}}}_i^{N_d-1,N_d}\right]\hfill \\ \hfill & \begin{array}{cc}\hfill \text{where}& {\tilde{\mathbf{z}}}_i^p={\tilde{\phi }}^p\left(x_i^p\right)\text{with}\hfill \\ \hfill & {\tilde{\phi }}^p\left(x_i^p\right)=\left[I\left(x_i^p<{\tau }_1^p\right),I\left({\tau }_1^p\le x_i^p<{\tau }_2^p\right),\dots ,I\left({\tau }_{N_p-1}^p\le x_i^p\right)\right]\hfill \\ \hfill \text{and}& {\tilde{\mathbf{z}}}_i^{p_1,p_2}=\text{vec}\left({\tilde{\phi }}^{p_{1,2}}\left(x_i^{p_1},x_i^{p_2}\right)\right)\text{with}\hfill \\ \hfill & {\tilde{\phi }}_{k,l}^{p_1,p_2}\left(x_i^{p_1},x_i^{p_2}\right)=I\left({\tau }_{k-1}^{p_1}\le x_i^{p_1}<{\tau }_k^{p_1}\&{\tau }_{l-1}^{p_2}\le x_i^{p_2}<{\tau }_l^{p_2}\right)\hfill \\ \hfill & k\in \left\{1,\dots ,N_{p_1}\right\},l\in \left\{1,\dots ,N_{p_2}\right\}\hfill \\ \hfill & \text{defining}{\tau }_0^{p_1}={\tau }_0^{p_2}=-\infty \text{and}{\tau }_{N_{p_1}}^{p_1}={\tau }_{N_{p_2}}^{p_1}=\infty \hfill \end{array}\hfill \end{array}}$ In this equation, I is an indicator function, yielding 0 or 1 when its argument is false or true, respectively. N_p refers to the number of intervals for variable x^p. Summarizing, every observation $x_i^p$ of a variable x^p is expanded into a binary vector ${\tilde{\mathbf{Z}}}_i^p$ containing only a single 1 indicating the interval to which the value of x^p belongs. Two-way interactions (x^p₁, x^p₂) yield a binary matrix containing only a single 1. The vec in the fourth line in the equation indicates vectorization of this matrix to obtain a vector. Finally, transformations for several variable values are simply concatenated to obtain one sparse vector ${\tilde{\mathbf{Z}}}_i$, a row of the matrix $\tilde{\mathbf{Z}}$. Which main effects and interactions are considered initially can be set by the user.

Reweighting. To obtain the simplification shown in Fig. 1, iterative reweighting is performed. It is a repetitive execution of the lpICS formulation in Eq. (1), where $\mathbf{D}\tilde{\mathbf{w}}$ has been replaced by $\tilde{\chi }\mathbf{D}\tilde{\mathbf{w}}\cdot \tilde{\chi }\in {\mathbb{R}}^{N_{ds}\times N_{ds}}$ is a diagonal matrix containing weights ${\tilde{\chi }}_{i,i}$, i = 1, …, N_ds. The reweighting values are computed based on the current solution of the problem $\tilde{\mathbf{w}}$ as ${\tilde{\chi }}_{i,i}=1∕\left({ϵ}_1+a\left|{\left(\mathbf{D}\tilde{\mathbf{w}}\right)}_i\right|\right)$. The reweighting can be summarized as:

 
_________________________________________________________________________________ 
Algorithm 2 Iterative reweighting in lpICS. 
_________________________________________________________________________________ 
  1:  Given:  ˜ w, a 
  2:  repeat 
  3:      ˜ wold = ˜ w 
  4:      ˜ χi,i =         1__________ 
ϵ1+a|(D˜wold)i|, ∀ i 
  5:      ˜ w ← solve Equation 1 with D˜w replaced by ˜ χ D˜w 
  6:  until AVG(|˜w − ˜ wold|) < ϵ2 or #iter == 10 
 __________________________________________________________________________________    

ϵ₁ and ϵ₂ are small values of 5e-4 and 1e-8, respectively. They help to avoid singularities and numerical issues. The value a can be selected by the user based on crossvalidation results on the training data. The loop continues until convergence, with a maximum of 10 iterations. Note also the difference with the loop in Algorithm 1 which surrounds the call to Algorthm 2. In the formercase, the aim was to remove variables, whereas here, the goal lies in the simplification of the model for every variable. Yet, in extreme cases, this simplification of the variable might also result in its removal.

Elastic net ICS (enICS)

The dimensionality of Z can be very high due to the binary expansion. As a consequence, the execution time is high, even though some solvers are able to take advantage of the sparse structure. To improve the execution time for larger problems, ICS supports a slightly different formulation based on elastic net, enICS. As in lpICS the sparsity is not induced on w, but on the difference Dw. However, this is done implicitly. Consider the formulation, in terms of the notation of lpICS: (3)${\displaystyle \underset{\tilde{\mathbf{w}},b}{min}{\parallel \tilde{\mathbf{Z}}\tilde{\mathbf{w}}+b-\mathbf{y}\parallel }_2^2+ϵ{\parallel \mathbf{D}\tilde{\mathbf{w}}\parallel }_2^2\text{s.t.}\parallel \mathbf{D}\tilde{\mathbf{w}}{\parallel }_1\le t.}$ In this equation, the symbols have the same meaning and dimensions as before. Additionally, ϵ is a parameter fixed at 0.05. It provides the stability of elastic net to the LASSO formulation, but other than that, it can be ignored. The hyperparameter t defines the tradeoff between sparsity (model simplicity) and performance. Note that the modelling is posed as regression rather than classification. Instead of minimizing the errors with respect to the classification border, it explicitly aims at the minimization of the deviation from the label values 1 and −1.

The advantage of using this approach becomes apparent when using a slightly different data representation converting it to a standard elastic net formulation, the actual enICS: (4)${\displaystyle \underset{\hat{\mathbf{w}},b}{min}{\parallel \hat{\mathbf{Z}}\hat{\mathbf{w}}+b-\mathbf{y}\parallel }_2^2+ϵ{\parallel \hat{\mathbf{w}}\parallel }_2^2\text{s.t.}\parallel \hat{\mathbf{w}}{\parallel }_1\le t.}$ Conceptually, $\hat{\mathbf{Z}}=\tilde{\mathbf{Z}}{\mathbf{D}}^{-1}$, although it is never calculated explicitly, but results from direct construction. Therefore, Eq. (4) directly optimizes the weight differences $\hat{\mathbf{w}}$. Zhou et al. (2015) have shown that the elastic net is equivalent to the dual formulation of a squared-hinge loss support vector machine. Due to its formulation, it can be efficiently solved in the primal (Chapelle, 2006) with a complexity related to the number of data points rather than the size of the expanded feature space, which is usually greater. It makes the problem complexity independent of the size of the feature space and allows for much faster execution. Furthermore, it automatically solves either the primal problem or the dual problem, depending on its characteristics e.g., in case of a large dataset. This approach is called SVEN (Support Vector Elastic Net).

Data transformation. It was mentioned that $\hat{\mathbf{Z}}$ is the result of direct construction from the data, similar to $\tilde{\mathbf{Z}}$. Hence, an enICS data transformation $\hat{\phi }$ can be defined as follows for a single observation x_i: (5)${\displaystyle \begin{array}{cc}\hfill & {\hat{\mathbf{Z}}}_i=\left[{\hat{\mathbf{Z}}}_i^1,{\hat{\mathbf{Z}}}_i^2,\dots ,{\hat{\mathbf{Z}}}_i^{N_d}\right]\hfill \\ \hfill & \text{where}{\hat{\mathbf{Z}}}_i^p={\hat{\phi }}^p\left(x_i^p\right)=\left[1,I\left(x_i^p\ge {\tau }_1^p\right),I\left(x_i^p\ge {\tau }_2^p\right),\dots ,I\left(x_i^p\ge {\tau }_{N_p-1}^p\right)\right]\hfill \end{array}.}$ Notice how the transformation uses the same bins as lpICS with the same bin thresholds τ. However, $\hat{\phi }$ uses more than a single 1 to encode a data point. If a value corresponds to a certain interval, all previous intervals are filled with 1. Hence, the weights no longer reflect the impact of an interval, but a difference with respect to the previous interval. Of course, in the end, one is still interested in the importance of the intervals themselves. This is easily calculated from the enICS weights by cumulative summing for each variable encoded in a reconstruction matrix R, such that $\tilde{\mathbf{w}}=\mathbf{R}\hat{\mathbf{w}}$.

Notice also how Eq. (5) no longer contains interaction terms. Indeed, SVEN is not applicable in that case since Eq. (3) becomes an arbitrary generalized elastic net. This can be converted to a standard elastic net, but additional constraints need to be added (Xu, Eis & Ramadge, 2013) which makes SVEN impossible. Therefore, enICS only supports main effects.

Reweighting. The reweighting procedure for enICS is similar to lpICS. Based on the current solution $\hat{\mathbf{w}}$ and the weighting parameter a supplied by the user, reweighting can be done as expressed in Algorithm 3. A key difference with Algorithm 2 lies in the calculation and application of the reweighting values $\hat{\chi }$ since they are applied on the data $\hat{\mathbf{Z}}$, rather than in the (lpICS) regularization on $\mathbf{D}\tilde{\mathbf{w}}$. Hence, they are the reciprocal of the lpICS reweighting values $\tilde{\chi }$, neglecting ϵ₁. The value of ϵ₂ remains at 1e-8.

 
____________________________________________________________________________________________ 
Algorithm 3 Iterative reweighting in enICS. 
____________________________________________________________________________________________ 
    1:  Given:  ˆ w, a 
    2:  repeat 
    3:      ˆ wold = ˆ w 
    4:      ˆ χi,i = a|ˆwi| 
    5:      ˆ w ← solve Equation 4 with ˆ Z replaced by ˆ Zˆ χ 
    6:  until AVG(|ˆw − ˆ wold|) < ϵ2 or #iter == 25 
____________________________________________________________________________________________    

Preselection

enICS was introduced as a way to improve the execution speed of ICS for larger datasets with a higher dimensionality. Instead of using another model, one can also perform feature selection beforehand (Billiet, Van Huffel & Van Belle, 2017). This reduces the dimensionality of the data. Yet, the aim is still to take ICS’s binning structure into account. Therefore, the so-called preselection is performed as a four-step procedure:

1. The data transformation in Eq. (2) is performed on the dataset X. Note that interaction can be included at this point as well. Using the computed $\tilde{\mathbf{Z}}$ and the known labels, one can solve a standard Support Vector Machine classification problem with linear kernel yielding coefficients ω.

2. A new data set X₂ is constructed. Main effects are modelled as ${\mathbf{X}}_2^p={\tilde{\mathbf{Z}}}^p{\omega }^p$, ∀p = 1, …, N_d. This represents the variable values in X by the coefficients of their intervals. Hence if only main effects are involved, X₂ ∈ ℝ^N×N_d has the same size as X. Each considered interaction increases the dimensionality of X₂ with one. These additional columns are given as ${\tilde{\mathbf{Z}}}^{p_1,p_2}{\omega }^{p_1,p_2}$, one for each interaction between variables x^p₁ and x^p₂. This approach can be considered as a data-driven discrete non-linear transformation. X₂ encodes data trends for each variable and interaction.

3. The goal remains to eliminate certain variables. To obtain this, any feature selection method could be applied on the trend-informed data X₂. Here, the elastic net was chosen. Similarly to the iterative reweighting discussed earlier, the selection is semi-automatic: the user can decide on the sparsity based on cross-validated performance estimation.

4. Finally, the elastic net indicates which columns of X₂ should be kept. Since every column corresponds to one main effect or interaction, these columns correspond to the set of effects and interactions that ICS should consider. This is an implicit input in the transformation function φ in the main ICS flow in Algorithm 1.

Setting up a problem

The toolbox groups all information about a specific ICS problem in a MATLAB struct problem. It has only two mandatory fields: xtrain, containing the training data and ytrain, the corresponding (binary) labels. The latter can contain any two numerical values, which will be transformed to the set {−1, 1}. Similarly, one can supply test data in the fields xtest and ytest. Checks on dimensionality correspondence are performed automatically, yielding error messages for incorrect input. Other fields that can be set include the variable types, their names, the groups (which main effects and/or interactions should be considered) and the algorithm to be used (lpICS or enICS). Furthermore, an options field allows to access more advanced settings such as whether the final model should be visualized, whether preselection should be performed, the maximum number of intervals for each variable, whether selection of the reweighting parameter should be automatic and if so, which threshold should be used to do so. The details of all the fields are discussed extensively in the ICS manual accessible with the toolbox at https://www.esat.kuleuven.be/stadius/icstoolbox.

Once a problem has been defined, the algorithm can be started by calling the buildICS function with the problem struct as argument. Snippet 1 shows a minimal and full example to get started with ICS. It assumes the existence of a file data.mat containing training matrices X_tr, y_tr and test matrices X_te, y_te. The matrices X_tr and X_te contain five variables (columns). Note that part of the code in the full setup is superfluous, since it corresponds to the defaults. It is only included as illustration.

 
 
                   Snippet 1: A minimal and full ICS setup example 
____________________________________________________________________________________ 
 
load  data . mat ; % contains  Xtr ,  Xte ,  ytr  and  yte 
 
% %%%%  minimal  setup  %%%%% 
%set  training  data  fields 
problem1 . xtrain  = Xtr ; 
problem1 . ytrain  = ytr ; 
 
% %%%%  full  setup  %%%%% 
%set  data  fields 
problem2 . xtrain  = Xtr ; 
problem2 . ytrain  = ytr ; 
problem2 . xtest  = Xte ; 
problem2 . ytest  = yte ; 
 
%the  numbers  refer  to  their  respective  columns  in  Xtr   and  Xte 
%x1   and  x3:  continuous,  x2   and  x4  :  categorical ,  x5  :  binary 
problem2 . var  = struct (` cont ',[1,3],` cat ',[2,4],` bin ',[5],` ord ',[]) ; 
problem2 . varnames  = {` foo ',` bar ',` rock ',` paper ',` scissors '}; 
%main  effects  for  variables  x1  ,  x2   and  x5  , %interactions  between  variables  x1   and  x3  ,  x2   and  x4 
problem2 . groups  = {1, 2, 5, [1,3], [2,4]}; 
%using  lpics 
problem2 . mainmethod  = ` lp '; 
 
%options 
problem2 . options . show  = true ; 
problem2 . options . preselect  = true ; 
%the  maximum  number  of  intervals  for  all  variables  is  set  to  42 
%in  practice ,  this  only  refers  to  continuous  and  ordinal  variables 
%categorical  and  binary  variables  revert  to  the  defaults . 
problem2 . options . maxbins  = 42; 
%manual  tuning  for  preselection  and  reweighting 
problem2 . options . auto  = false ; 
%suggested  weights  will  result  in  an  AUC  of  75%  of  the  maximum  cross 
%validation  AUC  for  preselection  and  80%  for  reweighting 
problem2 . options . cutoff  = [0.75, 0.8]; 
 
% %%%%  compute  the  ICS  model  for  both  problems  %%%%% 
ICS1  = buildICS ( problem1 ) ; 
ICS2  = buildICS ( problem2 ) 
____________________________________________________________________________________ 
 
  
    

Interactive weight selection

Users can interact with the training procedure at two points. Firstly, to make a sparsity-performance trade-off for preselection. Secondly, to balance the model complexity with regard to performance and risk calibration during the iterative reweighting. Only the latter is discussed here as they are very similar.

The interface for weight selection is shown in Fig. 2. It presents five graphs with the weight as independent variable. The results shown are based on cross-validation on the training data.

The top row shows performance measures. Figure 2A contains the Area under the ROC curve (AUC) as general performance characteristic. It equals 1 in case of perfect classification and 0.5 for guessing. The other two graphs in the top row represent a summary of the calibration curve relating the estimated risk with the empirical risk, namely as slope and bias (Fig. 2B , resp. Fig. 2C). In an ideal case, it has a unit slope and zero bias.

The bottom row indicates model complexity. Figure 2D lists the number of variables for each weight, whereas Fig. 2E displays the total number of intervals. Often, these two are closely related, but this might differ e.g., when interactions are present in the model.

Each graph shows a star representing the estimated optimal weight based on the user-defined cutoff. However, in this case it is clear that the model can be improved manually. The selection of the user is indicated with a green circle. Here, it indicates a model with a lower estimated complexity, higher AUC and better calibration values.

The user can select the preferred weight from the evaluated discrete set by means of a slider. A special value is −1. It corresponds to no reweighting.

Visualization

Visualization is an important factor. The toolbox supplies both a colour and grayscale representation. The visualization of the model itself consists firstly of a representation of all effects. An example for a main effect with three intervals is displayed in Fig. 3. The values of the thresholds τ are marked at the interval boundaries. The number of points for each interval is indicated by a colour code and by numerical values. In this case, if the value under consideration is higher than 0.91, two points will be added to the total score. Similarly, this can be done for other effects not displayed here. The same colour scheme is used for all effects, that is, the same colours represent the same values. In case of interactions, the representation is a grid rather than a row.

A second element of the model is the colour-coded risk profile, as shown in Fig. 4. It maps the scores (sum of the points) on top to estimated risks at the bottom. For example, if the sum of the contributions of all effects in the model yield a final score of three, the risk of belonging to the target class would be 75%.

The toolbox provides a listing of all effects and options for its visualization by means of a graphical interface (Fig. 5). It is automatically called at the end of model training (unless disabled), but can also be called on an already trained model with the function showICS. The Figure shows the model for the ionosphere dataset included with MATLAB. Complete code to generate it is available as a demo file in the toolbox. The two buttons on top display the assessment module (discussed below) and the risk profile. The middle part displays all selected effects. Here, three effects are part of the final model, as can be seen in the listbox. To its right, one can indicate if numerical values should be overlaid and whether to use colour coding or grayscale (e.g., for printing). The model details for one or more effects can be displayed via the buttons. Finally, the bottom part contains the common colour legend for all effects.

Assessment

The assessment module can be launched from the visualization interface (see Fig. 5). Normally, it shows assessment results for the training data only. If test data has been included at model setup, test results are summarized as well. Additionally, a model can be applied on new data with the function applyICS. Its outputs can be given as additional inputs to showICS to include results on the new data in the assessment.

The assessment window is displayed in Fig. 6. It evaluates the model at two levels. Firstly, the two graphs at the top show an alternative representation of the risk profile (Fig. 6A) and the calibration curve (Fig. 6B). The latter evaluates how well the predicted risk matches the empirical risk. Numerically, this is summarized in the lower right quadrant with the slope and bias of the curve, for training and test data. Secondly, the classification performance is evaluated using the Receiver Operating Characteristic (ROC) curve (Fig. 6C). It shows the tradeoff between sensitivity and specificity with an indication of the chosen tradeoff point. Moreover, the Area Under the Curve (AUC) is a robust indication of performance. Hence, it is included in the summary in the lower right quadrant, together with the classification accuracy (percentage of correctly classified samples). Finally, the number of intervals is included as a measure for the model complexity. The combination of these measures allows to thoroughly assess the quality of the obtained model.