Computational Methods for the Analysis of Complex Data

Programa de doctorado “Matemáticas” PHD DISSERTATION COMPUTATIONAL METHODS FOR THE ANALYSIS OF COMPLEX DATA Author MaRemedios Sillero Denamiel Supervisors Prof. Dr. Rafael Blanquero Bravo Prof. Dr. Emilio Carrizosa Priego Prof. Dra. Pepa Ramírez Cobo some text to complete some text to complete some text to complete some text to complete some text to complete some text to complete some text to complete some text to complete some text to complete A mis padres y mi hermana. A mis abuelos. Agradecimientos Y al final todo llega. Qué gran verdad. Siempre había visto este momento como algo muy lejano, intentando disfrutar al máximo del camino y ya, hoy, se cierra esta maravillosa etapa. Recuerdo perfectamente el día en el que mi director Emilio me brindó la oportunidad de trabajar en un proyecto de transferencia, junto con mi también director Rafa, y mi compañera y amiga Asun. Esto cambió todos mis planes, y de qué manera. Gracias por la confianza que depositaste en mí. Emilio, Rafa, Pepa, qué suerte teneros como directores. Me iniciasteis en la investigación, un camino que requiere mucha dedicación y constancia, pero que gratifica más aún. Gracias por enseñarme tanto, por vuestro tiempo y apoyo durante estos años. No siempre ha sido fácil, pero vosotros me habéis ayudado a continuar. Emilio, gracias por transmitirme tu pasión por tu trabajo sin darte cuenta, por tu dedicación y por animarme siempre que lo he necesitado (y siempre que no, también). Rafa, gracias por tu paciencia infinita cuando las cosas no salían, encontrando una solución para todo y siempre con una sonrisa. Pepa, has estado al pie del cañón, y así me lo has hecho saber. Hemos compartido mucho más que esta tesis, eres espectacular. Gracias por siempre estar. Las personas con las que trabajo son, sin duda, lo mejor de este trabajo. Poder viajar y vivir en otros países mientras aprendes de grandes profesionales es un regalo. Loli, cómo me gusta Copenhague, tu segunda ciudad. Gracias por hacerme sentir como en casa durante mis estancias, aun estando a tres mil kilómetros de mi familia. Todavía nos quedan muchos cafés juntas, y muchas exposiciones por descubrir. Mike, thank you for always welcoming me with a smile, for your patience teaching me and for the moments of chatting in the canteen while having coffee and tea. También me gustaría agradecer al Departamento de Estadística e Investigación Operativa de la Universidad de Sevilla, al IMUS y a sus equipos de Administración, por facilitarme siempre todo lo que he necesitado. Y, por supuesto, a los proyectos y a los IP de los proyectos que han financiado esta tesis, MTM2015-65915-R (Ministerio de Economía y Competitividad), PID2019-110886RB-I00 (Ministerio de Ciencia, Innovación y Universidades), FQM-329 y P18-FR-2369 (Junta de Andalucía), PR2019-029 (Universidad de Cádiz), Fundación BBVA y EC H2020 MSCA RISE NeEDS Project (Grant agreement ID: 822214). Cuántos momentos compartidos con mis compañeros del IMUS. Más que compañeros, IV V amigos. No los he podido tener mejores. Habéis sido cómplices en los buenos momentos y un apoyo fundamental en los no tan buenos. ¡Me siento tan orgullosa de formar parte de esta familia! No sé dónde acabaré, pero le habéis dejado el listón demasiado alto a mis futuros compañeros. Y, aunque suene a despedida, esto no queda aquí, aún nos queda mucho por disfrutar juntos. Gracias por tanto. Mis amigas de toda la vida, las que siempre están. Lorena, Eva, Virginia, sois lo más. No tengo palabras de agradecimiento a tantos años de amistad de la buena. Daniela, qué suerte que nos cruzáramos en el camino, eres única y no sabes cuánto nos enseñas. Marisa, gracias por hacer las veces de hermana mayor y por todos tus consejos. Vicente, al final aquellas tardes de estudio han valido la pena. Siempre habéis creído en mí, incluso mucho más que yo misma. Esta tesis también es vuestra. Mamá, papá, qué os voy a decir. No hay día que no me mostréis lo orgullosos que estáis de mí, pero yo sí que estoy orgullosa de vosotros, y de la familia que somos. Todo lo que he conseguido y consiga es y será gracias a vosotros y a vuestro apoyo incondicional. Rocío, mi hermana y confidente, siempre contaremos la una con la otra. Abuela Aurelia, eres fortaleza y un pilar fundamental para mí. Abuelo Higinio, abuela Francisca, abuelo Andrés, sé que desde arriba estáis celebrándolo conmigo. Os quiero. Gracias a todos, de corazón. Resumen Esta tesis combina las disciplinas de Investigación Operativa y Estadística con el fin de desarrollar nuevos métodos computacionales para extraer información de datos complejos. En este estudio, datos complejos se refiere a conjuntos de datos con un número elevado de muestras y/o variables, con diferentes tipos de variables, con estructuras de dependencia entre las variables, recogidos de diferentes fuentes (heterogéneos), posiblemente con clases desbalanceadas, con diferentes costes de clasificación incorrecta o caracterizados por valores extremos (datos de cola pesada), entre otros. La complejidad de los datos y las nuevas exigencias de los usuarios (modelos interpretables, modelos sensibles a costes de errores de predicción o modelos eficientes en tiempo de ejecución) implica un reto desde la perspectiva científica. Las principales contribuciones de esta tesis se engloban en tres marcos teóricos diferentes: Regresión, Clasificación e Inferencia Bayesiana. Repecto al primero, consideramos modelos de regresión lineal, donde una variable respuesta continua se pronostica a partir de un conjunto de variables predictoras. Por un lado, buscando soluciones interpretables en datos heterogéneos, proponemos una nueva versión del Lasso en la que se controla el rendimiento del modelo en los grupos de interés. Por otro lado, aplicamos técnicas de optimización matemática para proponer un modelo de regresión lineal diseñado específicamente para conjuntos de datos con variables predictoras categóricas y jerárquicas. En lo que se refiere a Clasificación, en esta tesis se ha explorado en profundidad el clasificador Naïve Bayes. Este método se ha adaptado para obtener una solución sparse (es decir, expreseada en términos de un subconjunto de variables predictoras). Además, el método se ha modificado para tratar con datos en los que se debe tener en cuenta los diferentes costes de clasificación incorrecta. En ambos problemas, se presentan nuevas estrategias para reducir los tiempos de ejecución. Finalmente, la última contribución de esta tesis es relativa a la inferencia Bayesiana. En particular, se considera un enfoque de estimación Bayesiano semiparamétrico para estimar la distribución Elíptica, la cual generaliza la distribución normal multivariante permitiendo colas más pesadas. La estructura de esta tesis es la siguiente. El Capítulo 1 revisa los conceptos teóricos necesarios para desarrollar los capítulos posteriores. En concreto, se revisan dos áreas de investigación principales: el aprendizaje automático (sparse y sensible a costes de errores de predicción) y la estadística Bayesiana. VI VII En el Capítulo 2 se propone un método basado en el Lasso en el que se han añadido restricciones cuadráticas para controlar los errores de predicción en individuos de interés. Este modelo de regresión, sparse y con restricciones, se define mediante un problema de optimización no lineal. El método resulta de especial interés cuando se trabaja con muestras heterogéneas donde los datos provienen de diferentes fuentes, como resulta habitual en muchos contextos biomédicos. El Capítulo 3 estudia modelos de regresión para variables predictoras categóricas con estructura jerárquica. El modelo es flexible en el sentido en que el usuario decide el nivel de detalle en la información que el modelo debe utilizar, atendiendo a razones de privacidad y confidencialidad. Para modelar el equilibrio entre el rendimiento del modelo y su complejidad, se define un problema cuadrático convexo con variables enteras y restricciones lineales. En el Capítulo 4, se introduce una versión sparse del clasificador Naïve Bayes, caracterizada por las tres propiedades siguientes. Primero, la selección de variables se realiza teniendo en cuenta la estructura de correlación de las variables predictoras. Segundo, se pueden usar diferentes medidas de rendimiento al seleccionar los subconjuntos de variables. Finalmente, el modelo permite incluir restricciones de rendimiento en los grupos de individuos de mayor interés. Mediante el diseño de una búsqueda inteligente, el método consigue tiempos de ejecución competitivos. El enfoque introducido en el Capítulo 2 también se explora en el Capítulo 5 con el fin de mejorar la predicción del clasificador Naïve Bayes en las clases de mayor interés para el usuario. A diferencia de la versión tradicional del clasificador basado en dos etapas (estimación primero y clasificación después), la nueva metodología integra ambas a la vez. El método se formula como un problema de optimización donde se maximiza la función de verosimilitud con restricciones para las tasas de clasificación de los grupos de interés. El Capítulo 6 aborda la estimación estadística desde la perspectiva Bayesiana de una clase general de distribuciones de cola pesada, las distribuciones Elípticas. Esta familia de distribuciones se han propuesto a nivel teórico en problemas asociados a datos extremos, típicos en contextos financieros o en teoría de riesgo. El enfoque que adoptamos es semiparamétrico basado en los procesos Dirichlet. Finalmente, el Capítulo 7 cierra esta tesis con conclusiones generales y futuras líneas de investigación. Abstract This PhD dissertation bridges the disciplines of Operations Research and Statistics to develop novel computational methods for the extraction of knowledge from complex data. In this research, complex data stands for datasets with many instances and/or variables, with different types of variables, with dependence structures among the variables, collected from different sources (heterogeneous), possibly with non-identical population class sizes, with different misclassification costs, or characterized by extreme instances (heavy-tailed data), among others. Recently, the complexity of the raw data in addition to new requests posed by practitioners (interpretable models, cost-sensitive models or models which are efficient in terms of running times) entail a challenge from a scientific perspective. The main contributions of this PhD dissertation are encompassed in three different research frameworks: Regression, Classification and Bayesian inference. Concerning the first, we consider linear regression models, where a continuous outcome variable is to be predicted by a set of features. On the one hand, seeking for interpretable solutions in heterogeneous datasets, we propose a novel version of the Lasso in which the performance of the method on groups of interest is controlled. On the other hand, we use mathematical optimization tools to propose a sparse linear regression model (that is, a model whose solution only depends on a subset of predictors) specifically designed for datasets with categorical and hierarchical features. Regarding the task of Classification, in this PhD dissertation we have explored in depth the Naïve Bayes classifier. This method has been adapted to obtain a sparse solution and also, it has been modified to deal with cost-sensitive datasets. For both problems, novel strategies for reducing high running times are presented. Finally, the last contribution of this dissertation concerns Bayesian inference methods. In particular, in the setting of heavy-tailed data, we consider a semi-parametric Bayesian approach to estimate the Elliptical distribution. The structure of this dissertation is as follows. Chapter 1 contains the theoretical background needed to develop the following chapters. In particular, two main research areas are reviewed: sparse and cost-sensitive statistical learning and Bayesian Statistics. Chapter 2 proposes a Lasso-based method in which quadratic performance constraints to bound the prediction errors in the individuals of interest are added to Lasso-based objective functions. This constrained sparse regression model is defined by a nonlinear optimization problem. Specifically, it has a direct application in heterogeneous samples where data are VIII 4 Chapter 1. Introduction The huge assorment of datasets describing timely and relevant real-world situations requires revisiting and updating decision methods and combine them with data analysis techniques, to yield Data-Driven Decision Making. Contemporary data are characterized by many instances and/or variables, different types of variables (quantitative, categorical, ordinal, clustered), dependence structures among the variables or extreme values (heavy-tailed data). Also, datasets can be characterized by asymmetric conditions (non-identical population classes size, different misclassification costs), where different errors may have different consequences (as happens in the context of medical diagnosis). Complexities can also be caused by the hierarchical structure of the data when, for instance, national statistics need to work with the partition of the country into regions, provinces, municipalities or neighbourhoods; or when a retailing company aims to group their products to be able to predict at the aggregate level. Recently, the complexity of the raw data in addition to new requests posed by practitioners (interpretable models, fair models or models which are efficient in terms of running times or memory demand for prediction) entail a challenge from a scientific and technological perspective. The aim of this thesis is to propose novel computational methods for the extraction of knowledge from the above-mentioned complex data. In fact, this thesis bridges the disciplines of Operations Research and Statistics to lead to new approaches that outperform current methods. 1.1 Sparse and cost-sensitive statistical learning The challenge of interpreting information and learning from complex data is at the core of the current Statistical Science. The aim of Statistical Learning theory is the estimation of a function f(X)for predicting the response variable Ygiven the set of predictors X, a random vector of dimension p. We refer to regression problem when the response Yis quantitative, whereas the term classification is employed in the case of a categorical response. 1.1.1 Sparsity in regression and classification Current datasets are usually characterized by a large number of features (that is, the dimension of the feature space pis large). This fact may have negative consequences in terms of the comprehensibility of solutions and, thus, Statistics and Operations Research fields are continually adapting to tackle this matter [Friedman et al., 2001; Hastie et al., 2015]. For example, a clear interpretation of the solutions is of crucial importance for some specific medical or credit scoring related problems, where the interest is to find the key variables that determine if a patient is sick or healthy, or those for predicting a bank’s customer as potentially defaulting, see Hand and Henley [1997]. The search for more interpretable and parsimonious solutions, common in multivariate contexts such as regression [Cai et al., 2009; Lin et al., 2011; Benítez-Peña et al., 2021], clustering [Maldonado et al., 2015], time series analysis [Carrizosa et al., 2017b] 1.1. Sparse and cost-sensitive statistical learning 5 or visualization [Carrizosa and Guerrero, 2014], has recently led to the development of sparse multivariate techniques, see Hastie et al. [2015]. Excessive computational costs, redundant variables and noise are other drawbacks associated with high-dimensional data. Besides the identification of significant predictors, which provides a good interpretation of the model reducing the computational costs and noise, another fundamental criterion for evaluating the model’s performance is the accuracy of prediction. It could be said that looking for sparse models helps to obtain low prediction errors, since it avoids overparametrized models and thus, overfitting [Carrizosa et al., 2016]. For those reasons, sparsity is a desirable property that regression and classification models should satisfy. Regarding parametric regression contexts, a solution is said to be sparse if only a subset of coefficients are non-zero. In particular, in some parts of this dissertation, we will focus on linear regression models, where the response variable is expressed as a linear combination of the predictors. The conventional linear regression procedure to estimate the coefficients, the Ordinary Least Squares (OLS), is known not to be sparse and, as a consequence, new penalization techniques appeared in the literature, see Gui et al. [2017] and Li et al. [2020] for a detailed overview. An example is the best-subset selection [Garside, 1965], which achieves sparser solutions, but suffers from high variability and computational difficulties [Fan and Li, 2001]. In contrast, the ridge regression [Hoerl and Kennard, 1970], a continuous shrinkage method that adds to the objective function an l2−norm penalty over the coefficients to be estimated, achieves its better accuracy prediction through a bias-variance tradeoff. Nevertheless, ridge regression is known not to be able to render a parsimonious solution. To overcome those shortcomings, Tibshirani [1996] proposed the Lasso regularization technique, which includes an l1-norm penalization term instead and, thus, achieves both estimation and selection of relevant predictors simultaneously by construction. One of the advantages of the Lasso is that the entire path of solutions can be found thanks to the LARS algorithm [Efron et al., 2004]. In addition, it is well-known that, under some conditions, the Lasso enjoys good theoretical and statistical properties [Donoho et al., 1995; Friedman et al., 2001; Bühlmann and Van-De Geer, 2011]. To visualize the effect of the penalty term in the Lasso formulation, consider the wellknown prostate database [Stamey et al., 1989], which consists of the measurements of 8 predictors and one response variable (clinical measures) on 97 men who were about to receive a radical prostatectomy. If the goal is to minimize the overall mean squared error (MSE), the parameter vector can be estimated by a fitting procedure as OLS. The results obtained under the OLS are shown in the first two rows of Table 1.1, where 3/4of the total set has been used to fit the model (training set) and the remaining samples (testing set) to the assessment of the generalization error of the resulting model. The overall MSE, as well as the number of coefficients involved in the model, are presented. The third and fourth rows in Table 1.1 provide the results obtained under the Lasso for the prostate dataset. In this case, in comparison with OLS results, a sparser and therefore, a more interpretable solution has been obtained at 6 Chapter 1. Introduction Method Overall MSE Non-Zero Coefficients OLS Training set 0.344 8 Testing set 0.373 Lasso Training set 0.365 5 Testing set 0.408 Table 1.1: Results obtained using prostate dataset the expense of slightly worsening the MSE values. Many different variants of the Lasso have been proposed. For example, in Zou [2006] adaptive weights for penalizing different coefficients in the l1penalty are included as a way for fitting sparser models under more general conditions. Moreover, in the presence of highly correlated predictor variables (as is usual in microarray studies) or when predictors are structurally grouped (e.g. dummy variables), the Lasso sometimes does not perform well and, as a consequence, the elastic net of Zou and Hastie [2005] and the group lasso [Yuan and Lin, 2006; Simon et al., 2011] were proposed. They combine l2and l1penalties to try to select (or remove) the correlated or structured predictor variables together. See Hastie et al. [2015] for an extensive review about the Lasso problem and generalizations. In classification, sparsity is closely linked to the concepts of Variable Selection and Feature Selection [George and McCulloch, 1993; Zou and Hastie, 2005; Lin et al., 2011; Carrizosa et al., 2016], whose aim is to identify the relevant variables within a set of many predictors so that classification accuracy (the percentage of class labels predicted correctly) is not reduced. First works published on feature selection handled datasets of a few predictors [Blum and Langley, 1997], but the size of data got larger and a number of variable selection techniques were proposed. There exist two main groups of methods that select features: filters [Guyon et al., 2006; Saeys et al., 2007] and wrappers [Kohavi and John, 1997; Saeys et al., 2007]. Broadly speaking, whereas the selection of features in the former group is based on the intrinsic properties of the features in addition to the correlation with the target variable, the latter tries to find a subset of features for training the model, adding or removing features from the subset according to the results drawn from the model. For instance, the filter Correlation based Feature Selection (CFS), see Hall [2000], is based on the assumption that a good subset of attributes should be highly correlated with the response variable but, on the other hand, there should exhibit low dependency among them. Other example can be the wrapper Boruta [Kursa and Rudnicki, 2010], which is in principle designed using a Random Forest strategy [Breiman, 2001], but can be modified and adapted to any classifier. There also exists an alternative to the use of the filters and wrappers, which has been recently considered in the literature. These are the so-called embedded methods, which embed the feature selection process into the classifier construction. In particular, there are some specific examples of embedded methods proposed for achieving sparse versions of the well-known Support Vector Machine (SVM) classifier [Cortes and Vapnik, 1995], see Aytug [2015], Maldonado et al. [2017], Ghaddar and Naoum-Sawaya [2018], 1.1. Sparse and cost-sensitive statistical learning 7 Blanquero et al. [2019] and Benítez-Peña et al. [2019], whereas Blanquero et al. [2020] look for interpretability for random forests. Concerning the Naïve Bayes (NB) classifier [Hand and Yu, 2001], some works have addressed different strategies for variable reduction. For example, McCallum and Nigam [1998] and Feng et al. [2015] base their feature selection approaches on the univariate correlations between features and the class. In this sense, Tang et al. [2016b] and Tang et al. [2016a] aim to rank the features according to their capacity for classification or a specific feature selection criterion, respectively. See Chandrashekar and Sahin [2014] for a survey of the differences among well-known filter, wrapper and embedded methods. 1.1.2 Cost-sensitive procedures Most traditional classification approaches search for maximizing accuracy, and do not take into account the differences between types of misclassification errors. However, in many real-world problems, such as those mentioned previously (medical diagnosis, credit card fraud detection), it is more important to achieve better classification rates for the individuals of interest (ill people, defaulting customers), since the consequences of wrong predictions across the classes may be very different. Cost-sensitive learning methods [Elkan, 2001; Zadrozny and Elkan, 2001] take into consideration different cost matrices that describe the costs for misclassification [Turney, 2000]. These approaches turn out to be very convenient for unbalanced datasets, where the minority class may be the worst classified (and the most critical one). In the case of binary classification (positive (1) and negative (0) classes), the cost matrix has the structure in Table 1.2. Whereas c10 is the cost of a false positive (actual negative but predicted as positive), c01 is the cost of a false negative. In contrast, c00 and c11 are associated with correct predictions. It seems natural that the cost of labeling a new individual incorrectly should always be greater than the cost of labeling it correctly. actual negative actual positive predict negative c00 c01 predict positive c10 c11 Table 1.2: Cost matrix for binary classification In classification contexts, Yidentifies the class. Then, the observed dataset is formed by pairs (y, x)from an unknown joint distribution. Given the cost matrix, the classifier will assign a new individual xto the class that has the minimum expected cost. Thus, if p(j|x)is the probability of classifying xinto class j, a new individual xwill be assigned to class 1 if and only if (iif) p(0 |x)c10 +p(1 |x)c11 ≤p(0 |x)c00 +p(1 |x)c01, 8 Chapter 1. Introduction or equivalently, p(0 |x)(c10 −c00)≤p(1 |x)(c01 −c11). As can be observed, the optimal decision is unchanged if a constant is added to each column of the cost matrix, which implies a simpler matrix where c00 and c11 are equal to 0. Finally, as p(0 |x)=1−p(1 |x), the classifier will assign the positive class to the new individual iif p(1 |x)≥c10 c10 +c01 .(1.1) Likewise, this procedure can be extended when more than two classes are considered, where more rows and columns are adding to the cost matrix. The classification threshold in (1.1) can be used by those classifiers which produce probability estimates, thus becoming cost-sensitive. Otherwise, misclassification costs can be directly introduced into the classifier construction. See Bradford et al. [1998], Freitas et al. [2007], Carrizosa et al. [2008], Sun et al. [2009], Datta and Das [2015] and Lee et al. [2017] for more details and applications. As some examples, consider Datta and Das [2015], Carrizosa et al. [2008] and Lee et al. [2017], which focus on the SVM classifier. In Datta and Das [2015] the decision boundary shift is combined with unequal misclassification penalties. On the other hand, in Carrizosa et al. [2008] a biobjective problem, which simultaneous minimizes the misclassification rates, is performed. In Lee et al. [2017], the authors propose a new weight adjustment factor that is applied to a weighted SVM. In the context of decision trees, Freitas et al. [2007]; Ling et al. [2004] introduce tree-building strategies which choose the splitting criterion by minimizing the misclassification costs, whereas Bradford et al. [1998] performs the pruning of a subtree following the cost information. Cost-sensitive versions of neural networks for unbalanced data classification have also been studied in the literature [Cao et al., 2013; Zhi-Hua Zhou and Xu-Ying Liu, 2006]. Other approaches can be found, for example in Peng et al. [2014], where a new version of the so-called data gravitation-based classification model is proposed. Although the term cost-sensitive learning has been mainly exploited in classification contexts [He and Ma, 2013; Prati et al., 2015], it can also be extended to regression contexts, where the response variable is quantitative but the whole sample is splitted into different groups (or classes). In particular, one can think of a Lasso model whose objective function includes weighted penalties over the MSE of each group (quadratic penalties). Those penalties would play the role of the costs introduced in Table 1.2. In fact, new versions of cost-sensitive regression models have been recently proposed in the literature (see Ollier and Viallon [2017] and reference therein). Under the umbrella of cost sensitive learning, it is generally assumed that misclassification costs are given and known. Unfortunately, fixing precise values for such misclassification costs may be problematic in real-world applications. In addition, in this way only an indirect control on misclassification rates is obtained. The application of mathematical optimization tools, the 1.2. The Bayesian Paradigm 9 approach undertaken in this thesis, seems to be promising [Carrizosa and Romero Morales, 2013] and not fully explored: one overall criterion is to be optimized, while constraints are introduced in the model to demand admissible values for the efficiency measures under consideration. 1.2 The Bayesian Paradigm Different disciplines, such as Statistics, are required to obtain theoretical results for the effectiveness of sparse and cost-sensitive statistical learning methods. Bayesian decision theory [Berger, 2013], involves an assortment of inference techniques that can be used for estimating probabilities of interest associated with the properties of learning algorithms. This section is devoted to introduce the Statistical Learning through a Bayesian decision theory perspective. 1.2.1 Bayesian Statistical Learning To select the best function fto predict the response Yin terms of the predictor variables X, the loss function L(f(X), Y )is used. To characterize the performance of f, the expected value of L(f(X), Y ), which is called the risk, is calculated. Given the observed sample X=x, the expression R(f) = EX[EY|X(L(f(X), Y )|X)], is known as the posterior expected loss or risk. Then, the best function fis that that minimizes the posterior expected loss. Such solution is known as the Bayes estimator under the loss function L. In the case of quantitative Y, the most common choice for Lis the squared error loss, L(f(X), Y ) = (f(X)−Y)2, which results in the regression function f(X) = E(Y| X=x). Otherwise, when dealing with categorical output Y, a different loss function for penalizing prediction errors is needed. If, instead of a regression problem, we are interested in a classification problem with K possible classes in C, then a matrix L∈MK(R)is introduced. Whereas its diagonal will be zero, extra-diagonal values, which represent the price paid for wrong predictions, will usually be equal to one. Then, the risk is defined as R(f) = EX"K X k=1 L[Ck,ˆ C(X)]p(Ck|X)#, where Ckmeans the class k,k∈ {1, . . . , K}, and the estimate ˆ Cassumes values in C. Under the 0-1 loss function choice, ˆ C(x) = arg max c∈C p(c|X=x), 10 Chapter 1. Introduction or equivalently ˆ C(x) = Ckif p(Ck|X=x) = max c∈C p(c|X=x). This method, which classifies each instance xin the most probable class via the conditional distribution, is known as the Bayes classifier. The computation of p(Ck|x)may be cumbersome if the number of features pis large. However, the use of the Bayes theorem eases the previous computation since p(Ck|x) = π(Ck)p(x|Ck) p(x), where π(Ck)is the prior distribution for the class, p(x|Ck)is the likelihood function of the data and p(x)is the so-called evidence. Since the evidence is the same for all the classes, in practice, the interest is in computing the numerator. Among the assortment of current classification techniques, the NB classifier has played a prominent role because of its simplicity, tractability and efficiency, see Hand and Yu [2001]. The method is based on the assumption of conditional independence of the features to the class, p(x|Ck) = p(x1, . . . , xp|Ck) = p(x1|Ck). . . p(xp|Ck),(1.2) which notably simplifies the computation of the class probability. One of the advantages of the NB is that it usually estimates fewer parameters than other renowned classifiers, so it is less prone to make overfitting [Domingos and Pazzani, 1997; Hand and Yu, 2001]. As a consequence, a number of applications of the NB in real contexts can be found, for example, in medicine [Wolfson et al., 2015], genetics [Minnier et al., 2015], reliability [Turhan and Bener, 2009], risk [Minnier et al., 2015] or document analysis [Guan et al., 2014], among others. 1.2.2 Bayesian Inference As commented at the beginning of this section, Statistical Learning methods can be viewed from a Bayesian paradigm. Furthemore, when dealing with datasets of especial characteristics (for example, heavy tails in contexts as Economics and Finance), Bayesian statistical techniques have shown their potential in the literature [Owen and Rabinovitch, 1983; Andrews et al., 1993; Ramírez-Cobo et al., 2010; Fortunati et al., 2020]. The Bayesian inference scheme can be briefly described as follows. Let Xbe the random variable whose distribution depends on the unknown parameter, θ∈Θ, to be estimated. Given a sample x, the model or density function is denoted by f(x|θ), which, as a function of θ, L(θ), is called the likelihood function. In the Bayesian paradigm, the parameter θis treated as a random variable, and thus is associated with a (prior) distribution function. Once the sample is observed, the likelihood function helps to update the prior distribution into the posterior distribution, which is the final target of the inference process. Two main advantages of performing Bayesian data analysis versus the classical (frequentist) framework should be pointed 1.2. The Bayesian Paradigm 11 out. First, the available prior information about θcan be coherently incorporated into the statistical model, that is, experts can input prior knowledge in the process of modeling. Second, given the observed data, the posterior distribution, encompasses a more extensive information than that of a puntual estimation made in frequentist way. Fundamentals As commented before, the unknown parameter θis considered a random variable with prior distribution π(θ), which describes uncertainty related to this unknown parameter before data are observed. In addition, in the case the prior distribution is specified up to other parameters, they are called the hyperparameters. Once the data have been observed (X=x) the posterior distribution for θis denoted by π(θ|x). Thanks to Bayes rule, π(θ|x)is calculated dividing the joint distribution by the marginal distribution, m(x), that is, π(θ|x) = f(x|θ)π(θ) m(x)=f(x|θ)π(θ) RΘf(x|θ)π(θ)dθ.(1.3) If the posterior distribution remain in the same family as the prior distribution, we will say that we have a conjugate structure. In fact, computationally speaking, this case is the most convenient. However, since simple conjugate analysis has a limited modeling capability of real-life data, to achieve a closed form for the posterior is not a rule. Therefore, sophisticated techniques in the literature, such as Markov Chain Monte Carlo (MCMC) methods, are required for Bayesian computation, see[Gelfand and Smith, 1990] and the following section for further details. Posterior simulation The two main difficulties of Bayesian inference are the prior elicitation and the computation of the posterior distribution (1.3). However, in practice, it does not have a closed form and, therefore, advanced techniques are required. As commented previously, MCMC methods allow for sampling from an unknown distribution. Robert and Casella [2013] gives the following definition: AMarkov Chain Monte Carlo (MCMC) method for the simulation of a distribution fis any method producing an ergodic Markov chain (X(t)) whose stationary distribution is f. Therefore, these methods are based on the construction of a Markov chain, say θ(j), whose stationary distribution is equal to the (posterior) distribution of interest. One of the most popular MCMC algorithms is the Metropolis-Hastings algorithm. Assume that the unknown parameter is indeed a multivariate parameter, θ= (θ1, . . . , θp). Therefore, 12 Chapter 1. Introduction the objective is to find the posterior distribution of θ. Given a proposal or conditional distribution q(˜ θ|θ(j)), which generates a candidate at iteration j+ 1 based on the previous accepted value θ(j), the candidate is accepted with certain probability. A summary of the Metropolis- Hastings algorithm is described in Algorithm 1. Algorithm 1: Metropolis-Hastings 1. Set initial values θ(0) = (θ(0) 1, . . . , θ(0) p). 2. j= 1 3. Generate ˜ θ∼q(θ|θ(j)). 4. Calculate ρ=ρ(θ(j),˜ θ) = min (π(˜ θ)q(θ(j)|˜ θ) π(θ(j))q(˜ θ|θ(j)),1) 5. Set θ(j+1) =˜ θwith probability ρ θ(j)with probability 1−ρ 6. Increase j=j+ 1 and return to 3. This algorithm always accepts values with a higher likelihood ratio than the previous value. Note that for a given target distribution π, the proposal qmust satisfy that πis contained in the union of supports of all conditional distributions q(.|θ). Necessary conditions on the proposal distribution q, to ensure the convergence to the limiting distribution of the chain (π), can be studied in depth in Chapter 6 of Robert and Casella [2013]. The Gibbs sampler [Gelfand, 2000], which is used for obtaining samples from a joint density function, is a special case of Metropolis-Hastings algorithm when the proposal distribution qis fixed as the conditional distribution f(θi|θ−i), where θ−i= (θ1, . . . , θi−1, θi+1, . . . , θp), that are assumed to be known. Algorithm 2 describes how Gibbs sampling works. Algorithm 2: Gibbs sampler 1. Set initial values θ(0) = (θ(0) 1, . . . , θ(0) p). 2. j= 1 3. Generate θ(j+1) 1∼θ1|θ(j) −1. Update θ(j). 4. Generate θ(j+1) 2∼θ2|θ(j) −2. Update θ(j). 5. . . . 6. Generate θ(j+1) p∼θp|θ(j) −p. Update θ(j). 7. j=j+ 1. Go to 3. Obviously, at each step j,ρ(θ(j) i,θ(j) −i,˜ θ)=1, which implies that each update in the Gibbs algorithm is accepted. Discussion about how to fix updating order can be found in Gilks et al. [1994]. In order to solve complex problems, a mixture of different MCMC algorithms (hibrid methods) are usually implemented. For instance, Müller [1991] proposes using Gibbs sampler 2.1. Introduction 19 the predictor matrix; then, the Lasso solution can be defined as ˆ βLasso(λ) = arg min β 1 nky−Xβk2+λk(β1, . . . , βp)k1(2.1) where λ≥0is a tuning parameter and k.k1is the l1norm. To visualize the effect of the penalty term in the Lasso formulation, consider the third and fourth rows in Table 2.1, which provide the results obtained under the Lasso for the prostate dataset. However, the Lasso presents some limitations; in particular, the literature related to the Lasso has not undertaken the problem of fully controlling the accuracy prediction on certain individuals of interest. In the previous prostate database, assume for instance that we are interested in fitting a sparse regression model to the dataset where, apart from obtaining a small overall mean squared error, also the prediction error for the young individuals should not exceed a given threshold. In this work we propose a Lasso-based model that allows for such aim, namely the cost-sensitive constrained Lasso, denoted from now on as CSCLasso. The results obtained for the prostate database under the CSCLasso, whose definition and main properties shall be discussed in Section 2.2 and 2.3, are shown in the last two rows of Table 2.1. A threshold for the mean squared error over Group 1 is set equal to 0.357, which represents an improvement of 10% over the prediction error of the Lasso (0.397). Note that in the training set the MSE satisfies the imposed constraint, as expected. Also note that the improvement in Group 1 is at the expense of slightly increasing the prediction error over Group 2. In terms of sparsity, the CSCLasso model has needed an additional predictor variable comparing to Lasso in order to comply with the constraint. As it will be seen in Section 2.2, the novel approach is set up by adding convex quadratic constraints to the Lasso formulation, and aims to control the performance measure on certain groups of interest. Other approaches have considered constrained versions of the Lasso before, see for example James et al. [2020], Gaines et al. [2018], Torres-Barrán et al. [2018], Hu et al. [2015] and references therein. In such works, equality or/and inequality linear constraints are considered for imposing prior knowledge and structure onto the coefficient estimates. In our approach instead, quadratic convex constraints are formulated and thus, our approach and results generalize those previously obtained in the literature. Not only constrained versions of the Lasso can be found in the literature. Indeed, many different variants have been proposed. In addition to those introduced in Chapter 1, e.g. the elastic net or the group lasso, other extension is to consider ˆ βLasso(λ) = arg min β 1 nky−Xβk2+λkAβk1.(2.2) instead of (2.1), where Ais a fixed matrix (see Tibshirani and Taylor [2011]). If A= (0|Ip), then the Lasso objective function is obtained; however, other forms of Adifferent from the identity can be found in the literature, see for example Ollier and Viallon [2017]. In fact, 20 Chapter 2. A cost-sensitive constrained Lasso various choices of Ain (2.2) define problems that are already well-known in the literature as the fused lasso [Tibshirani et al., 2005]. This chapter is structured as follows. In Section 2.2, the cost-sensitive constrained Lasso (CSCLasso) is introduced and some key issues are discussed. Section 2.3 considers theoretical properties of the CSCLasso, as the existence and uniqueness of solution, limit behaviour (in terms of the penalty parameter) and consistency. Section 2.4 presents a detailed numerical analysis with both simulated and real datasets, and finally, some conclusions are provided in Section 2.5. Technical proofs are relegated to Appendix A. 2.2 The cost-sensitive constrained Lasso: definition and key aspects This section presents the cost-sensitive constrained Lasso, which, as will be seen, is defined through an optimization problem with constraints related to prediction errors for individuals of interest. In addition, some computational details, as well as different key aspects concerning the tuning parameters of our proposal, are presented. 2.2.1 Definition The proposed CSCLasso is a novel variant of the Lasso where we shall demand that the prediction errors for the groups of interest are below certain threshold values, min β 1 n0ky0−X0βk2+λkAβk1 s.t. 1 n1ky1−X1βk2−f1≤0, . . . 1 nLkyL−XLβk2−fL≤0. (2.3) In the previous formulation, (y0,X0)is the set of observations used to build the sparse model with overall minimum MSE, which can be the complete dataset (y,X), or a subset of smaller size. Additionally, let (yl,Xl),l= 1, . . . , L, define groups of interest (not necessarily disjoint), where the MSE predictions are to be controlled. Then, nlis the number of instances related to group l. Finally, f= (f1, . . . , fL)contains the different threshold values for the MSE on the different groups. The solution of optimization problem (2.3) will be denoted by ˆ βCSCLasso(λ). From the formulation (2.3) it is natural to wonder whether running a Lasso on just the groups of interest is more advantageous. However, if a single Lasso is run on the groups of interest, dramatically bad predictions can be obtained when the resulting model is applied to new observations outside those groups, which is not the case for our approach. The same issue arises when a different Lasso model is built on each group of interest, but, in addition, new 2.2. The cost-sensitive constrained Lasso: definition and key aspects 21 observations are not given with their group of origin. Contrary to what happens with our novel approach (2.3), the Lpredictions obtained through the Ldifferent estimated Lasso models may not be suitable to give a final prediction for such new samples. The proposed method can be formulated as a Lasso with weighted quadratic penalties in the objective function associated with the different groups, but finding real meaning to their parameters (one per group) to be chosen is not an easy task (see Carrizosa and Romero Morales [2001] and the references therein) and the full control over the accuracy prediction on certain individuals of interest would disappear. However, the parameters f= (f1, . . . , fL)involved in our model have a clear interpretation and, in addition, this formulation enables us to bound the prediction errors in the different groups of interest. As an example, in the illustration of the method in Section 2.1 related to the prostate dataset, whereas the training set was used in the objective function with A= (0|I8), the prediction error over the young population of the training set, (y1,X1), is controlled through a performance constraint (f1= 0.357). In a real application, once the Lgroups of interest are selected by the user, threshold values f1, . . . , fLhave to be fixed. Note that these thresholds will depend directly on the dataset in question and the considered groups of interest. As a first option, they could be fixed by the user according to her demand, but therefore unfeasibility problems may appear when solving the CSCLasso problem (2.3). For that reason, in Section 2.2.3 two procedures for determining such threshold values so that (2.3) is feasible are given. Next, some other aspects related to the formulation of the CSCLasso and its resolution will be discussed. 2.2.2 Computational details The CSCLasso problem as defined by (2.3) is a non-differentiable convex optimization problem with quadratic and convex constraints. However, if we rewrite the non-differentiable term in (2.3) as Aβ=u+−u−, where u+= (u+ 1, . . . , u+ p)and u−= (u− 1, . . . , u− p)are new vectors of positive auxiliary variables, a differentiable version for the CSCLasso problem (2.3) is obtained in a straightforward manner as 22 Chapter 2. A cost-sensitive constrained Lasso min β,u+,u− 1 n0ky0−X0βk2+λ p X j=1 u+ j+λ p X j=1 u− j s.t. 1 n1ky1−X1βk2−f1≤0, . . . 1 nLkyL−XLβk2−fL≤0, Aβ=u+−u−, u+,u−≥0. This previous smooth formulation for the CSCLasso eases its resolution notably, since efficient solvers for quadratically constrained programming problems, such as Gurobi [Gurobi Optimization, 2018], are available. In particular, the Gurobi R interface will be used in this work to obtain all numerical results. Another remark concerning the formulation of the CSCLasso is that, instead of using the sum of squared deviations, least absolute deviations could have been considered. Then, (2.3) would be reduced to a regression problem under linear inequality constraints, as those described in James et al. [2020], Gaines et al. [2018] and Hu et al. [2015]. Nevertheless, to cope the nondifferentiability of the absolute value function, a huge number of constraints and new auxiliary variables, which would depend on n, should have been added. Consequently, these constrained approaches are likely to face severe numerical difficulties in practice for large datasets. 2.2.3 The choice of threshold values As commented in Section 2.2.1, threshold values f1, . . . , fLcould be fixed by the user. If the user is too demanding, imposing very low MSE threshold values for (some of) the different groups, the optimization problem may become unfeasible. Although a try-and-error procedure may be used, it would be very helpful to have strategies yielding feasible solutions. Here we propose two procedures for determining f1, . . . , fLin such a way that (2.3) is feasible. First, we propose a choice of the threshold values so that they are close to the OLS results, fl= (1 + τ)MSEl(ˆ βols), l = 1, . . . , L, (2.4) where MSEl(β) = 1 nlkyl− Xlβk2,l= 1, . . . , L and τ≥0is a small parameter whose meaning is the percentage of worsening with respect to the OLS prediction error. For the numerical example in Section 2.1, we could have imposed the threshold for the MSE over Group 1 equal to 0.391, which is a 10% (τ= 0.1) more than MSE1(ˆ βols)=0.355. The choice (2.4) deals with the heterogeneity coming from the variability of the different groups 2.2. The cost-sensitive constrained Lasso: definition and key aspects 23 (MSElis different across groups). Nevertheless, when heterogeneity related to the importance of each group is also considered, the parameter τcan be replaced in (2.4) by τl,l= 1, . . . , L. Next, we shall compute the minimum value of τ,τmin, so as to (2.3) is feasible. That is, the minimum τso that there exists β∗satisfying max l MSEl(β∗) MSEl(ˆ βols)!−1≤τ, and, therefore, τmin will be given as τmin = max l MSEl(β∗) MSEl(ˆ βols)!−1. Such τmin can be found as the optimal value of the following linear problem with convex and quadratic constraints min β,z z s.t. z≥MSEl(β) MSEl(ˆ βols)−1,∀l= 1, . . . , L. (2.5) The feasible version of the CSCLasso optimization problem can be formulated as min β 1 n0ky0−X0βk2+λkAβk1 s.t. 1 n1ky1−X1βk2−(1 + τ)MSE1(ˆ βols)≤0, . . . 1 nLkyL−XLβk2−(1 + τ)MSEL(ˆ βols)≤0, (2.6) where τ≥τmin. Finally, note that if τis big enough, then solving (2.6) is equivalent to solve the unconstrained problem. Indeed, it is possible to find the value of τ,τmax(λ), such that both the constrained and unconstrained problems are equivalent τmax(λ) = max l∈{1,...,L} MSEl(ˆ βLasso(λ)) MSEl(ˆ βols)−1.(2.7) A second possible choice for the threshold values follows an analogous approach but, instead of considering the results of the OLS, we shall consider the mean squared error of the 24 Chapter 2. A cost-sensitive constrained Lasso Lasso, as in the numerical example introduced in Section 2.1. For each l= 1, . . . , L, fl= (1 −γ)MSEl(ˆ βLasso(λ)), l = 1, . . . , L, (2.8) where γ≥0is related to the desired percentage of improvement over the Lasso solution (γ= 0.1in the numerical example of Section 2.1). In this case, we will compute the maximum value of γ,γmax, in such a way that (2.3) is feasible under (2.8), and the linear problem associated with γmax is max β,z z s.t. 1−MSEl(β) MSEl(ˆ βLasso(λ)) ≥z, ∀l= 1, . . . , L. (2.9) Thus, another possible feasible version of the CSCLasso optimization problem can be formulated as min β 1 n0ky0−X0βk2+λkAβk1 s.t. 1 n1ky1−X1βk2−(1 −γ)MSE1(ˆ βLasso(λ)) ≤0, . . . 1 nLkyL−XLβk2−(1 −γ)MSEL(ˆ βLasso(λ)) ≤0, (2.10) where γ≤γmax. Note that the two choices previously described for selecting the threshold values are not unique. Indeed, instead of using the MSE, another statistical measure as the R-squared can be considered. Further details about how they perform in numerical applications are described in Sections 2.2.4 and 2.4. 2.2.4 The role of the tuning parameters The CSCLasso, as defined by (2.6) or (2.10), is stated in terms of two tuning parameters, λ and τor λand γ, respectively. The first one, λ, is related to the sparsity of the solution, and the second one is linked to the user’s demanding level, since the degree of requirement increases as τ→τmin (or γ→γmax). In this section we investigate how the solution of the CSCLasso changes when λand τjointly vary (analogous results are obtained if λand γare analyzed instead). With this purpose, consider again the experimental setting as in the example of Section 2.1 related to prostate dataset with A= (0|Ip)(Lasso objective function), but in this occasion assume that the prediction errors of both groups (the young and the elderly people) shall be controlled. 2.2. The cost-sensitive constrained Lasso: definition and key aspects 25 The interval of variation of the parameter λis set to Iλ= [0,30]. Moreover, according to (2.5), the smallest value of τsuch that the CSCLasso optimization problem (2.6) is feasible is τmin = 0.055. On the other hand, following (2.7), τmax = max λ∈[0,30] τmax(λ)=2.355, although we will enlarge the interval of variation of τto also visualize the unconstrained solution; such interval will be finally set as Iτ= [τmin, τmax + 2] = [0.055,4.355]. Figure 2.1 represents, via a heat map, the solution for ˆ βCSCLasso 1(λ)for the different values of (λ, τ)in a grid contained in Iλ×Iτ. Some conclusions can be drawn from the figure. Consider first the cases where τand λare big enough. Since, in this case, τ≥τmax, then, as commented at the end of the previous section, solving (2.6) is equivalent to solving the Lasso. Therefore, ˆ βCSCLasso(λ) = ˆ βLasso(λ) = 0will be the optimal solution, provided that λis big enough. Analogously, if τ≥τmax but λ is small, then ˆ βCSCLasso(λ) = ˆ βLasso(λ), which will be equal to zero or not depending on the importance of the variable. When τis small, the constraints are demanding and, even for large values of λ, it might happen that ˆ βCSCLasso 1(λ)6=ˆ βLasso 1(λ)=0, as it is the case. Figure 2.1: Heat map of ˆ βCSCLasso 1(λ)using prostate dataset Figure B.1 (see Appendix B for further results) represents the analogous heat maps concerning ˆ βCSCLasso 2(λ),...,ˆ βCSCLasso 8(λ). A similar discussion as with ˆ βCSCLasso 1(λ)is applicable to these figures. An interesting remark to be made concerns the importance of each variable: while variable 1is the only one selected for the Lasso, the CSCLasso returns a less sparse solution in this case, since predictor variables 1, 2, 4 and 5 turn out to be significant. 26 Chapter 2. A cost-sensitive constrained Lasso However, this is not the rule, since there are examples where the level of sparsity is higher for the CSCLasso, as will be shown in Section 2.4. 2.3 Theoretical properties In this section we discuss some theoretical results concerning the CSCLasso model. In Section 2.3.1, the existence of a unique optimal solution to Problem (2.3) is proven for a fixed value of λ≥0. Section 2.3.2 deals with the limit behavior of the solution when λapproaches infinity. Finally, some consistency properties of the CSCLasso solution are derived in Section 2.3.3 from the Sample Average Approximation theory (see Shapiro et al. [2009]). 2.3.1 Existence and uniqueness of solution For constrained versions of the Lasso in the literature, as in James et al. [2020], it is not possible to obtain the path of solutions and therefore approximations are made with the use of numerical algorithms. For the CSCLasso problem, a closed form solution of expression ˆ βCSCLasso(λ)is not available. However, an implicit characterization of the CSCLasso solution (with only one constraint) can be found as the following result states. Proposition 1. Consider the CSCLasso problem with one constraint, min β 1 n0ky0−X0βk2+λkAβk1 s.t. 1 n1ky1−X1βk2−(1 + τ)MSE1(ˆ βols)≤0,(2.11) where A= (0|Ip)and assume that X0and X1are maximum rank matrices. Then ˆ βCSCLasso(λ) = 1 n0X0 0X0+1 n1 η(λ)X0 1X1−11 n0X0 0y0+1 n1 η(λ)X0 1y1−1 21 n0X0 0X0+1 n1 η(λ)X0 1X1−1 b(λ) (2.12) where η(λ)is the Lagrange multiplier of the constraint and the component s,s= 0,1, . . . , p, of the vector b(λ)is given by bs(λ) =      λ, if ˆ βCSCLasso s(λ)>0, −λ, if ˆ βCSCLasso s(λ)<0, 0else. From the previous proposition, it is clear that a closed form solution is hard to be obtained, even in the simplest scenario. 2.3. Theoretical properties 27 Nevertheless, given a fixed value of λ, the CSCLasso problem can be solved using quadratic programming via any of the standard solvers available in the literature. As an example, Figure 2.2 depicts the path of solutions for the prostate example introduced in Section 2.1, for an assortment of values of λin a grid (see Section 2.3.2 for details). Each line represents a component of ˆ βCSCLasso(λ),ˆ βCSCLasso j(λ)with j= 1,...,8. It can be observed from the figure that, contrary to what happens in the Lasso path of solutions (top panel of Figure 2.2), the CSCLasso path of solutions is not piecewise linear (bottom panel of Figure 2.2). Such non-linearity (due to the quadratic constraints) hinders the application of an iterative algorithm to obtain the path of solutions as those given in papers James et al. [2020] and Gaines et al. [2018]. Also note from Figure 2.3 that as a consequence of the performance constraints, the solution is stabilized when λincreases, but does not shrink to 0, as with Lasso. This is detailed in Section 2.3.2. Even without having the expression of the general solution of (2.3), we next prove that, under full rank assumptions, the solution is unique. First, in order to simplify the formulation of (2.3), henceforth its feasible set will be denoted by B, which is convex and closed. This is also true (and the results which follow remain valid) if, on top of the performance constraints, one adds linear constraints modeling other aspects of interest (for example, the sign of a certain coefficient can be fixed to be positive or negative depending on the known relation between the corresponding predictors with the response variable). In the same vein, henceforth, (y0,X0)=(y,X)is used to minimize the overall MSE. In this way, the CSCLasso problem (2.3) is rewritten as min β∈B 1 nky−Xβk2+λkAβk1.(2.13) The following result guarantees that the solution of problem (2.13) is unique. Theorem 1. Consider Problem (2.13) where Xis assumed to be a maximum rank matrix and its feasible region Bis a convex and closed set in Rp+1. Then, Problem (2.13) has a unique optimal solution. 2.3.2 Asymptotic behaviour One of the key points when dealing with Lasso-type problems is the choice of the regularization parameter λ. In the case of the Lasso, such choice is straightforward since the entire path of solutions is known to be piecewise linear, shrinking to 0. In particular, it is known that there exists a value of λ,λ∗, such that the solution ˆ βLasso(λ) = 0is optimal for all values λ≥λ∗. The following result provides explicitly the value of such λ∗. Proposition 2. Consider the Lasso model (2.2). Define λ∗as the optimal value of the linear 28 Chapter 2. A cost-sensitive constrained Lasso Figure 2.2: Path of solutions under Lasso (top) and CSCLasso (bottom) for prostate dataset 2.4. Numerical experiments 35 Figure 2.4: Median overall MSE over the testing sets and NZ percentage under the choice p= 100 Next, we test how the solution of the CSCLasso behaves with respect to other global performance measures as the l2distance, false positive and negative rates, which are defined not in terms of prediction errors, but on the correct fitting of the generator process (see Yu and Liu [2016]). In particular, the l2distance is defined as kˆ β−βk2, where βis the vector of coefficients that generated the datasets (described at the beginning of this section), and ˆ β(λ) are the estimators. In addition, the false positive rate (FPR) and false negative rate (FNR) are calculated as follows: FPR =|j:βj= 0 & ˆ βj(λ) = 0 | |j:βj= 0 |, FNR =|j:βj6= 0 & ˆ βj(λ)=0| |j:βj6= 0 |, where j= 1, . . . , p. The median of these three measures as well as the median of the overall MSE (already shown in Figures 2.4-2.5 and Figures B.2, B.3 and B.4 in Appendix B), are presented in Table 2.3. For the choices where p= 20, the FPR values are not given since all the predictors have associated non-zero coefficients when the datasets were created. From this table it can be deduced that similar or even better results, comparing with those of the Lasso, are obtained in the majority of the cases across the four different measures. A final remark concerning the computational cost of the CSCLasso when comparing with Lasso is as follows. The median user time required to solve the problem with the largest dataset considered in this study (nk= 300 and p= 500) is 0.85 seconds when the Lasso is run on Intel(R) Core(TM) i7-7500U CPU at 2.70GHz 2.90GHz with 8.0 GB of RAM; whereas the CSCLasso requires 6.60 seconds. Nevertheless, to better understand how the computation time behaves depending on pvalue, a grid in this parameter has been inspected, while nkis set to 36 Chapter 2. A cost-sensitive constrained Lasso Figure 2.5: Median MSEkover the testing sets for k= 7,...,20 under p= 100 features and, nk= 150 (top) and nk= 300 (bottom). Each subgraph represents one group and the Y-axis shows the different percentages of improvement 300. Figure 2.6 depicts the logarithm of the user times in seconds obtained under Lasso and CSCLasso models when nk= 300 and pchanges. Then, under a reasonable computational cost, the desired results are achieved. Further analyses regarding the computational times are shown in the Appendix B (Figure B.5). 2.4. Numerical experiments 37 nk= 150 nk= 300 pOverall MSE l2distance FPR FNR Overall MSE l2distance FPR FNR 20 Lasso 0.667 4.149 - 0.050 0.666 4.146 - 0.000 Improv. 3% 0.691 4.108 - 0.000 0.673 4.099 - 0.000 Improv. 5% 0.695 4.104 - 0.000 0.683 4.089 - 0.000 Improv. 7% - - - - 0.693 4.077 - 0.000 Improv. 10% - - - - - - - - Improv. 15% - - - - - - - - Improv. 20% - - - - - - - - 100 Lasso 0.732 2.931 0.275 0.250 0.734 2.729 0.250 0.150 Improv. 3% 0.740 2.857 0.388 0.100 0.735 2.721 0.163 0.150 Improv. 5% 0.741 2.855 0.375 0.100 0.734 2.721 0.163 0.100 Improv. 7% 0.742 2.853 0.400 0.100 0.734 2.718 0.163 0.100 Improv. 10% 0.745 2.849 0.400 0.100 0.732 2.704 0.275 0.100 Improv. 15% 0.751 2.844 0.425 0.100 0.735 2.678 0.288 0.050 Improv. 20% 0.759 2.835 0.463 0.100 0.738 2.655 0.288 0.100 500 Lasso 0.772 2.778 0.135 0.300 0.744 2.750 0.065 0.300 Improv. 3% 0.772 2.779 0.040 0.350 0.744 2.658 0.063 0.250 Improv. 5% 0.772 2.775 0.050 0.300 0.745 2.653 0.075 0.250 Improv. 7% 0.772 2.769 0.042 0.350 0.746 2.647 0.077 0.250 Improv. 10% 0.772 2.760 0.056 0.350 0.748 2.632 0.104 0.250 Improv. 15% 0.771 2.754 0.063 0.300 0.752 2.640 0.129 0.250 Improv. 20% 0.774 2.735 0.069 0.350 0.760 2.607 0.148 0.200 Table 2.3: Median performance measures over testing sets for synthetic datasets Figure 2.6: A two-dimensional graph of the logarithm of the user times in seconds for nk= 300 as pincreases 38 Chapter 2. A cost-sensitive constrained Lasso 2.4.2 Leukemia dataset: a gene expression dataset The real stratified dataset described in Kouno et al. [2013] is explored here. The data contain information related to myeloid monocytic leukemia cells undergoing differentiation to macrophages. In particular, the dataset is formed by expression levels of 45 transcription factors (response and predictor variables) measured at 8distinct times (groups) of the differentiation process. As in Ollier and Viallon [2017], the aim is to predict the EGR2 transcription factor in terms of the other p= 44 factors. The sample size per group is equal to 120. Similarly as in Section 2.4.1, the Lasso was run and the overall prediction errors, individual prediction errors per group and percentage of non-zero coefficients are recorded. The records in Group 1 yield the best MSE using Lasso model. Therefore, we may be interested in obtaining an even better fitting for such data. The CSCLasso problem is solved with threshold values smaller than the Lasso error, which turns out to be 0.370. Table 2.4 shows the obtained median results for an assortment of improvement levels, namely, 5%,7%,10% and 15% or, equivalently, γ is equal to 0.05,0.05,0.07 and 0.15, respectively. From that table, it can be seen how the prediction error of interest (corresponding to Group 1) decreases with the improvement level, as expected. Similarly as in Section 2.4.1, the overall mean squared error does not exhibit significant changes, while the prediction errors in the uncontrolled groups do not exhibit the same behaviour. Some of them slightly improve (Group 6), others slightly worsen (as the Group 5 and Group 8) and others remain constant (as Group 2). Finally, in regards to the sparsity of the solution, for this dataset, less sparse solutions are obtained by CSCLasso in comparison with the Lasso ones. f1Overall MSE MSE1MSE2MSE3MSE4MSE5MSE6MSE7MSE8NZ Lasso - 0.620 0.370 0.417 0.902 0.496 0.480 0.535 0.496 0.685 53.33 Improv. 5% 0.352 0.623 0.357 0.417 0.918 0.500 0.555 0.512 0.497 0.719 73.33 Improv. 7% 0.344 0.626 0.348 0.418 0.915 0.504 0.565 0.514 0.497 0.722 66.67 Improv. 10% 0.333 0.630 0.335 0.418 0.911 0.510 0.574 0.513 0.498 0.728 66.67 Improv. 15% 0.315 0.636 0.331 0.415 0.903 0.523 0.591 0.512 0.502 0.737 73.33 Table 2.4: Median errors over testing set for gene expression dataset. Constraints imposed over Group 1 2.4.3 Communities and Crime dataset In this section, a real dataset from the UCI Machine Learning Repository [Lichman, 2013] will be analyzed. In particular, the so-called Communities and Crime Unnormalized Data Set shall be considered. The dataset is about communities within the United States and has already been inspected in the literature (see Redmond and Baveja [2002]). This dataset combines crime information from the FBI databases [U.S. Department, 1995] as well as socio-economic and law enforcement data from U.S. Department [1992a] and U.S. Department [1992b], respectively. The dataset is formed by p= 124 predictors, 23 of which present missing values, 2.5. Chapter summary 39 and n= 2215 instances, where the response variable measures the number of murders per 100K population. The predictor variables with missing values are not consider for the next experiments. As such, we finally consider p= 101 predictors. Additionally, for each instance (community), the region from which it comes is known. Thus, if we were interested in obtaining a good prediction in a certain region, say Midwest, we could control these communities by including a performance constraint. Table 2.5 shows the median errors over the testing set for Group 1, formed by the communities of Midwest, and over the rest of communities (Group 2). In terms of overall MSE and MSE over the two groups, similar conclusions as in Section 2.4.2 are drawn. Whereas different improvement levels are imposed, the MSE of interest (MSE1) is getting smaller but the overall prediction error is almost not affected by the constraint. Lastly, regarding the sparsity of the solution, an analogous behaviour as that observed in the case of simulated data is obtained: the solution becomes less sparse with the improvement level. f1Overall MSE MSE1MSE2NZ Lasso - 0.488 0.433 0.453 21.57 Improv. 5% 0.411 0.488 0.422 0.453 25.49 Improv. 7% 0.403 0.487 0.420 0.453 28.43 Improv. 10% 0.390 0.488 0.416 0.453 26.47 Improv. 15% 0.368 0.486 0.403 0.459 34.31 Table 2.5: Median errors over testing set for communities and crime dataset. Constraints imposed over Group 1 As previously commented, the groups of interest may overlap. As an illustration, assume that the interest is in controlling the prediction error in communities of Midwest or communities with a population density larger than or equal to the 75th percentile. Let Group 1 denote the communities from Midwest, while Group 2 represents the communities where the density of population is higher than the 75th percentile. For instance, if we aim to improve in a 7% the errors obtained by the Lasso model (equal to 0.513 and 0.442), then the CSCLasso results become 0.475 and 0.441, respectively. 2.5 Chapter summary In this work a new version of the Lasso regression model that strives to control the performance rates associated with individuals of interest is proposed. The method has a significant application in the context of heterogeneous data, where it is common that certain sources are more reliable than others, or simply the prediction on some groups of data are of higher interest, and thus a better fit is sought for some data. In order to control the individuals of interest, performance constraints are included in the regression model. This approach leads to a novel method (CSCLasso) which is not reported in the literature previously, up to our knowledge. Theoretical results concerning this novel methodology have been discussed and, in addition, the CSCLasso has been tested on six synthetic datasets with different properties, on a well-referenced real 40 Chapter 2. A cost-sensitive constrained Lasso stratified biomedical dataset and on a real social sciences dataset. The numerical section shows that, with a low computational cost, the accuracy prediction errors for the groups of interest are controlled. This is done at the expense of reducing sparsity (if the regularization parameter is kept fixed) or the overall accuracy. Chapter 3 On linear regression models with hierarchical categorical variables 42 44 Chapter 3. On linear regression models with hierarchical categorical variables In this chapter we study the mathematical optimization problem that trades off, in linear regression models, accuracy and model complexity, in the presence of categorical variables that have a hierarchical structure, with their categories arranged as a directed tree. In the literature, this kind of data appears in different fields of research, such as nested spatial data in Spatial Statistics [Gotway and Young, 2002], behavioral data in Retail Business Analytics [Griva et al., 2018], and economic activity data in Official Statistics [European Commission, 2008; Katz-Gerro and López Sintas, 2019]. 3.1 Introduction Let J0be the set of continuous and dummy predictor variables, whereas Jthe set of hierarchical categorical predictor variables. Then, consider the random vector (Y, X0,X), where X0 denotes the vector of the predictor variables in J0,Xdenotes the vector of categorical predictor variables in J, and Ydenotes the response variable. In the real-world dataset cancer-reg [Rippner, 2017] used in the numerical section, with individuals from the United States of America (U.S.), geography is a categorical variable with a hierarchical structure. According to the U.S. Department of Commerce Economics and Statistics Administration and the U.S. Census Bureau,geography can be coded using the states (51 in total), which is the highest level of granularity for which information is available in the dataset. This means that 51 coefficients need to be estimated for this variable, where individuals in the same state share the same coefficient in the linear regression model. The variable geography can alternatively be coded using the subregions, such as East-South Central,Middle Atlantic and New England, where each state belongs to exactly one of the 9 subregions. Consolidating individuals at the subregions, sharing the same coefficient, yields a lower level of granularity for geography, where, instead of 51, only 9coefficients need to be estimated and interpreted. The individuals can be further consolidated into 4 regions, namely West,South,Mid-West and North-East, where only 4coef- ficients would be associated to geography in the reduced linear regression model. Using these regions, one has the least granular representation of geography. This work is devoted to trading off accuracy of the linear regression model and its complexity, measured as a cost function of the level of granularity used to represent each of the hierarchical categorical variables. The categories of hierarchical categorical variable j∈ J can be arranged as a directed tree Tj, i.e., a directed graph with a root node, r(Tj), and a unique path from each node to r(Tj). In addition, let V(Tj)denote the set of nodes in the tree and L(Tj)⊂ V(Tj)the set of leaf nodes. See Figure 3.1 for the tree associated with the categories of geography, where the leaf nodes correspond to the states, going upstream we find the subregions and then the regions, which, in turn, are directly connected with the root node. Let (yi,x0 i,xi)be the vector associated with individual i, with x0 i= (x0 ij0)and xi= (xijv), where xijv is equal to 1 if individual ibelongs to category v∈ V(Tj)of variable j∈ J. If we were to use the most granular representation of the hierarchical categorical variables, we would need to use the categories associated with the 3.3. Numerical experiments 51 U.S. West Pacific WA ... OR Mountain ID ... MT South West-South Central OK ... LA East-South Central TN ... KY South Atlantic NC ... VA Mid-West West-North Central KS ... IA East-North Central IL ... OH North-East Middle Atlantic New England (a) S∗ 1when MSE∗(S∗ 1) = 0.408 and c= 44 U.S. West Pacific WA ... OR Mountain South West-South Central OK ... LA East-South Central South Atlantic Mid-West West-North Central KS ... IA East-North Central North-East (b) S∗ 1when MSE∗(S∗ 1) = 0.421 and c= 21 U.S. West Pacific Mountain South West-South Central East-South Central South Atlantic Mid-West West-North Central KS ... IA East-North Central North-East (c) S∗ 1when MSE∗(S∗ 1) = 0.427 and c= 14 Figure 3.3: Less granular representations for the geography variable in the cancer-reg dataset 52 Chapter 3. On linear regression models with hierarchical categorical variables of CRIM into two groups: those whose values are below (node M1,1) and above (node M1,2) the median. Second, the quartiles are used to subdivide M1,1(nodes Q1,1and Q1,2) and M1,2 (nodes Q1,3and Q1,4) into two nodes. This way we examine the thresholds of the continuous predictor variables required to predict the response variable. Figure 3.4: Tree associated with the variable CRIM in the housing dataset after being discretized CRIM M1,1 Q1,1Q1,2 M1,2 Q1,3Q1,4 When solving Problem (3.11) for the 37 values of cin the set {12,...,48}, we obtain the pareto frontier in Figure 3.5. The MSE of the model with the highest granularity for all hierarchical variables is 23.06. When we start reducing the granularity the MSE remains approximately the same. Actually, when cis reduced from 48 to 30, the accuracy is barely damaged but the complexity of the linear regression model is dramatically improved. Figures 3.6-3.7 show the subtrees S∗ jfor all j∈ J for the solution in Figure 3.5 that achieves the minimum AIC. In this solution, we can observe how variables INDUS,AGE and RAD are eliminated from the linear regression model, as their root node is the only one selected with a coefficient equal to zero. By contrast, we require the highest level of granularity for PTRATIO,Band LSTAT. For DIS, the linear regression model only needs to know whether the predictor variable is below the median. For the remaining predictor variables, leaf as well as non-leaf nodes are selected. 3.3.3 The synthetic data In this section we illustrate our approach on synthetic data. The data generating model is yi=X j∈J X l∈L(Tj) βS jlxijl +εi, i = 1, . . . , n, (3.15) where |J| = 2 and J0=∅. The values of the coefficients βS jl,l∈ L(Tj), are given in Figure 3.8. Note that the first two leaf nodes of T1have the same coefficient, and the same holds for the other two leaf nodes. Therefore, the tree can be pruned to avoid unnecessary splits, yielding 3.3. Numerical experiments 53 Figure 3.5: Pareto frontier for MSE versus the number of coefficients to be estimated in the reduced model for the hierarchical categorical variables in the housing dataset Figure 3.6: Less granular representations for the first six hierarchical categorical variables in the housing dataset for the solution in Figure 3.5 with MSE∗((S∗ j)j∈J ) = 23.37 and c= 32. Note that this is the solution that achieves the minimum AIC CRIM M1,1 Q1,1Q1,2 M1,2 ZN M2,1M2,2 Q2,3Q2,4 INDUS NOX M4,1M4,2 Q4,3Q4,4 RM M5,1M5,2 Q5,3Q5,4 AGE 54 Chapter 3. On linear regression models with hierarchical categorical variables Figure 3.7: Less granular representations for the last six hierarchical categorical variables in the housing dataset for the solution in Figure 3.5 with MSE∗((S∗ j)j∈J ) = 23.37 and c= 32. Note that this is the solution that achieves the minimum AIC DIS M7,1M7,2RAD T AX M9,1 Q9,1Q9,2 M9,2 P T RAT IO M10,1 Q10,1Q10,2 M10,2 Q10,3Q10,4 B M11,1 Q11,1Q11,2 M11,2 Q11,3Q11,4 LST AT M12,1 Q12,1Q12,2 M12,2 Q12,3Q12,4 the subtree in Figure 3.9. The same holds for T2. The error is taken εi∼N(0, σ2)for different values of σ2given below. We have n= 3000 individuals, evenly distributed across the different combinations of categories l1∈ L(T1)and l2∈ L(T2). The purpose of this section is twofold. First, we illustrate how our approach is able to recover the pruned tree underlying our synthetic data. Second, we carry out a study for assessment of the generalization error of the final chosen model. Let us consider σ2= 0.04 and solve Problem (3.11) for the 10 values of cin the set {2,...,11}. Figure 3.10(a) shows the pareto frontier for the number of coefficients to be estimated in the reduced model versus the MSE. For small values of MSE, the chosen nodes are the 8 green leaf nodes in Figure 3.9, which implies that our methodology is able to successfully detect the pruned tree underlying each hierarchical categorical variable in our data. Similar conclusions can be drawn when σ2= 0.16 (Figure 3.10(b)) and σ2= 0.36 (Figure 3.10(c)). To end the numerical section, we provide an estimation for the MSE and the complexity of the reduced model using a 10-fold cross validation approach, showing that our procedure works properly with the available individuals (training sample), but also for future individuals (testing sample). For each fold, the training set is used to solve Problem (3.11) and get S∗ j,j∈ J. Once the subtrees are found, and thus the reduced linear regression model, we calculate its MSE for the training and testing sets, which are plotted in Figures 3.11(a)-3.11(c) for the different values of σ2. As can be observed, the MSE values for the training sets (red lines) are only slightly 3.4. Chapter summary 55 Figure 3.8: Trees associated with the two hierarchical categorical variables in the synthetic dataset together with βS jl,l∈ L(Tj) 0.2 0.2−0.3−0.3 (a) Tree T1for hierarchical categorical variable j= 1 0.4 0.4−0.20.1−0.10.2 −0.8 (b) Tree T2for hierarchical categorical variable j= 2 smaller than those for the testing sets (blue lines). Then, in view of results, we can conclude that our methodology generalizes well. 3.4 Chapter summary In this work a new methodology to deal with hierarchical categorical variables, i.e., categorical variables that can be measured at different levels of granularity, has been developed. Through a Mixed Integer Convex Quadratic Problem with Linear Constraints, we study the tradeoff between accuracy and model complexity. Our proposal has been tested on both real-world and synthetic datasets. The numerical section shows that much less granular representations for the hierarchical categorical variables can be found at the expense of slightly damaging the accuracy. 56 Chapter 3. On linear regression models with hierarchical categorical variables Figure 3.9: Pruned tree and less granular representation of the two hierarchical categorical variables in Figure 3.8 from the synthetic dataset 0.2−0.3 (a) Pruned tree for the first hierarchical categorical variable 0.4 −0.20.1−0.10.2 −0.8 (b) Pruned tree for the second hierarchical categorical variable 3.4. Chapter summary 57 Tree Figure 3.1 Optimal Tree Figure 3.3(a) Optimal Tree Figure 3.3(b) Optimal Tree Figure 3.3(c) U.S. West Pacific WA 0.058 0.065 0.086 0.109 AK 0.720 0.716 0.651 CA -0.149 -0.138 -0.114 HI -0.978 -0.992 -0.990 OR -0.015 -0.010 0.002 Mountain ID -0.302 -0.313 -0.199 -0.161 WY 0.143 0.143 NV 0.481 0.499 UT -0.459 -0.466 CO -0.300 -0.303 AZ -0.473 -0.459 NM -0.273 -0.263 MT -0.131 -0.144 South West-South Central OK 0.692 0.683 0.705 0.567 AR 0.686 0.678 0.722 TX 0.416 0.408 0.432 LA 0.355 0.345 0.378 East-South Central TN 0.529 0.525 0.575 0.631 MS 0.545 0.539 AL 0.330 0.330 KY 0.681 0.683 South Atlantic NC 0.099 0.098 0.306 0.368 DE 0.043 0.056 FL 0.284 0.282 GA 0.129 0.121 MD 0.325 0.337 SC 0.373 0.369 WV 0.364 0.359 DC 0.350 0.393 VA 0.535 0.538 Mid-West West-North Central KS 0.638 0.659 0.538 0.580 MN 0.302 0.327 0.201 0.230 MO 0.581 0.575 0.574 0.611 NE 0.069 0.073 0.069 0.100 ND 0.123 0.132 0.097 0.103 SD 0.019 0.019 0.016 0.014 IA -0.077 -0.072 -0.088 -0.057 East-North Central IL 0.201 0.211 0.291 0.336 IN 0.527 0.526 MI 0.239 0.243 WI 0.145 0.148 OH 0.386 0.389 North-East Middle Atlantic PA -0.099 -0.102 -0.076 -0.025 NJ 0.074 NY -0.183 New England ME 0.327 0.121 VT 0.241 MA -0.053 RI 0.102 CT -0.311 NH 0.183 Table 3.1: Coefficients associated with four different representations for geography variable in the cancer-reg dataset 58 Chapter 3. On linear regression models with hierarchical categorical variables Variable Name Description Type Discretized Predictor CRIM Crime rate by town Continuous Yes ZN Proportion of residential land zoned for lots greater than 25.000 square feet Continuous Yes INDUS Proportion of nonretail business acres per town Continuous Yes NOX Nitrogen oxide concentrations Continuous Yes RM Average number of rooms Continuous Yes AGE Proportion of owner units built prior to 1940 Continuous Yes DIS Weighted distances to five employment centres Continuous Yes RAD Index of accessibility to radial highways Continuous Yes TAX Full value property tax rate ($/$10.000) Continuous Yes PTRATIO Pupil-teacher ratio by town school district Continuous Yes B Black proportion of population Continuous Yes LSTAT Proportion of population that is lower status Continuous Yes CHAS 1 if tract bounds river; 0 otherwise Binary No Response MEDV Median value of owner-occupied homes (in $1000’s) Continuous No Table 3.2: The predictor and the response variables in the housing dataset 3.4. Chapter summary 59 Figure 3.10: Pareto frontier for MSE versus the number of coefficients to be estimated in the reduced model for the synthetic dataset for different σ2values (a) σ2= 0.04 (b) σ2= 0.16 (c) σ2= 0.36 60 Chapter 3. On linear regression models with hierarchical categorical variables Figure 3.11: Average MSE (10-fold CV) versus the imposed threshold cwhen σ2changes (a) σ2= 0.04 (b) σ2= 0.16 (c) σ2= 0.36 4.3. A sparse Naïve Bayes 67 the variables may make the independence assumption less violated, and the accuracy improved (on top of the fact that a model with less variables is more explainable). In order to illustrate how the violation of the independence assumption may affect the performance of the NB, consider the next numerical example. A sample of size 2000 of a random vector (X1, X2, X3, X4)is simulated for two classes from a multivariate Normal distribution in such a way that the random variables are independent conditioned to the classes except for X1and X2which are correlated according to a Pearson coefficient of 0.95. A Gaussian NB classifier (that is, Xi|Ck∼Nµi,k, σ2 i,kfor i= 1,2,3,4,k= 1,2, where µi,k were randomly selected in the interval [1,7]) was run using all possible subsets of features and the results are shown in Table 4.1. The accuracy when all the variables are used is equal to 78.28, a value that is improved if the set {X1, X3, X4}is considered (accuracy equal to 79.94). Table 4.1: Performance rate for all possible combinations of features in a multivariate normal simulated example. Combination X1 X2 X3 X4 X1, X2 X1, X3 X1, X4 X2, X3 X2, X4 X3, X4 X1, X2, X3 X1, X2, X4 X1, X3, X4 X2, X3, X4 X1, X2, X3, X4 of variables ACC 68.08 68.51 68.55 69.01 68.38 74.99 75.08 74.86 75.25 75.68 73.58 73.85 79.94 79.84 78.28 Having illustrated that using just a subset of the features may improve accuracy, we face the combinatorial problem of finding the adequate set of features to be used. The previous brute force procedure, where all possible combinations of features are examined, turns out infeasible in practice, especially for large databases. Instead, in this work we propose a variable reduction method in which only certain combinations of features are sampled and evaluated. Such combinations, as will be seen in Section 4.3, shall be chosen by considering the dependencies among the features. 4.3 A sparse Naïve Bayes As commented in Section 4.2.2, considering all possible combinations of features to determine the best one is hard from a computational point of view, especially for large datasets since a total of 2p−1sets should be evaluated. The aim of this section is to describe an efficient methodology to guide the search of the subset of features, by inspecting only some subsets selected in terms of the dependence among features. As a result, a sparse, computationally tractable NB is obtained. 68 Chapter 4. Variable selection for Naïve Bayes classification 4.3.1 Description of the method The variable reduction strategy proposed in this section is based on a clustering of features made in terms of their dependencies. As commented in Section 4.2.1, the key assumption of the Naïve Bayes is the independence of the features conditioned to the class. The novel method presented in this work aims to preserve the independence assumption without damaging the predictive power of the classic NB. In other words, our methodology helps to select variables that are as independent as possible while provides good classification accuracy. To do that, we consider a dependence measure between random variables Xand Y, which increases with the degree of dependence between the variables. First, consider for i, j ∈ {1, . . . , p}and k= 1, . . . , K, the dependence between feature Xiand feature Xjconditioned on class Ck. In order to have a unique, summarized measure of dependence between Xiand Xj, let M be the matrix whose elements (Mij) represent the maximum dependence among all classes, between Xiand Xj. Note that such a choice represents the worst case scenario. A number of dependence measures proposed in the literature can be selected: Pearson correlation coefficient, Spearman’s rank-order correlation coefficient, Hoeffding Dstatistic (see Hoeffding [1948]), the mutual information coefficient (MI) [Linfoot, 1957], the Maximal Information coefficient (MIC) [Reshef et al., 2011] or the distance correlation coefficient [Székely et al., 2007], among others. We tried using these different measures and similar results were obtained (see Section 4.4 and the Supplementary Material in Appendix C). Therefore, since the mutual information measure enables us to work with both continuous and categorical variables and has been widely used in the literature [Kinney et al., 2010; Sharpee et al., 2004], we will select this measure. This coefficient quantifies the information about one variable Xprovided by a different variable Y, and it is defined as I(X, Y ) = ZYZX p(x, y) log p(x, y) p(x)p(y)dxdy, in the case of continuous variables. In the categorical case, the previous formula can be rewritten in terms of sums. The previous dependence measure can be computed by the function mutinformation from the infotheo package of the Statistical software environment R [R Core Team, 2017]. An illustration of the matrix Mfor the real dataset Statlog (Australian Credit Approval) (australian) from the UCI Machine Learning Repository [Lichman, 2013] is represented by Figure 4.1. The dataset concerns credit card applications, and is formed by 14 variables and two classes (+/-). Moreover, to visualize the different correlation patterns of the real-life datasets used throughout this work, at the end of Appendix C, the associated matrices Musing the MI measure are represented via heatmaps. Next, with the aim of performing a cluster analysis in terms of the degree of association among the features, a dissimilarity matrix Hof dimension p×pis defined in terms of matrix 4.3. A sparse Naïve Bayes 69 Figure 4.1: Heatmap associated to matrix M(based on MI correlation) corresponding to the australian dataset Mby the elements Hij =M∗−Mij M∗,(4.5) where M∗= max i,j∈{1,...,p}Mij. Note that, under the previous definition, the elements of Hare bounded below by zero, where this value represents the maximum degree of dependence. Moreover, the upper-bound of the elements of His one, which represents the minimum degree of dependence. Therefore, the higher the values of Hij are, the less dependence exists between Xiand Xj, according to the selected dependence measure. Note that, as described in the first section of Appendix C, the results obtained are rather robust regarding the dependence measure. Once the dependence measure is set, the classifier’s performance measure to be maximized in the embedded Variable Selection strategy can be chosen, among the previously described measures in Section 4.2.1, according to the user’s convenience and the properties of the dataset. Generally, the ACC is selected, but in some cases, e.g. for unbalanced datasets or when there exist critical classes, our proposal replace ACC with AUC, precision or a certain Recall. Thus, the novel method turns out specially advisable for datasets where the classes are unbalanced and/or of different importance. The selection of the dependence and the classifier’s performance measures is the first step of our algorithm (see Algorithm 4). Once we have chosen a dependence measure, and the elements of the matrix Hare computed, we perform a hierarchical cluster analysis of features according to the dissimilarity matrix H(step 2 of the algorithm). 70 Chapter 4. Variable selection for Naïve Bayes classification In the obtained dendrogram, the vertical axis represents the degree of dissimilarity. The higher the value of the height is, the less dependent the variables are, according to the dependence measure. For example, the dendrogram corresponding to the australian dataset is given by Figure 4.2. Such a dendrogram has been obtained using the routine hclust of the Statistical software environment R. In this case, V5and V6are highly dependent, as well as V9 and V10. However, the rest of variables are almost independent, since they cluster at heights between 0.9and 1. Figure 4.2: Cluster dendrogram (based on MI correlation) corresponding to the australian dataset Once the dendrogram is built, a (not necessarily regular) grid of a specified number Cof cuts along the height is fixed. The basic idea underlying the variable reduction strategy is to examine at each cut (or threshold) of the grid several combinations of features, in such a way that only one feature is selected per cluster since all elements in a cluster are assumed to be strongly dependent. As an example, consider Figure 4.2 and assume that one of the Ccuts is c= 0.74 (horizontal line). Then, we consider that there are 12 clusters: 10 clusters formed by only one feature and the clusters {V5, V6}and {V9, V10}. And, therefore, four independent combinations would be selected at this threshold: (V1, V7, V11, V3, V13, V2,V5, V4, V14, V8, V12,V9), (V1, V7, V11, V3, V13, V2,V5, V4, V14, V8, V12,V10),(V1, V7, V11, V3, V13, V2,V6, V4, V14, V8, V12,V9)and (V1, V7, V11, V3, V13, V2,V6, V4, V14, V8, V12,V10). Note that the higher (lower) the value of the cut, the more likely we are to choose independent (dependent) variables. Although the previous strategy reduces the computational cost of the brute force approach, it still may be costly for large datasets that originate a complex dendrogram with many combinations per threshold. In addition, removing some of the features from the combinations may 4.4. Numerical Illustrations 71 lead up to sparser and more accurate models, since it might happen that the (independent) variables selected in the combinations have a very low predictive power. In order to strive to avoid such inconveniences, we propose a refinement of the strategy as follows. First, a maximum number Sof combinations per threshold is set (if the total number of possibilities for a given threshold c,nc(c), does not exceed S, then all of them will be considered). In the previous example, nc(0.74) = 4. We should point out that parameters Cand Sare used to alleviate the computational burden, since, as commented before, Cfixes the number of cuts along the height in the dendrogram, and Sthe maximum number of combinations per threshold to be evaluated. Therefore, the higher Cand Sare, the higher the computer time is. For this reason, we will fix reasonable values for these parameters in Section 4.4.3. Second, for each cluster of variables to be examined, a value qrepresenting the probability of selecting this cluster for extracting randomly a variable to be included in the combination is also set. If we fix q= 0.4, the previous four combinations become (V1, V7, V8, V13, V14),(V1, V4, V6, V8, V9, V12, V14),(V2, V8) and (V4, V6, V10, V11), respectively. The parameter qis directly related to the sparsity degree: the lower the value of qis, the less variables are inspected (the expected number variables to be considered is equal to q×p). The choice of the values {C, S, q}will be discussed in Section 4.4. Once the set of combinations of features to be evaluated is reduced, the NB would be implemented and, its performance and feasibility on the constrains considered, evaluated for each combination. This is summarized in step 3 of the Algorithm 4. Finally, the feasible combination yielding the highest performance measure (accuracy, AUC or whatever chosen measure) would be considered the best, taking also into account the whole set of variables in this comparison (step 4). For australian, if no constraint is imposed, the features selected by our model are (V1, V4, V6, V8, V9, V12, V14), which achieve an ACC of 86.76, whereas the whole set of variables returns 85.29. According to the results, it can be deduced that our model has kept the important features, using only a half of the total set. However, in this dataset, the positive class (the load is granted) is the most risky. Then, if we impose e.g. that Recall + >92, the combination of features (V2, V8)would be the selected one. A summary concerning the strategy for the sparse NB is given by Algorithm 4. 4.4 Numerical Illustrations In this section, the behaviour and performance of our approach is illustrated throughout an extensive empirical study, using both simulated and real datasets. In the first case, a synthetic data set is simulated in order to test how the performance and level of sparsity of the proposed sparse NB changes with the level of dependence among the features. Second, ten real datasets from the UCI Machine Learning Repository [Lichman, 2013], presenting different correlation patterns, different degrees of unbalancedness and some of them combining both continuous 72 Chapter 4. Variable selection for Naïve Bayes classification Algorithm 4: Pseudo-code of the sparse NB 1. Select the dependence and the classifier’s performance measures. 2. Perform cluster analysis and build the dendrogram. 3. Variable reduction strategy: set specific values for the parameters {C, S, q}and initialize F=∅. for c= 1, . . . , C do for s= 1,...,min{nc(c), S}do (a) Obtain the s−th combination of features. For each cluster only one variable is randomly selected with probability q, and none with probability 1−q. (b) Construct the classifier for the s−th combination of features. (c) Evaluate the selected classifier’s performance measure and if feasible, add it to F. end end 4. Variable Selection: select the combination of variables leading to the best performance, among those in F. and categorical variables, will be analyzed under the sparse NB described in Section 4.3. In the experiments, the performance rates of the classifier shall be estimated according to an 10 runs 10−fold cross validation procedure. At each fold, the dataset is split into three sets, the so-called training, validation and testing sets. A tenth of the dataset is used as testing set, and the remaining nine tenths are for training set 2 3×9 10and validation set 1 3×9 10. Steps 2, 3(a) and 3(b) of Algorithm 4 are implemented on the training set. The different classifiers built in this way are compared according to their performance results (step 3(c)) on the validation set. The classifier (combination of features) with the highest performance measure on the validation set is chosen, and its average performance rates on the testing set are reported. Special emphasis will be made on the performance behavior and sparsity of the solutions of the proposed method. 4.4.1 Parameters setting The probability distribution for the features conditioned to the class Xi|Ck, for i= 1, . . . , p, k= 1, . . . , K, needs to be selected. It is well-known in the literature that the performance of the NB classifier improves when features are categorized using any discretization method [Liu et al., 2002; Boullé, 2004; Boullé, 2006]. Therefore, instead of imposing a specific probability distribution (such as the Gaussian), we adopted the discretization method based on an entropy criterion (see Dougherty et al. [1995]) and used the mdlp routine [Fayyad and Irani, 1993] from the discretization package of R. Now, we discuss the choice of the parameters {C, S, q}and the performance criterion. Choice of the parameter C The value of C, which represents the number of cuts in the vertical axis of the dendrogram, is critical for a proper sampling. As a default value, we propose to select the points of the 4.4. Numerical Illustrations 73 grid where features are clustered. When the routine hclust of R is used to generate the dendrogram (as in this work), one has C=p−1, where pis the number of features. In addition, hclust specifies where to make the cuts. However, a large value of Cmay slow down the execution of the algorithm notably and, on the other hand, it may lead to overfitting. For this reason, Cwill be defined as min{p−1,100}. As will be seen next, in Section 4.4.2, such a choice yields a right balance between the performance and the computational time for the considered datasets. Choice of the parameter S Regarding the value of S, which represents the maximum total of combinations for each cut, we tested several possible values for this parameter (see Appendix C), and settled on the final choice S= 25. Note that under the previous choices of Cand Sa total of max{25 ×(p− 1),25 ×100}combinations of features will be evaluated under the proposed sparse NB in contrast to the total number of possible combinations, equal to 2p−1. Choice of the parameter q Small values of qare associated with more sparsity (since fewer variables would be included in the combinations to be examined). Therefore, qshould be selected in such a way that it provides a compromise between the classifier’s performance and the sparsity of the solution. In particular, in Appendix C, different experiments to evaluate how the choice of this parameter affects the results can be found. Here, the selection of this parameter has been addressed according to the dependence matrix M, which is defined in Section 4.3. In particular, when 20% of the matrix elements are higher than 0.1(that is, from moderate to strong dependence cases), we fix q= 0.4(which implies a sparser solution). Otherwise (few dependent features), q= 0.6will be set. 4.4.2 Simulation study In this section we analyze how the performance of the sparse NB varies as dependence among features increases. In particular, we simulate data following Witten et al. [2014] and according to the model y=Xβ+εwith p∈ {100,200,300,400,500}. The errors ε1, . . . , εnare iid from aN(0,2.52)distribution. The observations (rows of X) are iid from a Np(0,Σ) distribution, where Σis a p×pblock diagonal matrix, with elements as follows: Σij =           1if i=j, ρif i≤p 4, j ≤p 4, i 6=j, ρif p 4+ 1 ≤i≤p, p 4+ 1 ≤j≤p, i 6=j, 0otherwise 74 Chapter 4. Variable selection for Naïve Bayes classification We explored various values of ρ, ranging from 0.1to 0.9. Furthermore, βi∼Unif[0.9, 1.1] for 1≤i≤ bp 4cand βi∼Unif[-1 3−0.1,−1 3+ 0.1] otherwise. In other words, there are two sets of p 4and 3p 4correlated features, respectively, and all the features are associated with the response. Finally, two classes are defined according to the sign of yn,n= 1,...,2000. The results in Table 4.2 have been obtained using the Mutual Information dependence measure (MI), S= 25 and values of qfixed as in Section 4.4.1. Moreover, the performance measure considered for these simulated experiments is the accuracy, and its average performance rates as in (4.3) on the testing set are reported in Table 4.2. Table 4.2: Average accuracy and sparsity (10 runs 10-fold CV) for simulated datasets. pMethod ρ= 0.1ρ= 0.5ρ= 0.7ρ= 0.9 ACC Sparsity ACC Sparsity ACC Sparsity ACC Sparsity 100 Sparse NB 88.20 66 93.30 37 93.20 24 95.05 21.30 Classic NB 89.50 100 87.05 100 87.10 100 86.95 100 200 Sparse NB 90.10 113.5 93.40 52.50 95.70 30.50 96.55 20.10 Classic NB 90.10 200 87.30 200 87.65 200 86.80 200 300 Sparse NB 89.85 168.60 92.40 79.70 94.70 37.30 95.95 18.20 Classic NB 90.15 300 87.00 300 87.15 300 87.85 300 400 Sparse NB 91.90 216.30 91.45 94 92.35 46.80 93.65 17 Classic NB 89.65 400 87.25 400 87.25 400 87.50 400 500 Sparse NB 91.90 283.20 90.70 102.70 93.85 29.60 91.90 5.90 Classic NB 90.15 500 87.05 500 87.45 500 87.55 500 Some conclusions can be drawn. On the one hand, in terms of sparsity levels, the sparse NB returns better results in the presence of moderate to strong dependence cases. For datasets where the dependences among features are weak, ρ= 0.1, our sparse strategy is able to remove around one third of the total number of variables whereas, in some cases, the ACC is slightly reduced with regards to the classic NB. While ρincreases, our proposal is able to significantly reduce the number of variables considered, achieving better ACC results than the classic NB, as the curse of dependency is alleviated. On the other hand, Figure 4.3 reports the logarithm of the average user times (in seconds) when the sparse NB is run on Intel(R) Core(TM) i7-7500U CPU at 2.70GHz 2.90GHz with 8.0 GB of RAM. The X-axis shows the ρvalues whereas each line represents the number of variables of the dataset (p). Overall, for weak dependences among features, the behaviour of runing time is monotonous respect to the number of variables, but this changes when ρincreases. 4.4.3 Datasets and benchmark approaches The so-called Breast Cancer Wisconsin (Diagnostic) Data Set,Wine Data Set,Mushroom, Waveform Database Generator Data Set (version 2), ISOLET Data Set,Multiple Features Data Set,SPECTF Heart Data Set,German Credit,Page Blocks Classification Data Set and Statlog (Australian Credit Approval) shall be considered. They are described in Table 4.3, whose first three colums report the dataset name, the number of instances and the class split. The number of 4.4. Numerical Illustrations 75 Figure 4.3: Scalability continuous variables (L) and categorical variables (L0) are presented in the last two columns. Three of the ten datasets, SPECTF,german and page blocks, are unbalanced datasets, due to the very different sizes of the classes. Table 4.3: Datasets description Name Instances Class split in % L L’ breast cancer 569 63(Benign)/37(Malignant) 30 0 wine 178 33(Class 1)/40(Class 2)/27(Class 3) 13 0 mushroom 8124 51.8(edible)/48.2(poisonous) 0 22 waveform 5000 33.33(Class 0)/33.33(Class 1)/33.33(Class 2) 40 0 ISOLET 7797 26 equiprobable classes (0.04) 617 0 Multiple Features 2000 9 equiprobable classes (0.11) 649 0 SPECTF 267 79(Abnormal)/21(Normal) 44 0 german 1000 70(Class 1)/30(Class 2) 7 13 page blocks 5473 90(Negative)/10(Positive) 10 0 australian 690 55.5(Negative)/44.5(Positive) 6 8 We aim to compare the novel method with alternative, well-known strategies for feature selection. In this study, we focus on techniques which perform hard variable selection and, in consequence, feature weighting approaches as in Jiang et al. [2019] have not been considered here. Specifically, we selected one filter and one wrapper that are well referenced in the literature and that can be easily adapted to the NB classifier to make a fair comparison. Our choice was the filter CFS and the wrapper Boruta, both introduced in Section 1.1.1. These methods are widely spread and can be computed by the routines cfs and Boruta, from R packages FSelector and Boruta, respectively. In order to adapt the wrapper Boruta to the NB classifier, we have used the function filterVarImp in R package caret as the function that returns the importance of the attributes, instead of the default getImpRfZ, which is based on the Random Forest classifier. It is important to highlight that the time limit is not an input parameter of the cfs and Boruta routines and therefore, differences in the computational 76 Chapter 4. Variable selection for Naïve Bayes classification costs were found (to be discussed later). Apart from the previous feature selection methods, that can be applied to any classifier, there are works that specifically deal with variable reduction for the NB. In particular, Boullé [2007] proposes a straightforward Bayesian modern-style approach, the MAP Approach for Variable Selection (noted SNB(MAP)), where the conditional probabilities are formulated according to p(Ck|x)∝π(Ck) p Y i=1 p(xi|Ck)ai, k = 1, . . . , K. (4.6) In Eq. (4.6), the values {ai}p i=1 are either 1or 0, depending on whether or not feature iis included in the model. Then, the posterior distribution of the different models (resulting from different choices of {ai}p i=1) is evaluated using a shrinkage prior so that parsimoniuos models are favoured. In the same paper, a search heuristic that performs a fast forward backward selection is described and therefore, it has been implemented in this work to run the different experiments. However, when the number of variables increases, note that the time required to run this method is excessive. For that reason, a time limit of eight hours for the two biggest considered real datasets (ISOLET and Multiple Features) was fixed. Finally, we have also compared with the Lasso approach for classification (see Vincent and Hansen [2014]), whose goal is precisely the same: obtain good classification performance while selecting few features. The routine fit in R package msgl has been used. Next, we will break the results down depending on the datasets are balanced or unbalanced. As we will show below, if necessary and motivated by the properties of the dataset, our proposal can be easily adapted in terms of the performance criterion to be optimized and the required constraints on groups of interest. To make a fair comparison, we do not impose any additional constraints and therefore only the performance criterion will change accordingly throughout these sections. However, an illustrative example where constraints are imposed is also included at the end of Section 4.4.5. 4.4.4 Results for balanced datasets For comparison purposes, consider the same parameters setting than in Section 4.4.2, where for waveform dataset, qis equal to 0.6and, for the remaining balanced databases, q= 0.4. We next analyze the performance and sparsity of the method, as well as the running times. The average accuracy, number of variables in the selected combinations and the CPU time in seconds for 1 fold-CV execution are shown by Figure 4.4. Moreover, a comparison between the sparse NB with the above-mentioned feature selection methods is made. The results under the classic NB, CFS, Boruta, SNB(MAP) and Lasso methods are also shown. Several conclusions can be drawn at this point. Note that the performance rates under the sparse NB are comparable to the classic NB using between a half and one third of the variables, except for waveform. As commented before, the novel approach is intended to 84 Chapter 5. Constrained Naïve Bayes with application to unbalanced data classification The approach introduced in Chapter 2, one overall criterion is optimized while constraints to demand admissible values for the efficiency measures under consideration are introduced in the model, is explored for improving the NB performance in the classes of most interest to the user. It will be seen that unlike the traditional NB, which is a two-step classifier (estimation first and classification next), the novel approach integrates both stages. In particular, maximum likelihood estimates are replaced here by constrained maximum likelihood estimates, where the constraints control the Recall values of the classes of interest. 5.1 Introduction In this chapter we propose a novel way of controlling misclassification rates, that do not call for using misclassification costs which may be hard to choose and are not usually given [Sun et al., 2007, 2009]. In particular, a new version of the NB is obtained by modeling performance contraints where the Recall (proportion of instances of a given class correctly classified) for the classes of interest is forced to be lower-bounded by certain thresholds. In this way, the user is allowed to assign different importance to the different classes according to her preferences. For example, in the previously considered breast cancer dataset, it may be desiderable to increase the Recall for the Malignant class, which is equal to 88.41. As it will be shown in Section 5.3, for this case such rate can be increased up to 91.88. Other example where performance constraints are useful is when fair classification is a requirement as a social criterion, and then the sensitive groups should be protected to avoid the discrimination against race, or other sensitive data [Romei and Ruggieri, 2014]. Acceptable values for the Recall of groups at risk could be fixed via the proposed method in this work. The problem of cost imbalance has been addressed in the literature from two different perspectives: Data-Level techniques and Algorithm-Level approches, see Leevy et al. [2018]. Whereas the former include data sampling methods and feature selection, the latter encompass cost-sensitive and hybrid/ensemble methods which adapt the base classifier to overcome the imbalance. Particularly, our approach can be seen as a cost-sensitive method. Cost-sensitive approaches have already been considered in the literature for well-known classifiers, such as those introduced in Section 1.1.2. However, there is a lack of methodology that allows the user to have full control on the different performance measures of interest at the same time, which is the what we explore in this work. This chapter is organized as follows. In Section 5.2 the notation is introduced and the proposed version of constrained NB (CNB from now on) is described. Section 5.3 illustrates the usefulness of our novel approach. Eight real databases with different sampling properties are thoroughly analyzed, and a detailed discussion concerning the Recall values of the proposed approach compared with the classic NB is given. Some conclusions are considered in Section 5.4. 5.2. The constrained Naïve Bayes 85 5.2 The constrained Naïve Bayes In our approach, the estimation is performed by solving a constrained maximum likelihood estimation problem, constraints being related with thresholds on the Recall values for different classes. The aim of this section is to describe the associated optimization problem. As a result, a computationally tractable classifier that allows the user to control its performance is obtained. 5.2.1 Preliminaries on NB classification: notation Consider again the random vector (Y, X), where X= (X1, . . . , Xp)contains pfeatures and Y identifies the class label. Assume that we have a classification problem with Kclasses. Then, for each class k∈ {1, . . . , K}, let πkdenote the prior probability of the class, and assume that Xj|(Y=k)has a probability density function fθjk (x), where θjk ∈Θjk. For k= 1, . . . , K, define θk= (θ1k, . . . , θpk). Let x= (x1, . . . , xp)be a new observation. The aim is to label it on one of the Kclasses. Then, under the 0-1 loss function, Bayesian Decision Theory establishes that xis classified in the most probable class according to the conditional distribution. The estimation of the associated parameters may be cumbersome if the number of features pis large. However, the use of the Bayes theorem, in addition to the assumption of independence (conditioned to the class) ease the previous estimation process. As it is well known, the latter assumption implies that the joint density function can be expressed as f(x1, . . . , xp, k) = πkf(x1, . . . , xp|k) =πkfθk(x) =πk p Y j=1 fθjk (xj), and thus the estimation process is reduced to estimate the parameters of each marginal distribution. Then, the NB classifier performs by assigning xto class ksatisfying πk p Y j=1 fθjk (xj)≥πl p Y j=1 fθjl (xj)∀l= 1, . . . , K. (5.1) Given a training sample of size n1,(k1,x1),...,(kn1,xn1), then θ= (θ1,...,θK)is estimated in NB via maximum likelihood [Hogg et al., 2005], and therefore computed as the solution of the optimization problem: max θ n1 X i=1 log fθki(xi)(5.2) Therefore, the classic NB can be seen as a two-step classifier, where the model parameter is first estimated as ˆ θfrom a training sample, and then (5.1) is applied under θ=ˆ θ. 86 Chapter 5. Constrained Naïve Bayes with application to unbalanced data classification 5.2.2 A novel formulation with performance constraints In order to calibrate the performance of a classifier, many measures have been defined in the literature, see Sokolova and Lapalme [2009]. In particular, the so-called Recallk, for k= 1, . . . , K, is defined as the sample fraction of individuals in class kwhich are correctly classified. Given a validation sample of size n2, where n2=Pkn2,k and n2,k is the size of class kin such a validation sample, (k, x(k) 1),...,(k, x(k) n2,k ), then the Recall for class kcan be expressed as functions of ˆ θ, Recallk(ˆ θ) = 1 n2,k n2,k X i=1 Ck(ˆ θ,x(k) i), k = 1, . . . , K, (5.3) where Ck(ˆ θ,x(k) i) =      1if the individual x(k) iis classified in class k 0otherwise. (5.4) Unlike the classic NB, based on a two-step approach, the CNB proposed in this work integrates the performance of the classifier (according to expression (5.3)) within the estimation step. In particular, the pursued aim is to estimate θas the solution of an optimization problem where the objective function is given using a training sample of size n1as in (5.2) and, to prevent overfitting, constraints on (5.3) are imposed on an independent sample (validation set) of size n2=PK k=1 n2,k, max θ n1 X i=1 log fθki(xi) s.t. 1 n2,k n2,k X i=1 Ck(θ,x(k) i)≥αk, k = 1, . . . , K (CNB) In the previous CNB optimization problem, αk∈(0,1) is a threshold, a lower-bound value close to 1, for k= 1, . . . , K, which is fixed by the user according to her requirements about the classification in the different classes. From the point of view of optimization, we assume that the function fθkiis smooth with respect to the parameter θki. Regarding the constraints, they are not smooth and therefore, gradient methods cannot be applied in order to solve Problem (CNB). This fact makes the resolution of (CNB) to be slow, especially for large datasets. However, a proxy version of (CNB) can be written in a more tractable way if the constraints are reformulated in terms of smooth functions as e Ck(θ,x(k);λ) = K Y l=1,l6=k F(ykl(θ,x(k)); λ),(5.5) 5.2. The constrained Naïve Bayes 87 where F(y;λ) = 1 1+e−λy is the sigmoid function and ykl(θ,x) = πk p Y j=1 fθjk (xj)−πl p Y j=1 fθjl (xj).(5.6) On the one hand, from the definition of the sigmoid function, it can be seen that limλ→∞ e Ck(θ,x(k);λ) = Ck(θ,x(k)), since for large λvalues F(ykl(θ,x(k)); λ)will only take the values 0 or 1 depending on the sign of ykl(θ,x(k)). Then, λis a hyperparameter big enough so that Cand e Care as close as possible. On the other hand, the reason why we use the product function to define e Cis explained below. Note that if any class lhas associated a density much greater than class k, then ykl will take a large negative value which makes F(ykl(θ,x(k)); λ)close to 0 and therefore e Ck(θ,x(k);λ)will also be close to 0. From the previous discussion, a differentiable version of the CNB problem is obtained as max θ n1 X i=1 log fθki(xi) s.t. 1 n2,k n2,k X i=1 ˜ Ck(θ,x(k) i)≥αk, k = 1, . . . , K. (SCNB) The smooth formulation (SCNB) can be solved using efficient solvers for nonlinear constrained programming (see, e.g. Birgin and Martínez [2008]). From now on, we refer to (SCNB) as our optimization problem. Two important remarks need to be made at this point. The first one regards the feasibility of the (SCNB). In a real application, threshold values α1, . . . , αKhave to be fixed. As a first option, they could be fixed by the user according to her demand, but it might be the case that (SCNB) is unfeasible. For that reason, we propose a procedure for determining the thresholds in such a way that (SCNB) is always feasible. If we consider a dataset with Kdifferent classes, let θ∗be the model parameter associated with (5.2) and k0be the critical class or the class where the method performs the worst. Suppose that the aim is to improve the Recall for such class k0, say αk0=1 n2,k0 n2,k0 X i=1 ˜ Ck0(θ∗,x(k0) i)+∆, with ∆>0. Then, in order to know the maximum threshold τfor the other classes k6=k0, k ∈ {1, . . . , K}, the next optimization problem can be solved: 88 Chapter 5. Constrained Naïve Bayes with application to unbalanced data classification max θ,τ τ s.t. 1 n2,k0 n2,k0 X i=1 ˜ Ck0(θ,x(k0) i)≥1 n2,k0 n2,k0 X i=1 ˜ Ck0(θ∗,x(k0) i)+∆ 1 n2,k n2,k X i=1 ˜ Ck(θ,x(k) i)≥τ, ∀k6=k0. This way we search the estimates θsuch that in the relevant class k0the Recall is improved in at least ∆with respect to the Recall in the traditional Bayes estimate and maximize the minimum Recall in the remaining classes. The second remark concerns the solutions of (SCNB), which are not maximum likelihood estimates any more, but maximum constrained likelihood estimates instead. On the contrary, the problem yields a solution with the highest sample likelihood fulfilling the constraints on performance on the independent sample. Up to our knowledge, this is a breaking approach that has never been considered in NB models. 5.3 Numerical results In this section, eight data sets from the UCI Machine Learning Repository and KEEL open source [Alcalá-Fdez et al., 2011, 2009] diverse, in both in the number of classes, sizes and imbalance ratio shall be analyzed. The description of the datasets can be found in Section 5.3.1 and the numerical experiments and obtained results will be considered in Section 5.3.2 and 5.3.3, respectively. 5.3.1 Datasets The datasets breast cancer,SPECTF,page-blocks,abalone,yeast,Satimage, RCV1 and letter will be considered. From all the available versions of the datasets, we have chosen those described in Table 5.1. The colums report the dataset name, the number of instances and features and finally, the class split of the eight considered datasets (page-blocks, abalone,yeast,Satimage and RCV1 are significantly more unbalanced than the other three). 5.3.2 Design of experiments Probability distributions setting and resolution of the optimization problem As comented in Section 5.2.1, a probability model needs to be selected for the features conditioned to the class. If the feature is continuous, in this chapter we will assume the normal 5.3. Numerical results 89 Name Intances Features Class split (%) breast cancer 569 30 (Benign,Malignant) (63, 37) SPECTF 267 44 (Abnormal,Normal) (79, 21) page-blocks 5473 10 (text,horiz. line,graphic,vert. line,picture) (89.8, 6, 0.5, 1.6, 2.1) abalone 4177 8 (1-5,6-10,11-15,16-29) (4.52, 28.39, 6.25, 60.83) yeast 1484 8 (CYT,EXC,ME1,ME3,MIT,NUC,POX,VAC) (32.42, 2.45, 3.08, 11.41, 17.09, 30.04, 1.40, 2.10) Satimage 6435 36 (1,2,3,4,5,7) (23.82, 10.92, 21.10, 9.73, 10.99, 23.43) RCV1 18758 21531 (C15,CCAT,E21,ECAT,GCAT,M11) (23.70, 20.12, 5.73, 9.54, 22.43, 5.61) letter 20000 16 From Ato Z(26 classes) Equally distributed Table 5.1: Datasets description distribution. From the point of view of the optimization, (SCNB) will be solved using solvers for smooth optimization. In particular, auglag and mma functions from R package nloptr will be used in this work to obtain all numerical results. Estimation of the performance rates The performance of the proposed classifier will be estimated using a classic 10 Monte-Carlo cross validation, since in this chapter we deal with highly unbalanced datasets. The dataset will be split into three sets, the so-called training set, validation set and testing set. One-third of the dataset is used as testing set, and the remaining two-thirds for training set 2 3×2 3and validation set 1 3×2 3. As explained in Section 5.2, the objective function will be optimized on the training set while the constraints will be evaluated on the validation set. Once the SCNB problem is solved, Recall values are estimated on the testing set. It must be highlighted that at each run, the training sample is built in a stratified way so that the proportion of samples per class is similar to the proportions depicted by Table 5.1. Finally, regarding the hyperparameter λ, after an extensive simulation study considering a wide grid of values, the choice λ= 23is set in the experiments since it provides a good match between Cand ˜ Cas in (5.4) and (5.5). Pre-processing for large datasets Problem SCNB turns out computationally costly for large datasets as the considered RCV1 dataset. As it is common in the literature (see Leevy et al. [2018] and references therein), we suggest to pre-process such datasets in a way that irrelevant variables are removed in a first step previous to the resolution of (SCNB). Specifically, in this work the importance of the predictor variables composing RCV1 were measured using the R function information.gain from 90 Chapter 5. Constrained Naïve Bayes with application to unbalanced data classification Figure 5.1: Boxplot of the importance of the variables for RCV1 dataset. FSelector. Figure 5.1 shows the boxplot of the importance associated with the variables of the RCV1 dataset. As it can be observed, most of them have an associated importance close to 0 and, then, only 392 of the total are going to be kept when solving (SCNB). The choice of thresholds In order to select the threshold values αkin Problem (SCNB), the classic NB classifier (5.2) was first run. Table 5.2 shows the Recall estimates for each class. In particular, from Table 5.2 and in order to improve the rate of the classes where the classic NB performs the worst, the set of thresholds to be tested in the numerical experiments shall be given by Table 5.3. Note that, according to the previous values, not only better rates for the classes with the worst associated Recall are imposed, but also, admissible values for the rest of classes are fixed. 5.3.3 Results The estimated rates are reported in Tables 5.4-5.11. The first row shows the results for the classic NB, when no thresholds are imposed. The first column shows the imposed thresholds for the Recall of each class, whereas the column and thresholds in bold correspond to the class where the classic method performs the worst. For example, in Table 5.5, it is required that the Recall of Normal class is at least 80, while over the Abnormal class the threshold varies from 67 to 70. The remaining columns provide the average Recall values measured on the test set. 5.3. Numerical results 91 Name Recall Classic NB breast cancer (Benign,Malignant) (96.15, 88.41) SPECTF (Abnormal,Normal) (65.40, 91.29) page-blocks (text,horiz. line,graphic,vert. line,picture) (90.61, 68.06, 67.14, 94.14, 40.11) abalone (1-5,6-10,11-15,16-29) (92.58, 61.23, 22.09, 53.97) yeast (CYT,EXC,ME1,ME3,MIT,NUC,POX,VAC) (2.37, 62.73, 73.57, 19.81, 54.13, 48.79, 53.33, 31.11) Satimage (1,2,3,4,5,7) (79.79, 89.59, 89.60, 66.70, 73.52, 74.21) RCV1 (C15,CCAT,E21,ECAT,GCAT,M11) (87.55, 5.53, 74.99, 23.16, 75.12, 88.14) letter From Ato Z(26 classes) First row in Table 5.11. The letter Sis the worst classified. Table 5.2: Average Recall of classic NB (10 Monte-Carlo cross validation) breast cancer: 85 (Benign), 90/91.5/93/94.5 (Malignant) SPECTF: 67/68.5/70 (Abnormal), 80 (Normal) page-blocks: 80 (text), 58 (horiz. line), 57 (graphic), 84 (vert. line), 43/44.5 (picture) abalone: 70 (1-5), 50 (6-10), 24/26/28/30 (11-15), 40 (16-29) yeast: 5/10/15 (CYT), 50 (EXC), 60 (ME1), 10 (ME3), 40 (MIT), 30 (NUC), 40 (POX), 20 (VAC) Satimage: 70 (1), 80 (2), 80 (3), 68/70 (4), 60 (5), 60 (7) RCV1: 60 (C15), 6/7/8 (CCAT), 60 (E21), 10 (ECAT), 60 (GCAT), 60 (M11) letter: 24.5/26/27.5 (S), Twenty-five percent less than the results of the classic NB (The rest of classes) Table 5.3: Tested thresholds As expected, the results under the constrained NB version differ from the results provided by the classic NB. For example, for the breast cancer dataset, the Recall values under the classic NB are 96.15 and 88.41, for Benign and Malignant class, respectively (Table 5.2). As commented in Section 5.1, it may be of interest to increase the Recall of the Malignant class. According to Table 5.4, if the minimum 94.50 is imposed for the Malignant class, the final rates change from 88.41 to 91.88. It is important to highlight two different facts concerning the previous results. First, note that a better rate for the Malignant class has been obtained, but at the expense of slightly decreasing the rate of the Benign class. Second, note that even though a rate equal to 94.50 was imposed, such value was not finally obtained, but a smaller one (91.88). This is not surprising, since the constraints are imposed for one sample, and tested on an independent set. From the results shown in Table 5.4-5.11, it can be concluded that the proposed approach allows the user to control the Recall values in such a way that the classes where the classic method performs the worst can be improved, even switching to the classes where the best rates are achieved (e.g., in breast cancer dataset) and changing thus the natural tendency of the classifier. Note that among the possible non-dominated solutions shown for each dataset, the 92 Chapter 5. Constrained Naïve Bayes with application to unbalanced data classification user could choose according to her interest and to what she is willing to lose in the less critical classes. Thresholds (Benign/Malignant) Recall Benign Recall Malignant Classic NB 96.15 88.41 85.00/90.00 94.62 88.70 85.00/91.50 94.10 90.00 85.00/93.00 92.39 91.59 85.00/94.50 91.45 91.88 Table 5.4: Average Recall values of SCNB (10 Monte-Carlo cross validation) for breast cancer Thresholds (Abnormal/Normal) Recall Abnormal Recall Normal Classic NB 65.40 91.29 67.00/80.00 67.09 88.72 68.50/80.00 67.37 88.08 70.00/80.00 69.63 87.08 Table 5.5: Average Recall values of SCNB (10 Monte-Carlo cross validation) for SPECTF Thresholds (text/horiz. line/ graphic/vert. line/picture) Recall text Recall horiz. line Recall graphic Recall vert. line Recall picture Classic NB 90.61 68.06 67.14 94.14 40.11 80.00/58.00/57.00/84.00/43.00 89.34 68.97 70.23 93.23 43.58 80.00/58.00/57.00/84.00/44.50 90.15 69.21 73.56 94.14 44.28 Table 5.6: Average Recall values of SCNB (10 Monte-Carlo cross validation) for page-blocks Thresholds (1-5/6-10/11-15/16-29) Recall 1-5 Recall 6-10 Recall 11-15 16-29 Classic NB 92.58 61.23 22.09 53.97 70.00/50.00/24.00/40.00 92.58 60.56 24.65 53.87 70.00/50.00/26.00/40.00 92.58 60.05 25.47 53.85 70.00/50.00/28.00/40.00 92.58 59.36 26.28 53.81 70.00/50.00/30.00/40.00 92.58 58.31 28.26 53.75 Table 5.7: Average Recall values of SCNB (10 Monte-Carlo cross validation) for abalone Thresholds (CYT/EXC/ME1/ME3/ MIT/NUC/POX/VAC) CYT EXC ME1 ME3 MIT NUC POX VAC Classic NB 2.37 62.73 73.57 19.81 54.13 48.79 53.33 31.11 5.00/50.00/60.00/10.00/40.00/30.00/40.00/20.00 3.88 61.82 74.29 24.34 56.75 47.45 53.33 35.56 10.00/50.00/60.00/10.00/40.00/30.00/40.00/20.00 6.18 58.18 75.00 24.91 56.50 47.52 55.00 35.56 15.00/50.00/60.00/10.00/40.00/30.00/40.00/20.00 14.87 58.18 75.71 28.68 58.63 39.36 53.33 33.33 Table 5.8: Average Recall values of SCNB (10 Monte-Carlo cross validation) for yeast 6.2. Preliminaries on elliptical and NI distributions 99 [1986], the maximum-likelihood estimators of the model parameters are obtained and in Kano et al. [1993]; Berkane et al. [1994], a pseudo-maximum likelihood approach is considered. On the other hand, the Bayesian perspective has also been explored for the estimation of elliptical distributions, see for example Osiewalski and Steel [1993]; Fang and Li [1999]; Branco et al. [2000]; Niekerk et al. [2015]. Finally, Maruyama and Seo [2003] undertakes a momentsmatching approach. The semi-parametric estimation approach has also been explored. For example, Liebscher [2005] used kernel density estimation to provide a classical, semi-parametric estimator. Also, Fortunati et al. [2020] undertakes the estimation of the elliptical distribution in semi-parametric fashion which results in a robust, efficient estimator. In this chapter we adopt the Bayesian paradigm and propose a semi-parametric approach based on Dirichlet processes [Ferguson, 1973]. The stochastic representation of the elliptical distribution discussed in Section 6.2 will play a key role for the inference approach. The remainder of the chapter is structured as follows. In Section 6.2, we formally define elliptical as well as NI distributions, and their main properties. Then, in Section 6.3 we illustrate how semi-parametric, Bayesian inference can be carried out for both NI and elliptical distributions. Our approach is illustrated in detail with simulated and a real dataset in Section 6.4. The chapter ends with some concluding remarks in Section 6.5. 6.2 Preliminaries on elliptical and NI distributions According to Cambanis et al. [1981] the elliptical contoured distributions is defined as follows. As previously, let Ydenote a k−dimensional random vector, µ∈Rkand let Σbe some k×kcovariance matrix. If the characteristic function of Y−µ,φY−µ(t)is a function of the quadratic form t0Σt, that is, if φY−µ(t) = φt0Σt, then, Yfollows an elliptical (contoured) distribution with location parameter µ, scale parameter Σand scalar function φ. Note that when φ(u) = e−u/2, then the normal multivariate distribution with parameters (µ,Σ) is recovered. The elliptical distribution admits a stochastic representation as proven in Schoenberg [1938]; Cambanis et al. [1981]. This is an advantageous result which allows to simulate in straightforward way from the elliptical distribution and provides a natural scheme for designing the estimation method (as will be detailed in Section 6.3). Prior to the stochastic representation, let us introduce some concepts in relation to the uniform distribution on the unit sphere in Rk. If Φk,k≥1is the class of all functions φ: [0,∞)→Rsuch that φ(||t||2),t∈Rk, is a characteristic function, then φ∈Φkif and only if φ(s) = Z∞ 0 Ωk(r2s)dF(r), s ≥0 100 Chapter 6. A Bayesian semi-parametric approach to normal/independent and elliptical distributions for some distribution function Fon [0,∞), where Ωk(||t||2),t∈Rk, is the characteristic function of a k−dimensional random vector Skwhich is uniformly distributed on the unit sphere in Rk. Remark 1. If Skis uniformly distributed on the unit sphere in Rk, then Sk=Z/|Z|, where Z∼N(0,Ik)and |Z|=qZ2 1+. . . +Z2 k. Theorem 3 (Schoenberg [1938]; Cambanis et al. [1981]).Ak×1random vector Yfollows an elliptical distribution with parameters µ,Σand φif and only if: Yd =µ+ASkU, (6.2) where Skis uniformly distributed on the (k-dimensional) hypersphere, Uis a non-negative, univariate random variable, independent of Skand Ais the unique, symmetric matrix, with non-negative diagonal elements such that AA = Σ. The derivation of the density function as in (6.1) from the stochastic representation (6.2) can be found in Cambanis et al. [1981] (section 4). For the purposes of this work, it is useful to introduce some reparametrizations of the elliptical distribution. Suppose that we have an elliptical variable defined as in Theorem 3. Then, we shall write: Yd =µ+ASk√Gwhere G=U2.(6.3) Given the stochastic formulation in (6.3), it is clear that an elliptical variable is alternatively defined by the location vector, µ, the scale matrix, Σ, and the distribution of the variable G, say FG(·). Therefore, from now on we shall refer to an elliptical variable satisfying (6.3) as Y∼EC(µ,Σ, FG). Note that we have G= (Y−µ)0Σ−1(Y−µ).(6.4) From Remark 1, we have that: Yd =µ+AZ/√W, (6.5) where W=|Z|2 G. For the general elliptical distribution as in (6.5), Wand Zare not necessarily independent. However, if this assumption is made, following Andrews and Mallows [1974], then Yis said to have a normal/independent (NI) distribution. A NI distribution is therefore defined by µ,Σ and the distribution of W, say FW(·), so we shall write Y∼NI(µ,Σ, FW)to represent a NI variable satisfying (6.5). In order to avoid identifiability problems is important to take into account that elliptical distributions are identified, up to a scale factor. If we consider an elliptical variable Y, defined as in (6.5), we can write A∗=bAand W∗=b2Wfor any non-negative constant, b, and then 6.3. Bayesian inference for NI and elliptical distributions 101 we have an equivalent stochastic representation as: Yd =µ+A∗Z/√W∗. Therefore, in order to ensure identifiability, throughout this work, we shall assume, for both the elliptical and NI models, that the diagonal elements of Σsum up to one, that is: k X i=1 Σii = 1,(6.6) where Σij is the (i, j)−th element of matrix Σ, for i, j = 1, . . . , k. Concerning the first two moments of elliptical distributions, it is known that if E(U)<∞, then E(Y) = µ. Also, if E(U2)<∞, then Σis proportional to the covariance matrix, E(Y−µ)(Y−µ)0=k−1E(U2)Σ. Finally, it is interesting to mention that if two random vectors are jointly elliptically distributed, then the conditional distribution of one given the other is elliptical. 6.3 Bayesian inference for NI and elliptical distributions In order to implement Bayesian inference for either the NI or the general, elliptical models, we first need to define prior distributions for the location vector, µand the scale matrix Σ. Here, we shall assume a normal prior distribution: µ∼N(m,V) where m= (m1, ..., mk)0is a constant vector and Vis a k×kvariance-covariance matrix. In practice, we set mto be a zero vector and V= 1000Ik, where Ikis the k×kidentity matrix, in order to provide a proper but relatively uninformative prior. In order to define a prior distribution for Σ, we need to recall the identifiability restriction in (6.6). To do this, we shall first define: Σ = diag(σ)Rdiag(σ), where diag(σ)is a diagonal matrix with diagonal elements σi=√Σii for i= 1, ..., k and R is a matrix with elements Rij =Σij pΣiiΣjj =Σij σiσj . Now, we can define an implicit prior for Σby setting prior distributions for σ2= (σ2 1, ..., σ2 k)0= 102 Chapter 6. A Bayesian semi-parametric approach to normal/independent and elliptical distributions (Σ11, ..., Σkk)0and R. Firstly, we assume that σ2has a Dirichlet distribution σ2∼Dir(aσ), where aσ= (aσ1, ..., aσk)0. In practice, we set aσi = 1 for i= 1, ..., k. For R, we follow an approach of Zhang et al. [2006] and use a parameter extended inverse Wishart distribution. Thus, we set M=D1 2RD1 2 where Dis a diagonal matrix. Then we assume that Mhas an inverse Wishart prior distribution, that is M ∼ IW (k+ 1,Ik). Under this structure, it can be shown that the implicit prior for Ris the marginal uniform distribution of Barnard et al. [2000]. Finally, it is necessary to define prior distributions for the distributions, FW, in the case of the NI model or FG, in the elliptical case. Instead of using a parametric model, here we shall use Dirichlet process mixtures of gamma distributions as in Hanson [2006]. We outline the prior set up and the posterior inference procedures for the NI and elliptical models respectively in the following subsections. 6.3.1 Inference for the NI distribution Assume that we observe a sample y1,...,ynfrom a NI distribution with stochastic representation as in (6.5), Yd =µ+AZ/√W, where Zand Ware independent. As Wis non-negative, a natural approach, which we shall take here, is to follow Hanson [2006] and model Wusing a Dirichlet process mixture of gamma distributions. Therefore, we assume that W|λ, η ∼Ga (λ, η) λ, η|F∼F(6.7) F∼DP(αF0), where DP(·)represents a Dirichlet process with paramaters α > 0and F0(·,·)being a prior mean distribution for F(·,·). Also following Hanson [2006], we shall assume that F0is such that the associated density function is a product of exponential distributions f0(λ, η|aλ, aη) = aλe−aλλaηe−aηηfor λ, η > 0,where aλ, aη>0, 6.3. Bayesian inference for NI and elliptical distributions 103 and we shall place hyperprior distributions α∼Ga(aα, bα) aλ∼Ga(bλ, cλ) aη∼Ga(bη, cη), where aα, bα, bλ, cλ, bη, cη>0. In the practical examples later, we set aα=bα= 2 and bλ=cλ=bη=cη= 0.001. Under this model, we have that integrating out over W, the implied distribution of Yis an infinite mixture of Student’s t distributions, all with location parameter µand scale parameter proportional to Σ. Given this prior formulation and supposing that we observe a sample, y1, ..., yn, then inference can be carried out using a Gibbs sampling scheme based on successively sampling the conditional posterior distributions. The relevant formulae are given below. Firstly, the conditional posterior distributions of µ,σ2,D,Rare as follows: µ|Σ,w,y∼N(m∗,V∗)where m∗=V∗V−1m+nΣ−1wywhere wy=1 n n X i=1 wiyi V∗=V−1+n¯wΣ−1−1where ¯w=1 n n X i=1 wi f(σ2|µ,R,w,y)∝ n Y i=1 Nyi|µ,1 wi Σf(σ2) f(R,D|µ,σ2,w,y)∝ n Y i=1 Nyi|µ,1 wi Σf(R,D) where, throughout, Σ = diag(σ)Rdiag(σ),w= (w1, ..., wn)0,N(yi|µ,1 wiΣ) refers to a multivariate normal density function with parameters µand 1 wiΣevaluated at yi, whereas f(σ2)and f(R,D)are the prior distributions. The conditional posterior of σ2can be sampled via a Metropolis-Hastings pass, using for example an adaptive logit sampler as in Director et al. [2017]. In order to sample the distribution of D,R, we follow Fang and Li [1999] and use a Wishart candidate for the composed matrix M=D1 2RD1 2. The conditional posterior distribution of Wiis directly evaluable as Wi|µ,Σ,yi∼Ga k 2,(yi−µ)0Σ−1(yi−µ). Now, the conditional posterior distributions of the Dirichlet process parameters follow from Hanson [2006]. Firstly, we can generate values from the conditional posterior of λi, ηigiven λ−i,η−i, aλ, aη, wiwhere λ−i= (λ1, ..., λi−1, λi+1, ..., λn)0and similarly for η. We have 104 Chapter 6. A Bayesian semi-parametric approach to normal/independent and elliptical distributions then that: (λi, ηi)=(λj, ηj)with probability pjfor j6=iin 1, ..., n, where pj∝Ga(wi|λj, ηj)and Ga(·|·,·)is a gamma density function. With probability pi, where pi∝αaλaη wi(wi+aη)(aλ−log wi+ log(wi+aη))2, it follows λi∼Ga(2, aλ−log wi+ log(wi+aη)) ηi|λi∼Ga(λi+ 1, wi+aη). Secondly, the hyperparameters are: aλ|λ∼Ga  bλ+m∗, cλ+ m∗ X j=1 λ∗ j  aη|η∼Ga  bη+m∗, cη+ m∗ X j=1 η∗ j , where there are m∗unique values, say {(λ∗ 1, η∗ 1), ..., (λ∗ m∗, η∗ m∗)}. Finally, values from the distribution of αcan be generated following Escobar and West [1995] by introducing a further latent variable, θ, such that: θ|α, n ∼Beta(α+ 1, n) α|θ, m∗∼(Ga(aα+m∗, bα−log θ)with probability aα+m∗−1 aα+m∗−1+n(bα−log θ) Ga(aα+m∗−1, bα−log θ)otherwise. 6.3.2 Inference for the elliptical distribution Now assume that we observe a sample y1,...,ynfrom (6.3). Then, we shall assume a Dirichlet process mixture prior: G|λ, η ∼Ga (λ, η) λ, η|F∼F F∼DP(αF0), where we shall use the same prior modeling for α, F0as previously in (6.7). In order to undertake inference, we first use (6.5), so that Yd =µ+AZ/√W 6.3. Bayesian inference for NI and elliptical distributions 105 where W=|Z|2/G. In addition, note that Xd =|Z|2∼χ2 k. Also, we have that: f(wi|λi, ηi) = Z∞ 0 gifXi(wigi)fGi(gi|λi, ηi)dgi =1 Beta k 2, λi wi 2 k 2−1ηλi i 2k 2wi 2+ηik 2+λi .(6.8) This implies that if we define Bi=Wi/(Wi+ 2ηi), then: Bi|λi, ηi∼Beta k 2, λi. Now, we have that the conditional distribution of Gi/xiis Gi|λi, ηi, xi∼Ga(λi, ηixi)so that, conditional on Xi=xi,Wihas an inverse gamma distribution, that is Wi|λi, ηi, xi∼Inv Gamma(λ, ηixi).(6.9) Recalling that Zi∼N(0,Ik), we can calculate the posterior conditional distribution of Zivia Bayes theorem by combining this with (6.9, 6.8) to give: f(zi|wi, λi, ηi) = (wi+ 2ηi)2N(zi|0,Ik)Inv Gamma (wi|λi, ηixi) 2ηiBeta bi|k 2, λi. where bi=wi wi+2ηi,xi=|zi|2and N(·|·,·),Inv Gamma(·|·,·)and Beta(·|·,·)represent the normal, inverse gamma and beta density functions respectively. Finally, we have that the density of Yiis: f(yi|µ,Σ, λi, ηi, wi) = w k 2 i|Σ|−1 2f(zi=√wiA−1(yi−µ)) As previously, given sample data y1,y2, ..., yn, then we can derive the relevant conditional posterior distributions needed to set up a Gibbs sampler as follows. Firstly, the conditional posterior distribution of µis f(µ|w,Σ,λ,η,y)∝ n Y i=1 f(yi|µ,Σ, λi, ηi, wi)N(µ|m,V). Values of µcan be sampled using, for example a random walk sampler centred on the current value of µ. 106 Chapter 6. A Bayesian semi-parametric approach to normal/independent and elliptical distributions In a similar way, we have: f(σ2|µ,R,w,λ,η,y)∝ n Y i=1 f(yi|µ,Σ, λi, ηi, wi)f(σ2) f(R,D|µ,σ2,w,λ,η,y)∝ n Y i=1 f(yi|µ,Σ, λi, ηi, wi)f(R,D) where in these formulae, we recall that Σ = diag(σ)Rdiag(σ). In a similar way to the previous section, we can sample values of σ2using a Metropolis Hastings pass via an adaptive logit sampler and we can sample from the distribution of R,Dfollowing the approach of Fang and Li [1999]. The conditional posterior distribution of Wihas a much simpler form which can be sampled directly. We have: Wi|µ,Σ,λ,η,yi∼Ga k 2,gi 2 where giis calculated as in (6.4). Finally, the distributions for the Dirichlet process parameters have the same form as for the NI model, but replacing the values wiwith gifor i= 1, ..., n, throughout and can thus be sampled as previously. 6.4 Numerical illustrations In this section we show how the proposed methodology for estimating the NI and elliptical distribution performs on several simulated and real datasets. With respect to synthetic data, we simulate data from both the NI and elliptical (but non-NI) models and apply the fitting algorithms described in Section 6.3. 6.4.1 Simulated data from a NI distribution In this example, we consider a sample of size n= 500 of a bi-dimensional NI distribution (Student-t with 4 degrees of freedom) with parameters µ= (10,1) and Σ = 0.6 0.4 0.4 0.4!. The algorithms were run a total of 250000 iterations with a burnin period of 10000 iterations. Figure 6.1 depicts the convergence of the algorithms by showing the evolution of the average parameters µ1and Σ11. Figure 6.2 shows the fits to the histograms of both components provided by the algorithm to fit an elliptical distribution (solid line) versus that to adjust a NI distribution (dashed line). And, Figure 6.3 shows the contours plots of the fitted elliptical (left panel) and NI (right panel) distributions. Finally, Table 6.1 depicts the starting points for initializing the model parameters, 95% credible intervals when an elliptical and NI distributions are used to fit the data, according to the algorithms proposed in Section 6.3. 6.4. Numerical illustrations 107 Figure 6.1: Averaged µ1(top) and Σ11 (bottom) from the total number of runs of MCMC for the NI simulated dataset. Left panels: NI fit. Right panel: Elliptical fit. 6.4.2 Simulated data from an elliptical distribution In this section we consider an elliptical distribution with G∼Weibull(5,2). The sample size n, and the parameters µand Σare defined as in Section 6.4.1. Figure 6.4 shows the scatterplot of the data generated from the non-NI elliptical model. It can be observed that the elliptical dataset has a doughnut form. Figure 6.5 reports the plots of the two components of Z(left and right panel, respectively) versus the values of W. We can visualize how both variables, Z and W, seem to be quadratically dependent. Thus, these data follow an elliptical but non-NI distribution. 108 Chapter 6. A Bayesian semi-parametric approach to normal/independent and elliptical distributions Figure 6.2: Fit to the histogram of the simulated NI data (t4) by an elliptical distribution (solid line) and a NI distribution (dashed line). Left panel: first component. Right panel: second component. Figure 6.3: Contour plots from fitted joint distributions for the simulated NI data (t4). Left panel: Elliptical fit. Right panel: NI fit. Having run a total of 250000 iterations and 10000 iterations for the burning period and having reached the convergence of the algorithms, Figure 6.6 shows the fits to the empirical distribution function of both components (black line) by the elliptical distribution (green line) as well as the NI distribution (red line). It is clear that these data cannot be modeled well via a NI distribution. Chapter 7 General conclusions and future work 116 118 Chapter 7. General conclusions and future work In this PhD dissertation, the disciplines of Operations Research and Statistics constitute the common thread to propose novel solutions in real-life situations regarding with complex datasets. The challenge of learning and interpreting information from complex data is at the core of the current Statistical Science. New computational methods to deal with current datasets, by means of the powerful tools Mathematical Optimization and Bayesian Data Analysis, have been presented through the chapters forming this dissertation. In recent years, datasets have been usually characterized by a large number of features. This fact may have negative consequences in terms of the comprehensibility of solutions and, thus, the search for more interpretable solutions has recently led to the development of sparse multivariate techniques. Specifically, chapters 2 and 3 demonstrate different tradeoffs between sparsity and predictive performance of linear regression models, whereas Chapter 4 seeks sparse Bayesian classification models. Besides, when dealing with real-world applications where there exist groups at risk, it is more important to improve performance rates for such groups. To this aim, in this thesis, we apply mathematical optimization tools, which seem to be suitable. Then, in chapters 2 and 5 we propose two new constrained optimization models: one overall criterion is optimized, while constraints on the efficiency measures under consideration for the individuals of interest are introduced. Additionally, when dealing with datasets with extreme instances (heavy-tailed data), Bayesian statistical techniques have proven useful, as Chapter 6 also confirms. The works presented in this PhD dissertation can be extended. Next, we discuss some open problems. •In Chapter 2 a novel version of the Lasso in which quadratic performance constraints are added to Lasso-based objective functions is introduced. A possible extension for this work could be to change the objective function. For example, for the sake of dealing with strongly correlated predictors, it may be of interest to change the objective function by that of the elastic net. Another non-straightforward extension could be to address classification problems (via the logistic regression) instead of regression problems. In this case, we would not address a quadratic formulation. •A new methodology to deal with hierarchical categorical variables in linear regression models has been studied in Chapter 3. A number of extensions to this work are the following. Firstly, when the number of categories is large, instead of solving Problem (3.11) considering all the categories at once, a sequential pruning can be used instead. The main idea is to consider subtrees in Tjand try to compress their categories solving Problem (3.11) sequentially. Another option to deal with large number of categories is to cluster them based on a dissimilarity, see Carrizosa et al. [2017a]; Cerda et al. [2018] and references therein. Secondly, our methodology can be extended to generalized linear models [Tibshirani, 1996], where, instead of predicting the response variable as in (3.1), a non- 119 linear relationship between the response variable and the predictors is through a linkage function. However, this extension makes the optimization problem highly nonlinear and its resolution is very challenging. •In Chapter 4 a novel feature selection methodology characterized by three main features is presented. Namely, (1) sparsity is achieved taking into account the correlation among the features, (2) different performance measures can be used to guide the selection of features and (3) performance constraints on groups of interest can be included. It has been explored in the case of the NB classifier due to its tractability and good performance, but other classifiers could have also been tested instead. •Chapter 5 is devoted to introduce a constrained version of the classic NB classifier with the aim of improving classification rates in the different classes, avoiding the use of misclassification costs. Although distributional assumptions have been considered throughout this chapter, a possible extension to this work is to consider non parametric estimation for the density function for continuous attributes via kernel density estimation. •Finally, elliptical distributions are examined in Chapter 6. In particular, semi-parametric inference has been carried out for both elliptical and NI distributions, using Dirichlet process mixture models. To continue this work we have to undertake model selection. For model selection purposes, one possibility would be to try to calculate the marginal likelihoods for the various models considered and then to use the Bayes factor to select an optimal model, following Basu and Chib [2003]. However, other simpler approaches which are appropriate for mixture models can also be explored [Spiegelhalter et al., 2002]. Additionally, a more extensive computation study with applications in Finance, where heavy-tailed distributions are commonly used, is planned to be carried out. List of Figures 2.1 Heat map of ˆ βCSCLasso 1(λ)using prostate dataset.............. 25 2.2 Path of solutions under Lasso (top) and CSCLasso (bottom) for prostate dataset ...................................... 28 2.3 Path of solutions under Lasso (top) and CSCLasso (bottom) for prostate dataset when λincreases ............................. 29 2.4 Median overall MSE over the testing sets and NZ percentage under the choice p= 100 ...................................... 35 2.5 Median MSEkover the testing sets for k= 7,...,20 under p= 100 features and, nk= 150 (top) and nk= 300 (bottom). Each subgraph represents one group and the Y-axis shows the different percentages of improvement . . . . . 36 2.6 A two-dimensional graph of the logarithm of the user times in seconds for nk= 300 as pincreases.............................. 37 3.1 Tree representation of the variable geography in the cancer-reg dataset . . 49 3.2 Pareto frontier for MSE versus the number of coefficients to be estimated in the reduced model for the hierarchical categorical variable geography in the cancer-reg dataset .............................. 50 3.3 Less granular representations for the geography variable in the cancer-reg dataset....................................... 51 3.4 Tree associated with the variable CRIM in the housing dataset after being discretized..................................... 52 3.5 Pareto frontier for MSE versus the number of coefficients to be estimated in the reduced model for the hierarchical categorical variables in the housing dataset 53 3.6 Less granular representations for the first six hierarchical categorical variables in the housing dataset for the solution in Figure 3.5 with MSE∗((S∗ j)j∈J ) = 23.37 and c= 32. Note that this is the solution that achieves the minimum AIC 53 3.7 Less granular representations for the last six hierarchical categorical variables in the housing dataset for the solution in Figure 3.5 with MSE∗((S∗ j)j∈J ) = 23.37 and c= 32. Note that this is the solution that achieves the minimum AIC 54 120 List of Figures 121 3.8 Trees associated with the two hierarchical categorical variables in the synthetic dataset together with βS jl,l∈ L(Tj)....................... 55 3.9 Pruned tree and less granular representation of the two hierarchical categorical variables in Figure 3.8 from the synthetic dataset . . . . . . . . . . . . . . . . 56 3.10 Pareto frontier for MSE versus the number of coefficients to be estimated in the reduced model for the synthetic dataset for different σ2values . . . . . . . 59 3.11 Average MSE (10-fold CV) versus the imposed threshold cwhen σ2changes . 60 4.1 Heatmap associated to matrix M(based on MI correlation) corresponding to the australian dataset ............................ 69 4.2 Cluster dendrogram (based on MI correlation) corresponding to the australian dataset....................................... 70 4.3 Scalability..................................... 75 4.4 Average accuracy, sparsity and CPU time (10 runs 10-fold CV) for breast cancer,wine,mushroom,waveform,ISOLET and Multiple Features datasets ...................................... 77 4.5 Average performance, sparsity and CPU time (10 runs 10-fold CV) for SPECTF, german and page blocks datasets...................... 79 5.1 Boxplot of the importance of the variables for RCV1 dataset. . . . . . . . . . . 90 5.2 Scalability: X-axis represents the number of instances whereas each line the numberoffeatures. ................................ 94 5.3 Scalability: X-axis represents the number of features whereas each line the numberofinstances. ............................... 94 6.1 Averaged µ1(top) and Σ11 (bottom) from the total number of runs of MCMC for the NI simulated dataset. Left panels: NI fit. Right panel: Elliptical fit. . . . 107 6.2 Fit to the histogram of the simulated NI data (t4) by an elliptical distribution (solid line) and a NI distribution (dashed line). Left panel: first component. Right panel:secondcomponent. .........................108 6.3 Contour plots from fitted joint distributions for the simulated NI data (t4). Left panel: Elliptical fit. Right panel:NIfit. .....................108 6.4 Scatterplot of the non-NI, elliptical data . . . . . . . . . . . . . . . . . . . . . 109 6.5 The X-axis shows the Zvalues (Left panel: First component. Right panel: Second component) whereas the Y-axis the variable Wfor the simulated non- NIellipticaldataset................................110 6.6 Fit to the empirical distribution function of the simulated non-NI elliptical dataset (black line) by the elliptical distribution (green line) and the NI distribution (red line). Left panel: first component. Right panel: second component.111 122 List of Figures 6.7 Evolution of µ1(top) and Σ11 (bottom) from the total number of runs of MCMC for the real dataset. Left panel: NI fit. Right panel: Elliptical fit. . . . . . . . . 113 6.8 Fit to the histogram of the real dataset by an elliptical distribution (solid line), an NI distribution (dashed line) and a multivariate normal (dotted line). Left panel: first component. Right panel: second component. . . . . . . . . . . . . 114 B.1 Heat maps of ˆ βCSCLasso(λ) = ( ˆ βCSCLasso 1(λ),...,ˆ βCSCLasso 8(λ)) using prostate dataset ......................................136 B.2 Median MSEkover the test sets for k= 7,...,20 under p= 20 features and the two nkoptions. Each subgraph represents one group and the Y-axis shows the different percentages of improvement . . . . . . . . . . . . . . . . . . . . 137 B.3 Median MSEkover the test sets for k= 7,...,20 under p= 500 features and the two nkoptions. Each subgraph represents one group and the Y-axis shows the different percentages of improvement . . . . . . . . . . . . . . . . . . . . 138 B.4 Median overall MSE over the test sets and NZ percentage under the choice p= 20 (top) and p= 500 (bottom)........................139 B.5 Four perspectives of the logarithm of the user times in seconds for Lasso (bottom surface in the four graphics) and CSCLasso (top surfaces) models across a grid in nkand p..................................140 C.1 Variable selection process for breast cancer,waveform and wine, respectively, under the choices of S= 10 and q= 0.6...............146 C.2 Heatmaps .....................................151 List of Tables 1.1 Results obtained using prostate dataset.................... 6 1.2 Cost matrix for binary classification . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Results obtained using prostate dataset.................... 18 2.2 Median errors over testing sets for synthetic datasets . . . . . . . . . . . . . . 33 2.3 Median performance measures over testing sets for synthetic datasets . . . . . . 37 2.4 Median errors over testing set for gene expression dataset. Constraints imposed over Group 1 ................................... 38 2.5 Median errors over testing set for communities and crime dataset. Constraints imposed over Group 1 .............................. 39 3.1 Coefficients associated with four different representations for geography variable in the cancer-reg dataset ......................... 57 3.2 The predictor and the response variables in the housing dataset . . . . . . . . 58 4.1 Performance rate for all possible combinations of features in a multivariate normalsimulatedexample............................. 67 4.2 Average accuracy and sparsity (10 runs 10-fold CV) for simulated datasets. . . 74 4.3 Datasetsdescription................................ 75 4.4 Average performance and sparsity (10 runs 10-fold CV) for australian dataset using the sparse NB with different performance measures to select the setofvariables .................................. 80 5.1 Datasetsdescription................................ 89 5.2 Average Recall of classic NB (10 Monte-Carlo cross validation) . . . . . . . . 91 5.3 Testedthresholds ................................. 91 5.4 Average Recall values of SCNB (10 Monte-Carlo cross validation) for breast cancer ...................................... 92 5.5 Average Recall values of SCNB (10 Monte-Carlo cross validation) for SPECTF 92 5.6 Average Recall values of SCNB (10 Monte-Carlo cross validation) for page-blocks 92 5.7 Average Recall values of SCNB (10 Monte-Carlo cross validation) for abalone 92 123 124 List of Tables 5.8 Average Recall values of SCNB (10 Monte-Carlo cross validation) for yeast . 92 5.9 Average Recall values of SCNB (10 Monte-Carlo cross validation) for Satimage 93 5.10 Average Recall values of SCNB (10 Monte-Carlo cross validation) for RCV1 using 392 variables of the total . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.11 Average Recall values of SCNB (10 Monte-Carlo cross validation) for letter 93 6.1 Models comparison for the simulated NI (t4) dataset . . . . . . . . . . . . . . 109 6.2 Models comparison for the simulated elliptical (but non-NI) dataset . . . . . . 110 6.3 Models comparison for the real Mountain Top Radar dataset . . . . . . . 111 C.1 Accuracy of the sparse NB (10-fold CV) for breast cancer,waveform and wine datasets ................................145 C.2 Sparsity results (10-fold CV) for breast cancer,waveform and wine datasets ......................................147 131 min λ,qλ s.t. 2 nX0y=A0q, kqk∞≤λ. The constraint kqk∞is equivalent to |qs| ≤ λ,s= 0,1, . . . , p and the result follows. Proof of Proposition 3 The proof follows very closely that of Theorem 1. First of all, it shall be proven that h:β7→ E[(Y−X0β)2]is coercive. It is strictly convex on βsince its Hessian matrix, 2E[XX0], is positive definite due to Xis an absolutely continuous p-dimensional random variable: u0E[XX0]u=E[u0XX0u] = E[(X0u)2]>0, since P(X0u= 0) = 0. Moreover, λkAβk1is a convex function on βand, therefore, E[(Y− X0β)2+λkAβk1]is also a strictly convex function on β. On the one hand, the eigenvalues of the Hessian matrix are all positive and, in particular, the smallest eigenvalue, say γr, will be non-zero. On the other hand, using the spectral decomposition of a symmetric matrix, E[(Y−X0β)2] = β0E[XX0]β−2E[YX]β+E[Y2] = β0Q0DQβ−2E[YX]β+E[Y2]≥ ≥β0Q0DQβ− | 2E[YX]β|+E[Y2]≥γrkQβk2−kE[YX]kkβk+E[Y2] = =γrkβk2−kE[YX]kkβk+E[Y2], where, in the second-to-last step, the Cauchy-Schwarz inequality was used. As kβk → +∞, then E[(Y−X0β)2]→+∞, that is, the quadratic function h(β) = E[(Y−X0β)2]is coercive. The next step in the proof is to transform the original true problem (2.17) into an equivalent one with a feasible compact region B∗. Given β∗∈B, since h(β) = E[(Y−X0β)2]is coercive, there exists Rsuch that E[(Y−X0β)2]> E[(Y−X0β∗)2+λkAβ∗k1], for all βwith kβk> R. For that reason, the problem (2.17) can be reduced to the feasible compact region B∗=B∩{β:kβk ≤ R}, which implies that the optimal solution is reached. Finally, the uniqueness of solution is a consequence of the strict convexity of the objective 132 Appendix A. Proofs function. Proof of Theorem 2 For the sake of simplicity, βCSCLasso(λ)and ˆ βCSCLasso(λ)will be denoted henceforth by β and ˆ β, respectively. In addition, let us consider the nonempty compact set C=B∩ {β: kβk ≤ R}, where Ris chosen according to the proof of Theorem 2.3.1. Theorem 2 is a direct consequence of Theorem 5.3 in Shapiro et al. [2009] under some technical conditions, namely: C1. The expected value function E[(Y−X0β)2+λkAβk1]is finite valued and continuous on C. C2. 1 nPn i=1((yi−x0 iβ)2+λkAβk1)converges to E[(Y−X0β)2+λkAβk1]w.p. 1, as n→ ∞, uniformly in β∈C. Let us denote F(β,(Y, X)) = (Y−X0β)2+λkAβk1. Then, the previous conditions C1 and C2 are consequences of Theorem 7.48 in Shapiro et al. [2009] provided that A1. for any β∈C, the function F(·,(Y, X)) is continuous at βfor almost every (Y, X), A2. the function F(β,(Y, X)), with β∈C, is dominated by an integrable function, A3. the sample is i.i.d. Given (Y, X), the function (Y−X0β)2+λkAβk1is continuous at βfor any β∈C, and therefore A1 is fulfilled. The sample is i.i.d. by hypothesis, and thus A3 holds too. Finally, in order to prove A2, it is necessary to find a measurable function g(Y, X)>0such that E[g(Y, X)] <∞and, for every β∈C,|F(β,(Y, X)) |≤ g(Y, X)w.p. 1. Using the Cauchy- Schwarz inequality, one has, |F(β,(Y, X)) |=|(Y−Xβ)2+λkAβk1|= =|Y2−2YX0β+β0XX0β+λkAβk1|≤ ≤Y2+ (X0β)2+ 2 |YX0β|+λkAβk1= =Y2+|X0β|2+2 |YX0β|+λkAβk1≤ ≤Y2+kXk2kβk2+ 2kYXkkβk+λkAβk1. Let M1and M2be given by M1= max β∈CkβkM2= max β∈C|Aβ| 133 which are well defined due to the compactness of C. Therefore, gcan be chosen as g(Y, X) = Y2+M2 1kXk2+ 2M1kYXk+λM2, which is positive and, since E(kXk2)<∞,E(Y2)<∞,E(kYXk)<∞, its expected value is finite. In consequence, A2 holds and the proof is concluded. Appendix B Further results 134 136 Appendix B. Further results Figure B.1: Heat maps of ˆ βCSCLasso(λ)=(ˆ βCSCLasso 1(λ),...,ˆ βCSCLasso 8(λ)) using prostate dataset 137 Figure B.2: Median MSEkover the test sets for k= 7,...,20 under p= 20 features and the two nkoptions. Each subgraph represents one group and the Y-axis shows the different percentages of improvement 138 Appendix B. Further results Figure B.3: Median MSEkover the test sets for k= 7,...,20 under p= 500 features and the two nkoptions. Each subgraph represents one group and the Y-axis shows the different percentages of improvement 139 Figure B.4: Median overall MSE over the test sets and NZ percentage under the choice p= 20 (top) and p= 500 (bottom) 140 Appendix B. Further results To fully understand how the computation time behaves depending on nkand pvalues, a grid in both parameters have been inspected. Figure B.5 displays the logarithm of the user times in seconds obtained under Lasso and CSCLasso models when nkand pchange. The perspective drawn in the top left figure shows that Lasso model (bottom surface) is solved faster and in a smoother way. Besides, whereas smaller times are obtained for both methods when nkand p are small, the biggest times are associated to nk= 300 and p= 500. Figure B.5: Four perspectives of the logarithm of the user times in seconds for Lasso (bottom surface in the four graphics) and CSCLasso (top surfaces) models across a grid in nkand p