Renée El

Melhem Bât

Pascal Blaise

M Stéphane

ÉLECTROTECHNIQUE E E A Électronique

M Philippe Delachartre
email: philippe.delachartre@insa-lyon.fr

Bénédicte Lanza

Mme Sandrine
email: sandrine.charles@univ-lyon1.fr

Mme Sylvie

M Hamamache Kheddouci

M Stéphane Benayoun
email: stephane.benayoun@ec-lyon.fr

Stéphanie Cauvin

M Jocelyn Bonjour

M Christian Montes

Dr Jean- Marc Petit

Dr Pierre-Edouard Portier

Sylvie Dumas

Worldwide, highway accidents have important social and financial impacts. To reduce their frequency and gravity, crash prediction models (CPM) are used to identify hazardous roadway segments and to provide actionable clues about the associated risk factors. CPM are either parametric statistical models, in particular generalized linear models (GLM), or machine learning models with a large number of parameters without associated uncertainty estimates (e.g., ensemble of decision trees, supportvector machine . . . ). Simple parametric models tend to be more interpretable but less effective than highly flexible non-parametric models that work like black-boxes. When pondering high stake decisions, such as in the context of highway safety, field experts expect predictive models to be both effective and glass-box interpretable. The models must assist them in conceiving and deploying preventive or remedial safety actions.

As such, we contribute to enhancing the predictive performance of parametric models while maintaining their interpretability. In the first place, a well-chosen hierarchical structure can handle correlations among groups of observations and significantly improve the quality of the models' predictions and of their interpretation. We propose to learn it by leveraging the output of a post-hoc explainability framework (viz., SHAP) applied to a highly flexible black-box model (viz., XGBoost). In our first contribution, this hierarchical structure informs a Bayesian multilevel GLM. Moreover, in an effort to further improve the predictive performance of the model without deteriorating its interpretability, we propose to extend its linear functional form to account for major first-order interactions between explanatory variables. These interactions are learnt from the data by analyzing the results of a trained self-organized polynomial network.

In our second contribution, we exploit the hierarchical structure even better by replacing the GLM with a simulated annealing based multi-objective symbolic regression algorithm to automate feature engineering and feature selection. Thus, by computing a cluster-specific ranking of expansions of regularized linear models ordered by increasing complexity, we facilitate a dynamic interpretative process which makes it possible to discover effective, efficient and interpretable predictive models.

Experiments have been conducted on a highway safety dataset and on more than ten public datasets covering classification and regression tasks. They show promising results with our two contributions outperforming traditional glass-box interpretable models while getting close to the best non-parametric models. Finally, we illustrate the benefits of our approach by introducing, on a realistic case study, an application we designed for highway safety experts. iv

Résumé

Dans le monde entier, les accidents de la route ont des impacts sociaux et financiers importants. Pour réduire leur fréquence et leur gravité, les modèles de prédiction d'accidents (CPM) sont utilisés pour identifier les segments de route dangereux et fournir des indices exploitables sur les facteurs de risque associés. Les CPM sont soit des modèles statistiques paramétriques, en particulier des modèles linéaires généralisés (GLM), soit des modèles d'apprentissage automatique avec un nombre important de paramètres sans estimation d'incertitude associée (e.g., ensemble d'arbres de décision, machine à vecteurs de support . . . ). Les modèles paramétriques simples ont tendance à être plus interprétables mais moins performants que les modèles non paramétriques très flexibles qui fonctionnent comme des boîtes noires. Lorsqu'ils réfléchissent à des décisions à fort enjeu, comme dans le contexte de la sécurité routière, les experts métier s'attendent à ce que les modèles prédictifs soient à la fois performants et interprétables. Les modèles doivent les aider à concevoir et à déployer des actions de sécurité préventives ou correctives.

Dans ces travaux, nous contribuons à améliorer les performances prédictives des modèles paramétriques tout en conservant leur interprétabilité. En premier lieu, une structure hiérarchique bien choisie peut gérer les corrélations entre groupes d' observations et améliorer significativement la qualité des prédictions des modèles et leur interprétation. Nous proposons de l'apprendre en exploitant le résultat d'une méthode d'interprétabilité post-hoc (viz., SHAP) appliquée à un modèle boîte noire flexible (viz., XGBoost). Dans notre première contribution, cette structure hiérarchique informe un GLM bayésien multiniveaux. De plus, dans le but d'améliorer encore les performances prédictives du modèle sans détériorer son interprétabilité, nous proposons d'étendre sa forme fonctionnelle linéaire pour tenir compte des interactions majeures de premier ordre entre variables explicatives. Ces interactions sont apprises à partir des données en analysant les résultats d'un réseau polynomial auto-organisé.

Dans notre deuxième contribution, nous exploitons encore mieux la structure hiérarchique en remplaçant le GLM par un algorithme de régression symbolique multi-objectif basé sur le recuit simulé pour automatiser la sélection des variables explicatives et l'extraction de caractéristiques (viz., interactions, transformations de v variables explicatives). Ainsi, en calculant un classement spécifique à chaque cluster des expansions de modèles linéaires régularisés ordonnés par complexité croissante, nous facilitons un processus d'interprétation dynamique qui permet de découvrir des modèles prédictifs efficaces, efficients et interprétables.

Des expériences ont été menées sur un jeu de données de sécurité routière et sur plus de dix jeux de données publics couvrant des problèmes de classification et de régression variés. Les résultats obtenus sont prometteurs étant donné que nos deux contributions surpassent les modèles interprétables traditionnels et se rapprochent des meilleurs modèles non paramétriques boîtes noires. Enfin, nous illustrons les bénéfices de notre approche en présentant, sur une étude réelle de cas, une application que nous avons conçue pour les experts de la sécurité routière. vi life to its fullest, will remain engraved in me. I feel so lucky to have known you.

The best for the end, I would like to thank my love, Manon. You can't imagine how much I admire you for your passion, your perseverance, for being you. Since you came into my life, everything seems happier and simpler. With you by my side, I feel able to reach my most unattainable goals and live my wildest dreams. Je t'aime, du plus fort que je peux.

Context and objectives

Road safety is a socio-economic concern: according to the World Health Organization [START_REF]Road traffic injuries[END_REF], approximately 1, 35 million people are killed each year on roadways around the world. Related expenses average 3% of the gross domestic product of a country. As stated by the French Road Safety Observatory [START_REF]Road safety annual report[END_REF], these costs grow exponentially with the severity of the accidents. In 2019, a property damage only accident incurred expenses up to 5 258 euros while the average cost of an accident with at least one fatality was 3 429 000 euros 1 . Since the 1970s, road safety has become a major challenge for successive French governments. Many new safety policies have been engaged in, such as reductions of speed limits, the multiplication of fixed radars, the increase in the amount of speeding fines. At the same time, vehicles became more secure. Road safety was therefore continuously improved. However, since 2010, the positive trend was attenuated and the reduction in the fatality rate started to stagnate. For this reason, the French government, supported by the European Union, motivated new safety policies whose main issuers would be local authorities and road managers. Among them, APRR 2 group, subsidiary of Eiffage, finances, maintains and manages the infrastructure of a 2 323 km-long highway network (see Fig. 1.1) in return for toll collection. The objective of the company is to offer high-performance transport infrastructure and support road users at every stage of their travel, guaranteeing them the best conditions for traffic flow and safety.

To secure their network, field experts of APRR leveraged an extension of the User's 1 These amounts include the various costs incurred for territorial and local authorities (hospitalization, insurance, etc.) and for road managers (marking, road works, etc.)

2 Autoroutes Paris-Rhin-Rhône: https://aprr.com/en Safety on Existing Roads (SURE) approach, a comprehensive and global framework for network safety management initiated by the French government. First, they implemented a multi-criteria decision analysis method to compute a ranking of segments based on the values of different crash-related criteria (e.g., annual average daily traffic, fatality or injury rates, . . . ). Then, after studying crash reports and conducting on-site visits, field experts elaborated pluriannual action plans (operation, maintenance, . . . ) on segments where the safety improvement were potentially the greatest.

In the end, this method contributed to a significant reduction in fatality and severity rates. Nevertheless, SURE does not work on reducing the crash likelihood and does not explain the mechanisms affecting crash-prone areas. Therefore, decision makers have opted for the development of new intelligent systems inspired by the latest advances in Artificial Intelligence (AI). A collaboration has been initiated between LIRIS3 laboratory and Data New Road4 , a subsidiary of APRR dedicated to the harnessing of the company's historical data. The main objective of this collaboration was to build an AI system to estimate accurately future hazardous segments and to identify risk factors. A system that meets these requirements will be a fundamental element in the implementation of new proactive safety policies, and will therefore have a beneficial and directly measurable impact.

Crash prediction models

In the field of road safety, research has focused on obtaining a better understanding of the mechanisms that affect the probability of an accident. Given the absence of precise data related to human behavior (lack of attention, fatigue, alcohol, . . . ), crash prediction models (CPM) do not provide cause-and-effect relationships that would explain the occurrence of an accident. Instead, they emphasize the use of factors that could influence the likelihood of an accident over a given period (e.g., year, month) and in a given geographical area. CPM are trained on historical data to discover hidden patterns between the crash-related variable and several explanatory variables (e.g., traffic, speed limit, altitude, . . . ). They are mainly used to identify risk factors in order to steer the evolution of safety policies.

In their survey [START_REF] Lord | The statistical analysis of crashfrequency data: A review and assessment of methodological alternatives[END_REF], Lord and Mannering provide a broad perspective on the variety of data-related issues raised by crash count prediction: over-dispersion of count data, temporal and spatial correlations due to multiple measurements of a same location at different times, fixed parameters that cannot adapt from one roadway to the next, low sample-mean due to the sparsity of crashes, non-linear relationships between crash-frequencies and explanatory variables, etc. Most of theses issues are made more prominent with the use of parametric models, in particular generalized linear models (GLM). Indeed, GLM must undergo many transformations to adapt to the crash count prediction context (e.g., the choice of a non-normal likelihood distribution, the integration of random effects and hierarchical models, etc). Whereas non-parametric approaches (e.g., neural networks, tree-based algorithms, SVM...) will deal with most of these issues without the need for specific adaptations and will usually offer better predictive performances than parametric models. On the other hand, this improvement in predictive power comes at a cost. These models, which involve a large number of parameters, are often considered as black boxes that do not allow direct understanding of the mechanisms that led to a prediction. Thus, when designing CPM, both dimensions (viz., performance and interpretability) must be considered. However, while performance is easily measured with widely adopted metrics, quantifying interpretability remains more debated.

Model interpretability: concepts, taxonomies

In our study, we consider a predictive model to be glass-box interpretable when it quantifies explicitly (i.e., not by simulation) the marginal effects of the explanatory variables and also possibly of a few simple transformations of these variables (e.g., multiplicative interactions, log transforms, etc.). Thus, we stress the importance of favoring simple functional forms namely, expansions of linear models.

On the contrary, non-parametric models allow for very flexible, but often complex, functional forms. For these models, some sense of the marginal effects can be gained by studying ceteris paribus profiles. They represent the influence of an explanatory variable by assuming that all the other variables remain constant. In the case of flexible functional forms, the choice of the fixed values for all the variables but the one for which the marginal effects are to be observed has a great influence on the interpretation of the effect. Often, different choices (if not all) should be considered to provide a more faithful view of the marginal effects. Of course, the number of configurations of fixed values grows quickly with the number of explanatory variables. Some heuristics have been used to manage this complexity. For example, Beck et al. [START_REF] Beck | Improving quantitative studies of international conflict: A conjecture[END_REF] chose to "hold constant the other variables at two values: high and low probability" of the target.

To truly measure the marginal effects, one should, as stated in [START_REF] Hainmueller | Kernel regularized least squares: Reducing misspecification bias with a flexible and interpretable machine learning approach[END_REF], "estimate partial derivative with respect to each input variable and at each observation through repeated simulations". This strategy is computationally prohibitive when implemented by simulating data from a complex model. Among the oldest computationally convenient variants, we find the permutation based methods applied to variable importance, partial dependence plots or individual conditional expectation plots. However, Hooker and Mentch [START_REF] Hooker | Please stop permuting features: An explanation and alternatives[END_REF] have found that "when features in the training set exhibit statistical dependence, permute-and-predict methods can be highly misleading when applied to the original model". More refined approximation methods have been proposed, such as LIME or SHAP. They are also based on observing the effects of a perturbation of a given instance on the output of the model. These approaches are often referred to as model-agnostic post-hoc explanation methods. However, recent works, in particular [START_REF] Slack | Fooling lime and shap: Adversarial attacks on post hoc explanation methods[END_REF] tend to show that it can be difficult to assess their reliability and robustness. Also, Garreau and Luxburg [START_REF] Garreau | Explaining the explainer: A first theoretical analysis of lime[END_REF], in a theoretical study of LIME, verified by simulations, showed that, while this post-hoc approach "discovers interesting features, it might forget some important features and the surrogate model is not faithful".

Generalized Additive Models (GAM) [START_REF] Trevor | Generalized additive models[END_REF] can be located halfway between complex flexible functional forms and simple glass-box models. They give some access to the marginal effects by enabling plots of the expectation of the target given the values taken by each explanatory variable. Lou et al. [START_REF] Lou | Intelligible models for classification and regression[END_REF] designed a GAM variant based on boosted trees, the Explainable Boosting Machine (EBM). It is often almost as accurate as flexible black-box models. However, Chang et al. [START_REF] Chang | How interpretable and trustworthy are gams[END_REF] observed that different GAM algorithms, while offering comparable results in terms of accuracy, can lead to very different interpretations of the predictions.

Finally, we consider that simple linear models with the inclusion of some transforms of the original explanatory variables and the use of a regularization term (e.g., ridge or LASSO) to control the bias/variance trade-off remain promising among the intrinsically glass-box approaches. However, they come with key challenges. If one doesn't want the user to be solely responsible for deciding which functions of the original variables are to be considered, a strategy must be adopted to explore a combinatorial space of potential functional forms. Also, for the model to stay interpretable in presence of a large number of explanatory variables, the number of relevant coefficients shouldn't be too large. However, variable selection is not without drawbacks. In particular, it suffers from the risk of hiding the main effects behind more complex correlated terms.

Machine learning formalization for crash prediction

Data description

Throughout this work, we will refer to x as the vector of explanatory variables describing a roadway segment during a given period. Each explanatory variable can be either continuous, discrete or categorical. Continuous variables can assume an infinite number of real values within a given interval. As opposed to a continuous variable, a discrete variable can assume only a finite number of numerical values within a given interval, such as the number of interchanges on a roadway segment in our context. The categorical variable refers to a variable that takes values in a discrete unordered set of categories {c 1 , ...c n }. Sometimes, categorical variables have only two categories and are therefore called binary variables. For instance, in our crash prediction problem, the variable presence of tunnel indicates whether a tunnel is present or not on a roadway segment. The dependent variable, also known as the target variable, will be denoted y. In the context of crash count prediction, y is a discrete positive integer that quantifies the number of crashes observed on a roadway segment for a given time interval. Note that in the literature, the problem can be tackled with a binary categorical variable. In this case, the focus is not to predict a crash count, but to know whether or not a roadway segment will be accidental given the values of the explanatory variables. These models require the definition of a threshold that is used to differentiate crashprone configurations from normal ones.

Supervised learning

In a supervised context, the predictive model must discover how to associate an observation x to a label y, based on a set of known examples {x i , y i } N i=1 . To model these dependencies, the process goes through two phases:

• training phase: the model is trained to find patterns in the training data that map the vectors of explanatory variables x to the dependent variables y. To learn this mapping, the model goes through a stage called model training, where its parameters (e.g., weights, bias) are updated so the function it approximates captures the most relevant dependencies among the variables. Sometimes, models involve hyper-parameters i.e, parameters that are determined before training and remain fixed afterwards. Thus, an additional step called validation step can be performed in order to verify the correctness of their values and adjust them if necessary.

• testing phase: the trained model is evaluated on a set of unknown samples with appropriate performance metrics (e.g., mean squared error for regression tasks, f 1-score for classification tasks). By doing so, one obtains a less biased estimate of the predictive performance of the model which will allow, firstly, to validate the generalization capacity of the model and secondly, to compare the predictive performance of the selected model with those of others.

However, doing a single train-test split will give, most of the time, an overoptimistic estimate of the predictive performance of the model. One way to mitigate this is to use k-fold cross-validation. In k-fold cross-validation, the dataset is split into k subsets of data (known as folds). The model is trained on k -1 subsets, and then evaluated on the remaining subset that was not used for training. This process is repeated k times, with a different test subset each time. The average performance obtained on the k test subsets gives an estimate of the model's predictive performance. In practice, k-fold cross-validation can serve multiple purposes such as defining the best set of hyper-parameters for a machine learning model or comparing multiple predictive models and selecting the one that has the best generalization ability.

Application to the highway dataset

In our work for highway safety analysis, the objective is to predict efficiently the crash count observed on highway segments, given the value of spatial explanatory variables (e.g., number of lanes, speed limit, . . . ) and temporal explanatory variables (annual average daily traffic, percentage of heavy vehicles). To learn these relationships, we leverage supervised learning: first, we start by generating a dataset with samples representing the crash count and the values of explanatory variables observed on a segment for a particular year (see table 1.1). The process of generating the dataset will be described in more detail in section 4.1. Then, we train a predictive model on a subset of randomly selected samples from the overall dataset (say for example 80%) to learn a function that maps the explanatory variables to the crash count. To efficiently help field experts, this function will have to meet different challenges presented in the next section. Finally, we evaluate the performance of the resulting function on the remaining samples.

Challenges and research questions

If we assume the observed data are generated from an unknown data-generating process f such that y = f (x) + ε, with ε being a zero-centered Gaussian noise of variance σ 2 , our objective is to find a function f that is as close as possible to the true but unknown process which generated the data. Given the available training data {(x i , y i )} N i=1 , we estimate f so it minimizes an objective function such as the mean squared error i (y i -f (x i )) 2 .

Bias-variance trade-off

To measure the ability of f to generalize to points outside of the training data, we compute the expected prediction error (EP E) of the prediction y * = f (x * ) + ϵ, with x * being an unseen sample. As shown by Hastie et al. [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF], EPE can be decomposed as:

EP E = E[(y * -f (x * )) 2 ] = (f (x * ) -E[ f (x * )]) 2 + E[ε 2 ] + E[ f (x * ) 2 ] -E[ f (x * )] 2 = bias( f ) 2 + σ 2 + var( f ) (1.1)
We observe that EP E embeds three different sources of error: σ 2 is the irreducible error that cannot fundamentally be reduced by any model, the bias measures the average error of f (x * ), and the variance var gives a measure of the variation of the prediction f (x * ) from one training dataset to another. An ideal model minimizes both the bias and the variance. However, in a world with imperfect models and finite data, it is nearly impossible to do both simultaneously. On the one hand, complex models (e.g., deep neural networks) approximate precisely training data but are exposed to over-fitting. On the other hand, simpler models with high bias, such as linear models, fail to capture important regularities in the data and are therefore prone to underfitting. Thus, efforts should be made to reach optimum model complexity for which the expected prediction error is generally close to its lowest level (see Fig. 1.2). One way of doing this is to use regularization to find a good trade-off between bias and variance. Likewise, dimensionality reduction and feature selection also decrease the variance of complex models.

Finally, when designing predictive models in a high-stakes context, it is necessary to not focus only on mitigating the bias-variance trade-off, but also to consider model interpretability because the understanding of the mechanisms that lead to a prediction is as valuable as the prediction itself.

Performance and interpretability Broadly speaking, non-parametric complex models, which cover big hypothesis spaces, often obtain very good performance on many predictive tasks. However, as they involve a large amount of parameters, they're often considered as black-boxes [START_REF] Barredo Arrieta | Explainable artificial intelligence (xai): Concepts, taxonomies, opportunities and challenges toward responsible ai[END_REF]. To gain insights into their underlying behaviors, post hoc explanation tools must be used. Nonetheless, these tools have several drawbacks w.r.t. the trust and veracity of the explanations they deliver [START_REF] Rudin | Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead[END_REF][START_REF] Slack | Fooling lime and shap: Adversarial attacks on post hoc explanation methods[END_REF]. Indeed, the representation of the model can be inaccurate in parts of the feature space and may even be misleading in the under-represented parts [START_REF] Molnar | Interpretable Machine Learning[END_REF]. Moreover, even though these tools can reproduce accurately the predictions of the original model, they may use completely different features. A well-known example is the proprietary prediction model COMPAS, a recidivism-risk scoring model used widely throughout the U.S. criminal justice system for parole and bail decisions, which has been criticized as being dependent on race based on the interpretations of an explanation model [START_REF] Mattu | Machine bias[END_REF] whereas it has been demonstrated afterwards that COMPAS does not seem to "depend strongly on either criminal history or proxies for race" [START_REF] Rudin | The age of secrecy and unfairness in recidivism prediction[END_REF].

Moreover, when black-box models do not behave as expected (e.g., excellent accuracy during training, but very bad performance in production or with testing data), trying to gain insights into the model is required for troubleshooting (debugging). Post hoc explanations are commonly used. However, combining a black-box model and an explanation model (or tool) makes the process of troubleshooting very sensitive. As the explanation is not always correct, it can be difficult to tell whether or not the black-box model is wrong [START_REF] Rudin | Interpretable machine learning: Fundamental principles and 10 grand challenges[END_REF].

Furthermore, as stated by Rudin [START_REF] Rudin | Interpretable machine learning: Fundamental principles and 10 grand challenges[END_REF], "black-box models often predict the right answer for the wrong reason, leading to excellent performance in training but poor performance in practice". In psychology, this is known as the clever Hans phenomenon [START_REF] Pfungst | the horse of Mr. Von Osten.) a contribution to experimental animal and human psychology[END_REF], that comes from a horse (of the same name) claimed to perform arithmetic and other intellectual tasks, but actually was watching the reactions of his trainer. In the AI context, this phenomenon describes black-box models that can learn different relationships between variables than the ones of the true data generating process and, at the same time, provide very accurate predictions. Thus, efforts should be made to avoid this phenomenon in domains where high-stakes decisions can be made from the analysis of predictive models such as healthcare, road safety, or criminal justice.

On the other hand, researchers working in the field of interpretable AI observe that linear models can perform very well on tabular data, especially when the full data science process is considered (e.g., variable selection, data preparation, model selection, . . . ) [START_REF] Rudin | Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead[END_REF]. This could arise from the Rashomon Effect [START_REF] Breiman | Statistical modeling: The two cultures (with comments and a rejoinder by the author)[END_REF], which characterize problems where many accurate-but-different models exist to describe the same data [START_REF] Semenova | A study in rashomon curves and volumes: A new perspective on generalization and model simplicity in machine learning[END_REF]. In [START_REF] Semenova | A study in rashomon curves and volumes: A new perspective on generalization and model simplicity in machine learning[END_REF], the authors show empirically that when the Rashomon set (i.e. the set of almost-equally-accurate models) is large, "most machine learning methods tend to perform similarly, and also in these cases, interpretable or sparse (yet accurate) models exist".

Lastly, interpretable models provide a simple and intelligible picture of the relationship between the explanatory variables and the dependent variable. This not only ease the process of elaborating new actions, but also allow end users to be confident in what they plan.

Publications and outline

Publications

In this thesis, we try to meet these challenges by extracting, from data and without prior expert knowledge, the information necessary to build efficient and highly interpretable models. In our first proposal, we consider that data can often be partitioned so that refined predictive models can apply to different parts more efficiently and more meaningfully than a global model. We discover such a structure by clustering the instances based on the features' scores returned by the SHAP [START_REF] Scott | A unified approach to interpreting model predictions[END_REF] post-hoc analysis of a flexible black-box model. This clustering then informs a Bayesian multilevel -hierarchical -generalized linear model. Moreover, we propose to integrate nonlinearities by adding to its underlying functional form, the most important interactions between explanatory variables learnt by a self-organized polynomial network. This first contribution, described in section 3.2 and section 3.3, was published in: Crash prediction for a French highway network with an XAI-informed Bayesian hierarchical model, 2020 IEEE International Conference on Big Data (Big Data) [START_REF] Veran | Crash prediction for a french highway network with an xai-informed bayesian hierarchical model[END_REF] In our second proposal, we propose to exploit the hierarchical structure even better by computing a cluster-specific ranking of expansions of regularized linear models ordered by increasing complexity. We design a symbolic regression approach for the automatic discovery of a mapping from the original explanatory variables to a basis expansion which adds both sparsity and flexibility by including transforms of the variables. This exploration of the space of potential expansions of linear models is guided by simulated annealing. More precisely, given a predictive performance metric and a complexity metric, the meta-heuristic search conducts a multi-objective optimization to return a list of Pareto optimal models. Through a dynamic interpretative process, end users can navigate around these models and select the best model according to their needs, from the simplest one highlighting main effects to a more specific one depicting particular configurations related to few instances in the dataset. This contribution, described in section 3.4, was published in:

Interpretable hierarchical symbolic regression for safety-critical systems with an application to highway crash prediction. Engineering Applications of Artificial Intelligence, 117:105534, 2023. [START_REF] Veran | Interpretable hierarchical symbolic regression for safety-critical systems with an application to highway crash prediction[END_REF]

Outline

In chapter 1, we started by introducing the context and stakes of road safety analysis and then described the numerous challenges crash prediction models have to face. In the next chapter, we will first provide more details on the different concepts of model interpretability through a review of the state of the art before introducing how crash predictions and explanations have been tackled in the community. In chapter 3, we present our contributions. More precisely, after giving an overview of our methodology in section 3.1, we will describe in section 3.2 the first module of our methodology, which consists in finding a hierarchical structure in the dataset. In section 3.3, we will explain how we integrate this hierarchical structure to enhance both the interpretability and efficiency of GLM. We will also explain how non-linearities are discovered and afterwards embedded in the GLM to enhance, yet again, its predictive performance. In section 3.4, we will describe how the aforementioned methodology can be improved by replacing the GLM by a symbolic regression approach that allows us to capture more relevant non-linearities while providing sparse and highly interpretable models. Then, in chapter 4, we will present the experiments carried out on the highway network dataset and on thirteen public datasets covering different tasks (viz., regression and classification) which underline very promising results as our framework outperforms fully interpretable models (e.g., linear models, shallow decision trees) while getting close to non-parametric models. Finally, in section 4.6, we unveil its potential for guiding a dynamic interpretative process with a realistic case study on the highway network dataset which highlights the benefit of using multi-objective optimization to help field experts develop new safety policies.

Chapter 2

Related work

Model interpretability

Properties

AI systems are used to assist field experts in making high stake decisions which may indirectly affect humans' lives. Thus, understanding how does the system behave is of paramount importance. In many areas, knowing the predictions made by a model is not enough. There is an emerging need that these models, besides being efficient, compel to numerous properties regarding their explainability [START_REF] Barredo Arrieta | Explainable artificial intelligence (xai): Concepts, taxonomies, opportunities and challenges toward responsible ai[END_REF]:

• trustworthiness: the models will act as intended when facing a given problem • informativeness: to support decision making, they provide a great amount of information on the problem being tackled

• confidence: they communicate on the confidence of their working regime

• accessibility: they allow end users to get more involved in their development cycle

• interactivity: they allow straightforward interactions with end users.

The design of interpretable AI systems has been tackled with two distinctive approaches rallied around the field of explainable Artificial intelligence (XAI). On the one hand, there are models that are inherently interpretable as their complexity is restricted by their design (e.g., GLM, decision trees). On the other hand, complex models (e.g., deep Neural Network, ensemble of decision trees) act as black-boxes and do not provide, at least initially, any information on the relationship that connects the explanatory variables to the dependent variable. Knowledge is often obtained by using post hoc explainability techniques.

Intrinsic interpretability

Some models are inherently interpretable and indeed meet the above properties. They are often perceived as glass-box models as the relation between the explanatory variables and the dependent variable is explicit, thus simplifying the understanding of the marginal effects. Among them, we can find shallow decision trees, generalized linear models (GLM), generalized additive models (GAM).

Decision trees

Decision trees can be used to answer regression problems and classification problems. They are commonly found under the term Classification and Regression Trees (CART) introduced by Breiman [START_REF] Breiman | Classification and regression trees[END_REF]. These models aim at predicting the value of a dependent variable by learning simple decision rules from the features. They divide the input space by recursively splitting the data into subsets according to tests on features (see Fig. 2.1). The procedure to select the best split varies according to the problem at hand. For regression tasks, CART takes a feature and determines which cut-off minimizes the variance of the dependent variable. For classification tasks, CART seeks to minimize the Gini impurity of the class distribution defined as:

Gini = 1 - C i=1 (p i ) 2
where p i is the probability of an observation to be classified to one of the C classes. The algorithm continues the search-and-split in the new nodes until a stopping criterion is reached (e.g., minimum number of observations in a node before splitting, minimum number of observations in a terminal node). Finally, to predict the outcome of a new observation, CART averages all training instances in the terminal node reached by the new observation after browsing the tree.

Interpreting the prediction is straightforward as it combines, from the root node to the terminal node, simple decision rules. However, for deep trees, the high number of decision rules used for a prediction makes it hard to understand. Thus, in our study, we consider a CART model to be interpretable if its depth is, at most, equal to 5.

Generalized linear model

In a GLM [START_REF] Mccullagh | Generalized linear models[END_REF], the dependent variable Y is assumed to be generated from a particular distribution in an exponential family (e.g., normal, binomial, Poisson distributions). The mean µ of the distribution is connected to a linear predictor Xβ such that: where g is a link function introduced to map the non-linear transformed mean g(µ)

g(E[Y |X]) = g(µ) = Xβ (2.1)
to the definition space of the linear predictor, the real line. The choice of a relevant link function depends on the problem at hand. For instance, for a dependent variable generated from a binomial distribution, the preferable link function is the logit. For that model, Eq. 2.1 becomes:

logit(µ) = ln( µ 1 -µ ) = Xβ
which is also referred to as logistic regression. For a normally distributed dependent variable, the link function is simply the identity function id:

id(µ) = µ = Xβ
The coefficient β can be estimated by maximum likelihood estimation or Bayesian inference (see section 3.3.2 for a detail explanation of these techniques).

Generalized additive model

To model nonlinear but still interpretable relationships between the dependent variable and the explanatory variables, one can apply transformations (e.g., logarithm, categorization) to the original variables. Another option is to use GAM [START_REF] Trevor | Generalized additive models[END_REF], an extension of GLM where the linear predictor embeds unknown smooth functions of m explanatory variables:

g(E[Y |X 1 , ..., X m ]) = β 0 + m i=1 f i (X i ) (2.2)
To learn the nonlinear flexible functions f i , different approaches have been proposed:

Spline basis The original one is based on spline basis and makes use of the backfitting algorithm [START_REF] Trevor | Generalized additive models[END_REF], which works by iterative smoothing of partial residuals. However, it is difficult to estimate the degree of smoothness of splines f i . This parameter can be either set by the users or estimated with cross-validation during model training. However, this original approach is computationally expensive (O(n 3 ) with n the number of observations). Thus, several methods try to reduce the size of the basis used for smoothing [START_REF] Fahrmeir | Bayesian inference for generalized additive mixed models based on markov random field priors[END_REF][START_REF] Kim | Smoothing spline gaussian regression: more scalable computation via efficient approximation[END_REF], or find sparse representations of the smooth functions using Markov random fields [START_REF] Rue | Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations[END_REF].

Modern machine learning techniques In [START_REF] Tutz | Generalized additive modeling with implicit variable selection by likelihood-based boosting[END_REF], the authors proposed to train GAM with a likelihood-based boosting procedure. Lately, Lou et al. [START_REF] Lou | Intelligible models for classification and regression[END_REF], in a framework called Explainable Boosting Machine (EBM), proposed a tree-based extension of GAM that introduces a small number of two-dimensional interactions:

g(E[Y |X 1 , ..., X m ]) = β 0 + f i (X i ) + f ij (x i , x j ) (2.3)
To learn the nonlinear function f , they combine bagging [START_REF] Breiman | Bagging predictors[END_REF], gradient boosting [START_REF] Jerome H Friedman | Greedy function approximation: a gradient boosting machine[END_REF] (further explained in section 3.2.2) and automatic detection of pairwise interactions.

During training, features are trained one at a time in round-robin fashion to mitigate the effects of multicollinearity. When applied to tabular data, EBM outperforms standard GAM based on regression splines and is often almost as accurate as flexible black-box models. Moreover, EBM is also considered to be interpretable because each feature contributes to the prediction in an additive way (see Eq. 2.3) thus simplifying the understanding of its influence on the dependent variable.

Post hoc interpretability

Highly flexible models usually reach better predictive performance on some datasets mainly because they manage to capture nonlinearities in the data. However, this comes along with more complex design that makes them act like black-boxes. Posthoc explanation tools are required to provide some degree of interpretability. Different approaches have been proposed such as global model-agnostic explanation methods, simulation-based methods or local explanation methods.

Global model-agnostic explanation methods

These post hoc methods are used to describe the average behavior of a machine learning model. They quantify the main effects and the interaction effects of the Partial Dependence plots First introduced by Friedman [START_REF] Jerome H Friedman | Greedy function approximation: a gradient boosting machine[END_REF], partial dependence plots (PDP) show the dependence between the dependent variable and a set of features of interest, marginalizing over the values of the complementary features:

pd(x S ) = E X C [ f (x S , X C )] = f (x S , x C )p(x C )dx C
where f is the ML model, x S are the features of interest and X C are the complementary features here treated as random variables. As we marginalize over the other features, we get a function that depends only on features in S [START_REF] Molnar | Interpretable Machine Learning[END_REF]. We obtain the PDP by computing the integral for various values of x S . This plot provides insights into the nature of the relationship between the target and a feature (see Fig. 2.2). However, one main limitation of PDP lies in its assumption of independence between features which may result in the transmission of misleading information.

Sobol indices

In [START_REF] Ilya | Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates[END_REF], based on a finite expansion of a multivariate function f and under several constraints to maintain its unicity, Sobol was able to decompose the global variance V (Y ) into parts attributable to input features and combinations of features:

V (Y ) = n i=1 V i + n i n j>i V ij + • • • + V 12...n (2.4)
where

V i = V ( fi (x i )) = V (E[Y |x i ]) V ij = V ( fij (x i , x j )) = V (E[Y |x i , x j ]) -V i -V j . . .
From Eq. 2.4, one can estimate the main effect of a feature x i with the first-order Sobol indices, which is the amount of variance in the output explained by that feature:

S i = V (E[Y |X i ]) V (Y )

Simulation-based methods

These methods seek to obtain a global understanding of a black-box model, not by using global post hoc analysis tools, but by simulating its prediction function with a simpler one. Tan et al. [START_REF] Tan | Learning global additive explanations for neural nets using model distillation[END_REF] proposed to leverage model distillation techniques [START_REF] Buciluǎ | Model compression[END_REF][START_REF] Hinton | Distilling the knowledge in a neural network[END_REF] to learn global additive explanations of the form:

Ĥ(x) = h 0 + i h i (x i ) + i̸ =j h ij (x i , x j ) + i̸ =j j̸ =k h ijk (x i , x j , x k ) + . . . (2.5)
where h are either splines or bagged trees. During training, they seek to minimize the mean squared error between the prediction function of the black-box model and Ĥ(x).

If h are splines, they use penalized maximum likelihood and estimate the smoothing parameters with cross-validation. If h are bagged trees, cyclic gradient boosting [START_REF] Lou | Intelligible models for classification and regression[END_REF] is used. To understand the behavior of the model, they analyze the feature shapes by plotting the feature's contribution h i (x i ) against the domain of x i . Another approach, proposed by Lakkaraju et al. [START_REF] Lakkaraju | Faithful and customizable explanations of black box models[END_REF], use sub-space explanations to mimic the behavior of black-box models in classification problems. They leverage a two-level decision set representation, where the outer decision rules describe the subspaces (i.e. specific regions of the feature space), and the inner decision rules explain the decision logic of the black-box model within the related sub-spaces. Moreover, their approach allow end users to customize explanations by adding certain features of interest in the outer level representing the sub-spaces.

Local model-agnostic explanation methods

The objective of local explanation method is to learn local approximations to explain individual predictions. A broad variety of methods has been proposed and will be discussed in this section. First, we start by explaining how a contribution of a feature is computed on a simple linear model prediction and how it is extended to any model by means of Shapley values, a game theoretic approach. Then, we introduce LIME, a method to interpret individual model predictions based on locally approximating the model around a given prediction. Finally, we explain how Shapley values and LIME have been connected together to provide an efficient tool for local model-agnostic explanations.

Simple computation of contributions in a linear model setting A simple linear model prediction for an instance x is computed with:

f (x) = β 0 + p i=1 β i x i
where β i is the weight corresponding to variable i, and x i the value of variable i on this instance. The contribution ϕ i of the i-th feature can be computed by:

ϕ i ( f ) = β i x i -E[β i X i ]
where E[β i X i ] is the average effect of feature i [START_REF] Molnar | Interpretable Machine Learning[END_REF]. Thus, with a linear model, computing a feature's contribution to a prediction is straightforward. To measure how much a feature affects the prediction of more complex nonlinear models, a game theoretic approach based on the computation of Shapley values has been developed. [START_REF] Ls Shapley | Quota solutions op n-person games1[END_REF] measure how much feature i contributes to the prediction of any model f . Specifically, if S is a subset of the p features considered in the model, then the Shapley value ϕ i of a feature i is defined as:

Model-agnostic computation of a contribution with Shapley regression value Shapley regression values

ϕ i ( f ) = S⊆{1,...,p}\{i} |S|!(p -|S| -1)! p! (v(S ∪ {i}) -v(S)) (2.6)
where, for a given instance x, v x (S) is the prediction of feature values in S that are marginalized over features not included in S [START_REF] Molnar | Interpretable Machine Learning[END_REF]:

v x (S) = f (x 1 , ..., x p )dP x / ∈S -E X [ f (X)]
Local interpretable model-agnostic explanations (LIME) Ribeiro et al. [START_REF] Tulio Ribeiro | why should i trust you?" explaining the predictions of any classifier[END_REF] introduced LIME to approximate the model's behavior in the neighborhood of a given instance. For an instance x, they build a locally faithful and interpretable surrogate model g by solving:

ξ(x) = arg min g∈G L( f , g, π x ) + Ω(g) (2.7)
with L a loss function, f the model being explained, G a set of interpretable models (e.g., decision trees, linear models), and Ω a complexity measure. In Eq. 2.7, Ω can be perceived as a regularization term that constrains the surrogate model g to provide short explanations. Moreover, to ensure local fidelity, Ribeiro et al. define L as being the locality-aware square loss:

L( f , g, π x ) = z,z ′ π x (z)(f (z) -g(z ′ )) 2 (2.8)
where z ′ are binary vectors sampled at random from x ′ ∈ {0, 1} pviz., the interpretable representation of the instance x ∈ R p , p being the number of input features, which indicates if the feature i is present (x ′ i = 1) or missing (x ′ i = 0) in the surrogate model -and z are the original representations of z ′ . In Eq. 2.8, π x is a proximity measure defined by the authors as being the exponential kernel:

π x = exp(-d(x, z) 2 /σ 2 ),
with d a distance function and σ the kernel width.

LIME creates sparse linear explanations with a measure of faithfulness in the neighborhood of a given instance thanks to the introduction of the kernel π x . However, LIME highlights some limitations regarding the definition of this neighborhood. Defining a good neighborhood in an automatic manner with a proper kernel is still in development. For now, the selection of appropriate kernel settings is made heuristically. Also, as the neighborhood is computed based on points sampled uniformly at random, the analysis can be performed on unlikely data points. Moreover, recent work [START_REF] Slack | Fooling lime and shap: Adversarial attacks on post hoc explanation methods[END_REF] tend to show that it can be difficult to assess its reliability and robustness. In a theoretical study of LIME verified by simulations, Garreau and Luxburg [START_REF] Slack | Fooling lime and shap: Adversarial attacks on post hoc explanation methods[END_REF] showed that, while this post-hoc approach "discovers interesting features, it might forget some important features and the surrogate model is not faithful". Shapley additive explanation (SHAP) Lundberg and Lee [START_REF] Scott | A unified approach to interpreting model predictions[END_REF] unify existing approaches such as LIME [START_REF] Tulio Ribeiro | why should i trust you?" explaining the predictions of any classifier[END_REF], Shapley values [START_REF] Ls Shapley | Quota solutions op n-person games1[END_REF], DeepLIFT [START_REF] Shrikumar | Learning important features through propagating activation differences[END_REF], etc. SHAP quantifies the contribution of each original explanatory variable to each prediction. The authors proposed to represent Shapley value explanations as an additive feature attribution method. Let h x be an input mapping from simplified inputs x ′ , indicating feature presence, to original inputs x such that x = h x (x ′ ), and f a black-box model. An additive feature attribution model tends to reach:

g(z ′ ) ≈ f (h x (z ′ )) whenever z ′ ≈ x ′
Moreover, a model g in this class is a linear function of binary variables:

g(z ′ ) = ϕ 0 + M i=1 ϕ i z ′ i (2.9)
where the z ′ i variables indicate if the feature is observed (z ′ i = 1) or unknown (z ′ i = 0), ϕ i ∈ R are the feature attribution values, and M is the number of input features.

The authors state that there is a single unique solution that resolves this linear function while having the three following desirable properties [START_REF] Scott | A unified approach to interpreting model predictions[END_REF]:

• local accuracy: the explanation model g shall match the output of the original model f for the simplified input x ′

• missingness: a missing feature has no impact

• consistency: "if a model changes so that some simplified input's contribution increases or stays the same regardless of the other inputs, that input's attribution should not decrease" [START_REF] Scott | A unified approach to interpreting model predictions[END_REF] They prove that the coefficients ϕ of this linear function in fact corresponds to the Shapley values of the features:

ϕ i ( f , x) = z ′ ⊆x ′ |z ′ |!(M -|z ′ | -1)! M ! [ fx (z ′ ) -fx (z ′ \ i)]
with M the number of simplified input features, |z ′ | the number of non-zero entries in z ′ , and z ′ ⊆ x ′ represents all z ′ vectors where the non-zero entries are a subset of the non-zero entries in x ′ .

In the same paper, the authors introduce KernelSHAP, a model-agnostic approximation method that connects linear LIME and Shapley values. They prove that, given an appropriate weight function π x ′ (z ′ ), for an instance x, the Shapley values of the original features are the coefficients of the linear model g mimimizing the loss function defined in eq. 2.7:

L( f , g, π x ) = z ′ ∈Z [ f (h x (z ′ )) -g(z ′ )] 2 π x ′ (z ′ )
where:

π x ′ (z ′ ) = M -1 M |z ′ | |z ′ |(M -|z ′ |)

Model interpretability in the ML community

The machine learning research community offers nuanced perspectives about the merits of post-hoc explanations. As a representative example, Lipton [START_REF] Zachary C Lipton | The mythos of model interpretability: In machine learning, the concept of interpretability is both important and slippery[END_REF] suggests that post-hoc explanations should not be ruled-out as valid, although indirect, means of knowledge about the underlying data generating process. Lipton also underlines the potential risk of focusing on misleading information when relying on post-hoc explanations. Moreover, he considers that transparent linear models may not always be more interpretable than deep neural networks (DNN) because they often need heavily engineered features to obtain similar performances. Likewise, Poursabzi-Sangdeh et al. [START_REF] Poursabzi-Sangdeh | Manipulating and measuring model interpretability[END_REF] observe that practitioners can be affected by the information overload phenomenon [START_REF] Russell L Ackoff | Management misinformation systems[END_REF][START_REF] Keller | Effects of quality and quantity of information on decision effectiveness[END_REF] when the number of features becomes too large. Otherwise, Rudin [START_REF] Rudin | Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead[END_REF] emphasizes the importance of taking into account the whole data analysis process, including the preprocessing steps: "when considering problems that have structured data with meaningful features, there is often no significant difference in performance between more complex classifiers (DNN, boosted decision trees, random forests) and much simpler classifiers (logistic regression, decision lists) after preprocessing". Rudin points out that there is not necessarily a trade-off between accuracy and interpretability: performance gaps can be reduced iteratively through better data processing and model understanding. The latter is facilitated by the use of interpretable models.

Crash prediction models and their interpretability

In their survey, Lord and Mannering [START_REF] Lord | The statistical analysis of crashfrequency data: A review and assessment of methodological alternatives[END_REF] describe the different methods used for longterm crash frequency analysis. These methods can be classified in two categories: on one side, parametric statistical models explicitly associate the crash related variable to a vector of input explanatory variables. On the other side, non-parametric models use more complex design to better fit the data at the expense of their interpretability.

In this section, we describe both approaches, starting with the parametric statistical models.

Parametric statistical models

Poisson regressions

Crash-count data observations y i being positive integers, they have originally been modelled by Poisson regressions [START_REF] Jones | Analysis of the frequency and duration of freeway accidents in seattle[END_REF][START_REF] Sarath | Estimating truck accident rate and involvements using linear and poisson regression models[END_REF][START_REF] Miaou | The relationship between truck accidents and geometric design of road sections: Poisson versus negative binomial regressions[END_REF]. Poisson distribution is a special shape of the binomial distribution with a small probability of an event and a large but unknown number of trials. GLM have been used to link a linear predictor made of p explanatory variables to the rate λ i of a Poisson likelihood. As explained in section 2.1.2, an appropriate canonical link function is chosen to map the definition space of crash data (viz., discrete and positive only) to the real line covered by the linear predictor. For Poisson regressions, a log-link function is selected. The model can be formulated as:

y i ∼ P oisson(λ i ) log(λ i ) = β 0 + p j=1 β j x ij

Negative Binomial regression

The Poisson distribution has a unique free parameter, the rate λ i . The variance cannot be adjusted independently of the mean. Yet, the observed variance of crash count data often exceeds this amount. A negative binomial (NB) distribution can be proposed to better deal with this over-dispersion phenomenon. Given the probability of crash, the NB gives the probability of observing n crashes before the α-th non-crash. With λ representing the mean, its probability mass function can be parameterized as follows [START_REF] Walck | Hand-book on statistical distributions for experimentalists[END_REF]:

NB(n; λ, α) = n + α -1 n λ λ + α n α λ + α α
The variance of the NB is λ + λ 2 α . Therefore, the parameter α controls the amount of over-dispersion. When α → ∞, the NB likelihood approaches a P oisson(λ) distribution. Moreover, as shown in [START_REF] Walck | Hand-book on statistical distributions for experimentalists[END_REF], a NB(λ, α) corresponds to a P oisson(λ) where λ comes from a gamma distribution Gamma(a = α λ , b = α) whose probability density function is defined as:

Gamma(x; a, b) = a b x b-1 e -ax Γ(b)
for x > 0, Γ(b) being the gamma function.

The Poisson-gamma model -a continuous mixture of Poisson distributions with rates distributed as a gamma distribution -has been identified as a reference model by road safety experts (see for example the AASHTO Highway Safety Manual [START_REF]The highway safety manual[END_REF]). We find it in the core of many related work [START_REF] Miaou | Modeling traffic crash-flow relationships for intersections: dispersion parameter, functional form, and bayes versus empirical bayes methods[END_REF][START_REF] Lord | Modeling crash-flowdensity and crash-flow-v/c ratio relationships for rural and urban freeway segments[END_REF][START_REF] El | Comparison of two negative binomial regression techniques in developing accident prediction models[END_REF][START_REF] Lord | Examining the effects of site selection criteria for evaluating the effectiveness of traffic safety countermeasures[END_REF] where crash count regression is done with a linear model attached to the λ parameter of the NB through a log-link function.

Generalized linear mixed model

In some studies, correlation among observations can be taken into account on the temporal level. This arises in the context of panel data analysis, where crash data are considered to be collected through a series of repeated observations of the same groups (e.g., roadway segments) over a time period. To account for the differences between segments, GLM are accommodated to introduce random effects in addition to the usual fixed effects. The resulting model is known as a generalized linear mixed model (GLMM). In this configuration, Poisson regressions are adapted as follows:

y ik ∼ P oisson(λ ik ) log(λ ik ) = β 0 + p j=1 β j x ij + u k
where u k is the random effect for group k, supposed to be distributed according a zero-centered Gaussian distribution with variance σ 2 (σ > 0).

GLMM have found many applications in road safety analysis. Amongst them, Johansson et al. [START_REF] Johansson | Speed limitation and motorway casualties: a time series count data regression approach[END_REF] analyzed the effect of a lowered speed limit on the number of crashes on Swedish highways. In a study of crashes caused by median crossovers in Washington State, Shankar et al. [START_REF] Venkataraman | Evaluating median crossover likelihoods with clustered accident counts: An empirical inquiry using the random effects negative binomial model[END_REF] compared a standard NB and a NB modified to account for random effects relating to each site. However, they did not observe any benefits from using a NB with random effects.

Hierarchical model

GLM can also be refined into hierarchical models to take into account clusters of related observations. For example, crashes occurring in a given geographical region may possess specific characteristics while not differing entirely from crashes in other regions. A simple Poisson regression, by pooling all the observations together, would assume an invariant population and couldn't benefit from regional peculiarities. Otherwise, k clusters could be modeled with the addition of k -1 mutually exclusive binary variables to the linear model, but this would correspond to no pooling at all and the clusters would be assumed independent of one another. Contrariwise, hierarchical models, also known as multilevel models, offer partial pooling through an adaptive regularizing prior. Thus, in the following multilevel Poisson regression with k clusters, hyperpriors µ and σ will allow an adaptive shrinkage of the cluster-specific β jk towards a common mean: level 0 -model

y ik ∼ P oisson(λ ik ) log(λ ik ) = β 0k + j β jk x ijk level 1 -priors β jk ∼ N (µ, σ) level 2 -hyperpriors µ ∼ N (0, 100) σ ∼ HN (100)
with HN being the half-normal distribution. Gelman indicates in [START_REF] Gelman | Prior distributions for variance parameters in hierarchical models (comment on article by browne and draper)[END_REF] that a halfnormal distribution with high standard deviation is a non-informative but proper prior for the variance parameter σ.

In this way, Jones and Jørgensen [START_REF] Andrew | The use of multilevel models for the prediction of road accident outcomes[END_REF] design a multilevel model to predict the severity of an incident given the involved casualties (level 1), their respective vehicles (level 2) and the accident location (level 3). Ahmed et al. [START_REF] Ahmed | Exploring a bayesian hierarchical approach for developing safety performance functions for a mountainous freeway[END_REF] conceive a multi-level model to predict crashes on a mountainous freeway by modeling both the dry or snow seasons and spatial correlation between adjacent sites. Deublein et al. [START_REF] Deublein | Prediction of road accidents: A bayesian hierarchical approach[END_REF] propose a multilevel model to manage simultaneously 4 response variables (resp. injury accidents, light injuries, severe injuries and fatalities) through gamma updating of the parameters. Finally, Fawcett et al. [START_REF] Fawcett | A novel bayesian hierarchical model for road safety hotspot prediction[END_REF] predict future safety hotspots with a multilevel model where, first, the variance of the rate of a NB increases with the timestamp of an observation and, second, a global trend effect is altered by site-specific ones.

Intrinsic interpretability

The aforementioned models are interpretable by nature. Indeed, the relation between the explanatory variables and the dependent one is explicit. Nevertheless, as Poisson regressions and their variants use a log-link function to map the linear predictor to the crash-related variable, the effects of the coefficients β are not as obvious as those of a linear model. To mitigate this, Kweon and Kockelmam [START_REF] Kweon | Safety effects of speed limit changes: Use of panel models, including speed, use, and design variables[END_REF] proposed to use the incidence rate ratio which measures the percentage of change in the crash related variable when an explanatory variable is increased by a unit.

Moreover, analyzing the effects of coefficients shall not be based solely on point estimates but must come along with the knowledge of the uncertainty associated with them. Most of studies mitigate this by using Bayesian inference. By associating a probability distribution to each coefficient, they can measure the underlying uncertainty of coefficients.

Non-parametric Machine Learning models

Support vector machines SVM were originally introduced by Cortes and Vapnik [START_REF] Cortes | Support-vector networks[END_REF] for classification tasks. A version for regression has been developed by Drucker et al. [START_REF] Drucker | Support vector regression machines[END_REF]. Among the variants proposed for regressions, the epsilon-insensitive SVM (ε-SVM) is the most common. Given a vector y ∈ R n of dependent variables and X ∈ R n×p the associated explanatory variables, ε-SVM aims to find a function f (x i ) that deviates from y i by a value inferior to ε for each training point x i , i ∈ {1, ..., n} (see Fig. 2.3). In addition, f has to be as flat as possible. To achieve this, ε-SVM solves the following problem:

min w,b,ζ,ζ * 1 2 w T w + C n i=1 (ζ i + ζ * i ) (2.10) subject to y i -w T ϕ(x i ) -b ≤ ε + ζ i , w T ϕ(x i ) + b -y i ≤ ε + ζ * i , ζ i , ζ * i ≥ 0, i = 1, ..., n
where ζ i , ζ * i penalize observations lying outside the ε margin and C is a positive regularization constant that controls the penalty imposed by ζ i and ζ * i . In this problem, ϕ is a transformation function that maps x to a higher dimensional space where the problem resolution is linear. Without it, ε-SVM won't have the ability to model nonlinear relationship between explanatory variables and the dependent variable in the original space.

The optimization problem described above is computationally simpler to solve by constructing a Lagrangian function which introduce positive multipliers α i and α * i for each observation x i [START_REF] Drucker | Support vector regression machines[END_REF]:

min α,α * 1 2 (α -α * ) T G(α -α * ) + ε1 T (α + α * ) -y T (α * -α) (2.11) subject to 1 T (α -α * ) = 0 0 ≤ α, α * ≤ C, i = 1, ..., n
In eq. 2.11, G is a n × n positive semidefinite matrix whose elements are computed with the kernel function k:

G ij = k(x i , x j ) where k(x i , x j ) = ⟨ϕ(x i ), ϕ(x j )⟩ (⟨•, •⟩ being the inner product).
Here, the explicit mapping ϕ can be avoided as we can obtain G by computing G ij directly with the kernel function. This process is known as the kernel trick. By doing so, ε-SVM estimates f in the transformed input space and not the original one. To generate this transformation, various kernels have been proposed: linear degree of the polynomial), or based on radial basis function (k

(k(x i , x j ) = x i x j ), polynomial (k(x i , x j ) = (1 + x T i x j ) d with d the
(x i , x j ) = exp(-||x i - x j || 2 /σ 2 ).
Finally, to make a prediction, ε-SVM computes the contribution of each support vector (i.e. a vector x i for which α i or α * i is equal to zero):

f (x) = i∈S (α i -α * i )k(x i , x) + b
where S is the set of support vectors.

In the binary classification problem, the objective is not to find an hyperplane that contains a maximum number of observations within its ε margin, but to determine the best decision boundary between the two classes. In this case, the goal is to maximize the margin between the closest vectors of each class and the hyperplane.

In highway safety analysis, SVM have been applied to various problems. Li et al. [START_REF] Li | Predicting motor vehicle crashes using support vector machine models[END_REF] compare a NB and a SVM on crash-count predictions for rural roadway segments located in Texas, USA. They observe that SVM predict crash data more effectively and accurately than traditional NB. In a study on crash data collected at 326 freeway diverge areas, Li et al. [START_REF] Li | Using support vector machine models for crash injury severity analysis[END_REF] apply SVM to the task of predicting crash severities. When compared to an ordered probit model, they found that SVM produces better prediction performance. Lately, Dong et al. [START_REF] Dong | Support vector machine in crash prediction at the level of traffic analysis zones: assessing the spatial proximity effects[END_REF] study the effects on crash-count predictions of spatial correlations at different scales on a dataset including more than 57 000 crashes in the Hillsborough county of Florida.

Artificial neural networks

Artificial neural networks (ANN) are biologically inspired computational networks first introduced by McCulloch [START_REF] Warren | A logical calculus of the ideas immanent in nervous activity[END_REF]. ANN are composed of artificial neurons typically organized into a multi-layer structure. Neurons of one layer are connected to neurons of the next layer. An ANN generally consists of three types of layers: the input layer receives the external data, then the hidden layer(s) process the data and finally, the output layer produces a result (e.g. crash-count prediction). In the well-known configuration of multilayer perceptrons (MLP), neurons are fully connected i.e., every neuron in one layer is connected to every neuron in the next layer (see Fig. 2.4a). Each neuron has inputs and produces a single output by first computing weighted sum of its inputs and then applying a nonlinear activation function σ such as sigmoid or hyperbolic tangent (see Fig. 2

.4b).

When training an ANN, the objective is to adjust the connection weights to minimize the cost function that measures the overall performance of the network. One of the most common approach is backpropagation: first, it calculates the gradient of the cost function associated with a given state with respect to the weights, and then updates the weights in the opposite direction of the gradient, as it is the direction of the steepest descent. Several variants of gradient descent have been developed and are well explained in [START_REF] Ruder | An overview of gradient descent optimization algorithms[END_REF].

ANN have found numerous applications in highway safety analysis. Abdelwahab and Abdel-Aty [START_REF] Hassan | Development of artificial neural network models to predict driver injury severity in traffic accidents at signalized intersections[END_REF] compares a backpropagation neural network (BPNN) with an ordered logit statistical model to predict the severity of accidents at intersections. Chang [START_REF] Chang | Analysis of freeway accident frequencies: negative binomial regression versus artificial neural network[END_REF] compares a one hidden layer BPNN with NB regression for crash frequencies prediction. BPNN slightly outperforms the NB model with a difference of 0.6% of accuracy on testing data. Huang et al. Post-hoc interpretability SVM and ANN are opaque decision systems. The former, by introducing the kernel function, models complex relationships in a high-dimensional space. The later often involves a large number of parameters growing exponentially with the number of hidden layers. To reduce the complexity of the models, one can restrict the kernel of SVM to be linear or define a simpler architecture for the ANN (e.g. a perceptron with a single hidden layer with few neurons). However, trying to reduce complexity in this way will tend to deteriorate the predictive results. To balance performance and interpretability, all studies of the previous section opted instead for the use of posthoc explanation methods. Authors perform a sensitivity analysis of the black-box model: for each explanatory variable, while keeping all other variables unchanged, they record the effect on the output prediction of a perturbation of the current variable (see Fig. 2.5). However, as stated in [START_REF] Xie | Predicting motor vehicle collisions using bayesian neural network models: An empirical analysis[END_REF], the relationship between the current explanatory variable and crash frequency may vary due to correlated variables, making it difficult to interpret the sensitivity plots. Moreover, by generating simulated data while assuming the explanatory variables to be independent, sensitivity analysis will indiscriminately generate potentially misleading hypothetical predictions for unlikely data points. Finally, this method only provides global interpretations and does not allow domain experts to identify roadway segments where the predictive model behaves singularly. To mitigate this, some studies used local model-agnostic explanation methods. Mihaita et al. [START_REF] Mihaita | Arterial incident duration prediction using a bi-level framework of extreme gradient-tree boosting[END_REF] investigate with SHAP the impact of different features on arterial incident duration. They observe that the number of affected lanes and the hour of the day seem to be the most important features which increase the incident duration prediction. Parsa et al. [START_REF] Amir Bahador Parsa | Toward safer highways, application of xgboost and shap for real-time accident detection and feature analysis[END_REF] extend the use of SHAP by analyzing complex and nonlinear joint impacts of features on the output of a XGBoost model, in the context of real-time accident detection.

Chapter 3 Contributions

Overview

Introduction

In previous chapters, we stressed the need to be able to provide explanations in addition to predictions when actions are pondered based on the analysis of predictive models. To obtain such explanations, we saw that two approaches have been developed: intrinsically interpretable models that offer a direct understanding of their internal decision process, versus post-hoc analysis of black-box models. However, we have seen that using post-hoc explanation methods comes with some limitations. For data-driven high-stakes decisions, it is imperative that the explanations transparently and truthfully reflect what the model has learned. Moreover, as observed by Rudin [START_REF] Rudin | Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead[END_REF], simple and interpretable models can achieve predictive performance similar to black-box models on well-structured tabular data. For these reasons, we favored the use of simple interpretable models.

The focus of our work is to enhance the predictive power of simple models without deteriorating their high degree of interpretability. Our main contributions aim to achieve this goal in two steps. First, we introduce a supervised method to discover a partition of the original observations and build a hierarchical model above it. Second, we introduce two algorithmic approaches (viz., a polynomial neural network, and an extension of multi-objective symbolic regression) to discover highly discriminant non-linear transforms of the original variables. The former can handle correlations among groups of observations which usually lead to improvements in the quality of the models' predictions and of their interpretation. The latter, while remaining simple (e.g. first-order interactions), allow the models to capture more of the variability in the dependent variable.

Bayesian hierarchical generalized linear model

In our first proposal presented in Fig. 3.1, we describe how to extract, from data and without prior expert knowledge, the information necessary to build a hierarchical Bayesian model with first-order interactions between explanatory variables to efficiently solve regression or classification problems while preserving interpretability. First, in the hierarchical structure module described in section 3.2, we elucidate how Shapley values of the variables give rise to a clustering of the original observations that is likely to make sense in terms of the problem to be solved. This clustering then informs a multilevel Bayesian model (see section 3.3). Second, we explain how a self-organized neural network reveals the most important interactions between explanatory variables (viz., first-order interactions module). These interactions are then integrated to the functional form of the multilevel model (see section 3.3).

Interpretable hierarchical symbolic regression

In our second proposal (see Fig. 3.2), the discovered hierarchical structure is even better exploited by using symbolic regression to capture sparser, more complex but still interpretable relationships between the explanatory variables and the crash count. To do so, to the whole training dataset and to each cluster of the hierarchical structure, we apply a variant of the symbolic regression (SR) method to find expansions of linear models with effective and interpretable functional forms. To this end, presented in section 3.4, we designed a multi-objective simulated annealing algorithm to solve the SR problem. Thus, we can discover Pareto optimal predictive models1 with various trade-offs between accuracy and complexity. In our case, symbolic regression serves two purposes. First, it discovers global models, learned from the whole training 3.2 Supervised learning of a hierarchical structure

Introduction

As explained in section 2.2.1, a well-chosen hierarchical structure can handle correlations among groups of observations and significantly improve the quality of the models' predictions and of their interpretations. However, in the literature, this structure, which is dependent on the dataset, is the result of either expert knowledge or unsupervised clustering. The former is not always available, and the latter does not account for the task being solved (e.g., prediction of crash counts) and tries to group features without knowing if they're relevant for an outcome of interest [START_REF] Scott M Lundberg | From local explanations to global understanding with explainable ai for trees[END_REF]. Moreover, in our work, we would like to group roadway segments that share a similar crash severity for similar reasons. In other words, when computing clusters, we do not want to focus only on the distribution of the crash-related variable, but also want to account for the influence of each explanatory variable. Indeed, a similar crash severity between segments may be explained by different risk factors. Thus, clustering segments by considering only the dependent variable may result in segments in a same cluster that are in fact hazardous for different reasons.

To obtain such a hierarchical structure, we propose to learn it from the data in two stages. First, we analyze the results of a black-box machine learning model with a local post hoc explanation tool (viz., SHAP). More specifically, for each instance of the dataset, SHAP computes the contribution of each explanatory variable to the prediction. Then, we apply a clustering on these SHAP explanation instances. In other terms, we make use of an unsupervised clustering where samples are grouped together based on their explanations (viz., the SHAP explanation instances), thus alleviating the aforementioned shortcomings of unsupervised clustering.

In this section, we describe how we successfully retrieve a relevant hierarchical structure with an unsupervised clustering on labeled datasets2 . First, we start by presenting gradient boosting trees, the models we selected to discover a latent structure in the data. Then, we explain how the black-box effect induced by these models is mitigated by leveraging a local explanation tool (viz., SHAP) to measure the features' contributions to each observation. Finally, we explain how a hierarchical agglomerative clustering applied to these SHAP explanation instances help to discover a relevant structure made up of similarities of explanations between observations.

Training of a non-parametric ML model

Model selection

We are interested in highly flexible, efficient models with fast training time on tabular data. To answer the first two requirements, numerous models can be considered: nonlinear SVM, ensemble of trees, deep learning models, . . . However, training a nonlinear SVM is computationally infeasible on large dataset (between O(n 2 ) and O(n 3 ), n being the number of samples, for a standard implementation of SVM [START_REF] Abdiansah | Time complexity analysis of support vector machines (svm) in libsvm[END_REF]) which make them poorly suited to the preliminary task of discovering a hierarchical structure in the data. On their side, deep learning models are better adapted to image recognition, natural language processing or speech recognition [START_REF] Scott M Lundberg | From local explanations to global understanding with explainable ai for trees[END_REF].

In their experiments, Chen et al. [START_REF] Chen | Xgboost: A scalable tree boosting system[END_REF] observe that tree-based models consistently outperform standard deep learning models on tabular datasets where "features are individually meaningful and do not have a strong multi-scale temporal or spatial structures". Recently, Borisov et al. [START_REF] Borisov | Deep neural networks and tabular data: A survey[END_REF] compare algorithms based on gradientboosted trees with a broad variety of deep learning models on three middle and large size datasets covering two classifications and one regression task. They observe that gradient-boosted trees outperform deep learning models on all datasets, while having lower training time.

A balance of computational efficiency, ease of use, and high accuracy have made tree-based models the most popular non-linear model type. In 2021, 75% of most popular ML models used by data scientists and engineers are based on trees [START_REF] Kaggle | The state of ml and data science[END_REF]. Among ensemble of trees, random forest or gradient boosting trees are the most popular. However, the former tends to reduce only the variance with a bagging strategy while the latter, by combining multiple models and gradient descent on residuals, reduces both the variance and bias. Plus, gradient boosting trees often outperform random forests [START_REF] Borisov | Deep neural networks and tabular data: A survey[END_REF]. For these reasons, we make use of gradient boosting trees as a first step to discover the hierarchical structure.

Gradient tree boosting

Given a dataset {(x i , y i )} n i=1 , gradient tree boosting aims to build, through additive training, an ensemble of trees T that minimizes a loss function L (e.g., mean squared error for regression tasks, logistic loss for classification tasks). The algorithm starts by initializing the model with the most accurate constant α:

T 0 (x) = arg min α n i=1 L(y i , α)
Then, at each iteration m, the model is expanded in a greedy fashion such that:

T m (x) = T m-1 (x) + arg min tm n i=1 L(y i , T m-1 (x i ) + t m (x i ))
where t m is a decision tree. In other words, the model T m tries to correct the error of its predecessor T m-1 by adding a new decision tree, also called weak learner. However, choosing the decision tree that best reduces the loss function at each iteration has a high computational cost. Thus, assuming L is differentiable, the optimization process is simplified by leveraging gradient descent. The new decision tree t m is trained on {(x i , r im )} n i=1 with r im the residuals indicating locally the steepest direction. Specifically, r im are defined as:

r im = - ∂L(y i , T (x i )) ∂T (x i ) T (x)=T m-1 (x) , i = 1, . . . , n
Then, t m is added to the expansion such that:

T m (x) = T m-1 (x) + γ m (t m (x))
where γ m is selected by minimizing the loss function of the new expansion:

γ m = arg min γ n i=1 L(y i , T m-1 (x) + γ m (t m (x)))
The expansion process stops when a stopping criterion is reached (e.g., maximum number of iteration, satisfactory value of the loss).

Implementation

We select the XGBoost3 's implementation of gradient tree boosting. The major benefits of using this library lie in its high portability, scalability, and high performance as it provides state-of-the-art results on many problems [START_REF] Chen | Xgboost: A scalable tree boosting system[END_REF][START_REF] Borisov | Deep neural networks and tabular data: A survey[END_REF].

Moreover, the authors use three different methods to prevent overfitting [START_REF] Chen | Xgboost: A scalable tree boosting system[END_REF]. First, they penalize the complexity of the model when optimizing the loss function. Their definition of model complexity penalizes both the number of leaves in the trained trees, and the leaf values by means of a l2-regularization term. Then, they add shrinkage, a technique proposed by Friedman [START_REF] Jerome H Friedman | Stochastic gradient boosting[END_REF], to control the influence of new added trees during the update process. Finally, they use features subsampling to lower the number of features available to each new tree and thereby increase the variance between trees.

Local post hoc explanations

The aforementioned model is an ensemble of trees which consists of hundreds, if not thousands, of trees grown sequentially (in our experiments, the number of trees is set to 100, the default value of the XGBoost package). Even though the authors penalize the complexity of each tree while training, the resulting model is hardly interpretable. To mitigate this black-box effect, we compute local explanations with SHAP [START_REF] Scott | A unified approach to interpreting model predictions[END_REF] (see section 2.1.3), to understand how the model maps the original explanatory variables to an outcome (e.g., a crash count in the case of highway safety). In addition to retaining local faithfulness, SHAP presents important properties described in section 2.1.3. To each observation, SHAP associates a linear function g:

g(z ′ ) = ϕ 0 + M i=1 ϕ i z ′ i (3.1)
where M is the number of input features, z ′ i ∈ {0, 1} M and ϕ i ∈ R are their contributions which correspond to the game theoretic concept of Shapley values.

Among the different implementations of SHAP, we selected TreeSHAP, an efficient tree-based algorithm for fast and consistent computations of exact Shapley values [START_REF] Scott M Lundberg | Consistent individualized feature attribution for tree ensembles[END_REF][START_REF] Scott M Lundberg | From local explanations to global understanding with explainable ai for trees[END_REF]. Compared to KernelSHAP, TreeSHAP reduces the complexity of Shapley value computation from exponential to low order polynomial time [START_REF] Scott M Lundberg | Consistent individualized feature attribution for tree ensembles[END_REF]. Aside from ensuring the properties defined in section 2.1.3 (viz., local accuracy, missingness and consistency), TreeSHAP is also able to account for feature dependence. Moreover, in [START_REF] Scott M Lundberg | From local explanations to global understanding with explainable ai for trees[END_REF], the authors observed that TreeSHAP consistently outperforms alternative methods across a benchmark of 21 different local explanation metrics. Furthermore, in Lundberg and Lee [START_REF] Scott | A unified approach to interpreting model predictions[END_REF], the authors propose a force plot visualization (see Fig. 3.3) to materialize how much each contribution shifts the output relatively to the overall expected value of the XGBoost model. These contributions can be perceived as being forces that shift the output towards or away from the expected value of the model (hence the name of the plot): the higher the contribution, the larger the magnitude of the force is in Fig. 3.3. In practice, these plots help us understand visually how the black-box model behaves on single instances of our dataset. In the next section, we will describe how they also reveal a hierarchical structure in the data.

Discovery of a hierarchical structure

When applied to all observations, the SHAP forceplots can be clustered by similarities of their profiles. On the highway network dataset, we discover clusters of roadway segments which are similar based on their SHAP explanation instances (see Fig. 3.4): the left part of the figure characterizes roadway segments that are on average moderately hazardous, while the middle and far right parts represent highly hazardous segments.

To obtain such a structure in our tabular dataset, we use a hierarchical agglomerative clustering (HAC) of the observations based on the explanatory variables' contributions as provided by the SHAP analysis. We select the squared euclidean distance as a measure of distance between pairs of instances and the ward linkage criterion to measure the dissimilarity of groups of instances. Unlike other distance measures, one characteristic of the ward linkage criterion lies in the fact that it considers the inner-variability of clusters when computing distances between two clusters. More specifically, given A and B two clusters, it states that the distance d between these

d(A, B) = i∈A∪B |x i -m A∪B | 2 - i∈A |x i -m A | 2 - i∈B |x i -m B | 2 = n A n B n A + n B |m A -m B | 2
where | • | 2 is the norm of the vector, m J the center of cluster J and n J the number of SHAP instances in it. At the beginning, each cluster consists of a unique instance of SHAP and d is equal to zero. Then, at each iteration, clusters are merged together in a way that minimizes the increase in d. The procedure is repeated until a unique cluster is obtained, composed of the whole set of SHAP instances. Finally, we estimate the optimal number of clusters by detecting the greatest increase in the squared Euclidean distance between clusters when their number decreases (see Thorndike [START_REF] Robert | Who belongs in the family[END_REF] for the original presentation of this widely used method). Furthermore, in production, we do not have a straightforward way of knowing to which cluster a new observation belongs to. However, we need this information as our underlying predictive models require this prior knowledge. Thus, to learn how to associate an observation of the test dataset4 to a cluster, we propose to train a decision tree classifier to approximate a good mapping between the original explanatory variables and the cluster of SHAP instances. In section 4.5.1, we report on cross-validation measures showing, on various datasets, that this association is very accurate.

Representation of the hierarchical structure in highway safety

In this section, we suppose that a hierarchical structure has been discovered by using ten years of data, from January 1st, 2008 to December 31th, 2017. Four clusters have been identified and are presented in Fig. 3.5. To obtain additional information on these automatically discovered clusters, we propose to visualize the distributions of explanatory variables and the dependent variable. In Fig. 3.6a and Fig. 3.6b, we illustrate this for the observed crash counts and the annual average daily traffic (AADT) which allows us to understand that clusters do not have the same amount of historical crash counts, and present as well differences regarding the explanatory variables such as the traffic related one. If necessary, experts can fine-tune the number clusters with the help of a dendrogram to obtain clusters that better match with their domain knowledge. They may be interested in broader analysis (with less clusters) or more detailed analysis (with more clusters) that highlight particular group of roadway segments.

When reducing the number of clusters from 4 to 2, one remains identical (viz., Cluster 1 of Fig. 3.7) while the others are merged together. The difference between the two resulting clusters mainly lies in the Shapley values of the AADT, which explains most of the critical nature of cluster 1. When looking at the historical data, we observe that this cluster is more hazardous than cluster 0 (see Fig. 3.8a) and is composed of roadway segments with high traffic (see Fig. Conversely, when increasing the number of clusters, some particular hazardous configurations, explained differently by the SHAP analysis, appear on the highway network. For instance, cluster 0 of Fig 3 .5, a moderately hazardous cluster composed mainly of rural and mountainous segments, is now divided into cluster 1 and 4 in Fig. 3.9. The most hazardous segments of the previous cluster (see Fig. 3.6a) are now grouped into cluster 4. These segments are located in the mountainous part of the network (see Fig. 3.10b).

However, moving away from the initial clusters defined by the automatic process is not without risks, especially when the number of clusters is increased. Indeed, increasing the number of clusters will cut the explanation space in a way that some clusters will only be represented by a few samples which will bring a high source of uncertainty in the models. Moreover, we observed in our experiments that the performance is altered when the number of clusters moves away from the original value defined in the supervised way. We recommend to stay close to the original value as it represents a good trade-off between performance and interpretability. Finally, road safety experts confirm the relevance of these automatically discovered clusters. They also point out that the time saved can advantageously be spent on, for example, planning remedial actions.

Bayesian hierarchical generalized linear model

In our first proposal, the hierarchical structure discovered in the previous chapter is used to enhance the predictive performance and interpretability of generalized linear models (GLM). To account for such a structure, we use a Bayesian inferred hierarchical GLM. In section 3.3.1, we start by describing the multilevel structure of such models and then explain how Bayesian inference, with advances in probabilistic programming, infers efficiently the coefficients of these models (see section 3.3.2).

Moreover, as the standard formulation of GLM generally attaches a linear functional form to a parameter of the likelihood distribution, interactions between the explanatory variables are not taken into account. In section 3.3.3, we also describe how we mitigate this by first learning simple interactions from the analysis of the results of a special kind of polynomial neural network, and then by adding them to the linear functional form.

We finish this chapter by explaining how we evaluate and validate the models with a well-known technique, the Posterior Predictive Check.

Model description

To integrate the hierarchical structure, made of k clusters, into a GLM, we design the following multilevel model: level 0 -model

Y ik ∼ L(Y ik |µ ik , Ω) µ ik = E[Y ik |X ik ] = g -1 (η ik ) η ik = β 0k + N j=1 β jk X ijk level 1 -priors β jk ∼ N (µ, σ) level 2 -hyperpriors µ ∼ N (0, 100) σ ∼ HN (100)
According to this model, output Y ik , for observation i in cluster k, is generated from a likelihood distribution L parameterized with a set of specific parameters Ω and also µ ik , the expected value of Y ik conditioned on the observations. In section 2.2.1, we saw that, for crash prediction, L can be a N egativeBinomial(Y ik |λ, α) where λ is the expected number of crashes and α controls the amount of allowed overdispersion. Next, η ik is a linear transformation of the explanatory variables. The inverse link function g -1 is necessary to map the domain of η ik (viz., the real line) to the one of µ ik . For example, since the rate λ of a negative binomial must be positive, the exponential function is used for crash-count prediction. In case of binary classification, when the Bernouilli(p) likelihood is used, p being a probability, g can be set to the logit function. Indeed, the logit maps a parameter constrained between 0 and 1 onto the real line (the inverse link function g -1 is, in that case, the logistic function). Finally, as explained in the previous section, the cluster-specific coefficients β jk , for the N explanatory variables, depend on hyperpriors µ and σ, thus allowing an adaptive shrinkage to a mean common to all the observations.

Bayesian inference

Description

Bayesian Inference (BI) treats the parameters as random variable and associates a probability distribution to each parameter. Specifically, BI derives the posterior probability from the prior p(θ) and the likelihood p(y | θ):

p(θ | y) ∝ p(y | θ)p(θ) (3.2)
with θ the set of parameters in the models. Based on the knowledge of observed crash-count data, the parameters are updated according to eq. 3.2. However, the posterior distribution is usually not obtained in a closed-form distribution but by means of approximation techniques. In our work, we estimate the distribution with a sampling procedure based on a Markov Chain Monte Carlo (MCMC) algorithm. Despite the use of computationally expensive sampling algorithms, BI becomes more and more efficient thanks to the advances of probabilistic programming [START_REF] Development | Stan: A c++ library for probability and sampling, version 2.5[END_REF][START_REF] Salvatier | Probabilistic programming in python using pymc3[END_REF].

Note that we could have used the frequentist alternative based on Maximum likelihood estimation (MLE). MLE use iterative algorithm, such as the Newton-Raphson method or gradient descent, to estimate coefficients β j . MLE is closely connected to the Maximum A Posteriori (MAP) estimate, which is the mode of the posterior distribution after BI. Indeed, given a uniform prior, MLE and MAP estimates are identical. However, with MLE, each coefficient has a unique point estimate, thus the introduction of prior knowledge is not possible. Moreover, with BI, we also obtain uncertainty knowledge for each MAP estimate. For these reasons, we use BI to infer the parameters of our models.

Uncertainty knowledge

The analysis of the posterior distributions (see Fig. 3.11) allows one to measure the influence of each variable on the crash count. Due to the inherent partial pooling nature of the model, posterior distributions are dissimilar among clusters thus revealing various impacts of the same explanatory variable on the crash count. For instance, the posterior's mean related to the percentage of heavy vehicles in Fig. 3.11b is positive for cluster 0 and 3, but negative for cluster 1.

Moreover, we observe that the shapes of the posterior distributions are different from one cluster to another. This variability is linked to the size of the clusters: in general, the more samples, the more confidence in the estimates. For example, in Fig. 3.11a the posterior of the traffic related variable (log(AADT)) has a sharp probability distribution for cluster 0 and 3 but a flatter one for cluster 1 and 2. The latter highlights a greater uncertainty in estimating the coefficient associated with the traffic and calls for more vigilance when drawing conclusions for this risk factor.

Supervised learning of interactions with a polynomial network

Objectives

When first-order interactions between explanatory variables are integrated into a GLM, the relationship between a variable and the target may depend on the value of another variable. This can lessen the gap in predictive power between Bayesian inferred GLM and ML algorithms inherently able to capture complex nonlinear relationships. Moreover, these simple interactions are interpretable while the potentially highly entangled ones discovered by ML algorithms will often remain inaccessible to human understanding and increase the risk of overfitting.

Group Method of Data Handling algorithm

In our approach, important first-order interactions are discovered with a variant of the Group Method of Data Handling (GMDH) family of supervised algorithms [START_REF] Grigorevich | Polynomial theory of complex systems[END_REF]. This GMDH algorithm is a self-organized multi-layered structure of nodes (see Fig. 3.12). Each node generates its output z by applying a linear function with a covariation term to a pair of inputs (x 1 , x 2 ) taken among either the nodes of the previous layer or the original explanatory variables:

z = a 0 + a 1 x 1 + a 2 x 2 + a 3 x 1 x 2
Let n be the number of nodes of the previous layer and m be the number of explanatory variables. To build the next layer, for each of the n+m 2 polynomials, the a 0 , ..., a 3 parameters are set by minimizing through Ridge regression the least square error made by the polynomial when it approximates the target on a train dataset. Then, the fitted polynomials are evaluated on a validation dataset to select the top m constituting the new layer. When the score obtained on the validation dataset by the best node of the last layer added stops improving, the process terminates and the best node of the penultimate layer is the output of the network. Thus, this self-organized network discovers a polynomial that approximates the relationship observed on the training dataset between the explanatory variables and the target.

In our approach, we use this polynomial to select the most important first-order interactions between explanatory variables. If the coefficient of a term involving the product of two variables exceeds a given percentage δ of the magnitude of the biggest coefficient, δ being a hyper-parameter of our framework, we add to our Bayesian hierarchical model an interaction between these two variables.

Integration of interactions into the GLM

For a given observation i, let {int i1 , int i2 , ...int iM } be the set of first-order interactions selected by our GMDH-based methodology. To integrate them into the Bayesian hierarchical model, we modify the linear predictor η ik (see section 3.3.1) such that:

η ik = β 0k + N j=1 β jk X ijk + M m=1 β mk int im

Visualizing the effects of interactions

In our experiments, the GMDH polynomial highlights a major first-order interaction between speed limit and altitude. Triptych plots, introduced by Mc Elreath [83, p.234], are made to visualize such interactions. Thus, Fig. 3.13 depicts the bivariate relationship between annual average daily traffic (AADT) and predicted crash counts for cluster 1 from Fig. 3.5, depending on whether or not an interaction with the number of resting places is integrated into the Bayesian model. We observe that the slopes of the regression lines, constant and positive for both models, are steeper when considering an interaction. Thus, taking into account its interaction with the number of resting places, the positive influence of AADT on predicted crash counts gets more pronounced when the number of resting places increase.

Implementation details

We trained the polynomial network with GmdhPy5 , an open source Python implementation of the GMDH algorithm. For fast computation time, we restrict the maximum number of layers to 5. The number of selected best neurons is equal to the number of original features, its default value. Once the model is trained, the underlying polynomial is obtained by recursive parsing of the network, from the final node to the input layer containing the explanatory variables. Major first-order interactions are identified when the absolute value of their coefficients are superior to a threshold value defined as δ × max, with max being the coefficient with the highest amplitude, and δ = 0.01.

Model evaluation and validation

To use our model for point-estimate prediction, we must derive a Bayesian estimator from the posterior distributions. We use the mean of the posteriors which can be shown to minimize the mean squared error. In that way, we can compare our approach to black-box ML algorithms with standard quality metrics on a test dataset. Moreover, to check if the estimation of the posterior distributions converged well, we use the Posterior Predictive Check (PPC) graphical analysis method. Indeed, according to Gelman and Hill in [44, p. 158], PPC is a technique to "simulate replicated data under the fitted model and then compare these to the observed data". It allows one to look for systematic discrepancies between real and simulated data [START_REF] Gelman | Bayesian data analysis[END_REF].

To illustrate this, we compare in Fig. 3.14, for two examples of clusters, the histogram of observed crashes with that of samples drawn from the posterior distribution of crash-counts. Dashed green lines indicate the means of, respectively, the observed and replicated data. The two distributions being similar, our model fits adequately the data. Thus, the integration of latent structure and interactions do not bring skewed prior knowledge that would disturb the inference.

Conclusion

In this section, we have seen how a Bayesian learning of a hierarchical GLM can be enhanced by the introduction of a hierarchical structure and interactions identified by methods from the field of explainable artificial intelligence. In section 4.5, we report on results obtained on a 5-fold cross-validation, showing that these models improve the predictive performance of standard GLM (without multilevel formulation and interaction) and of the other simple interpretable models such as linear models or shallow decision trees. However, this first proposal underlines several limitations. First, if one does not want the user to be solely responsible for deciding which functions of the original variables are to be considered, a strategy must be adopted to explore a combinatorial space of potential functional forms. Also, for the model to stay interpretable in presence of a large number of explanatory variables, the number of relevant coefficient shouldn't be too large or ones can suffer from the information overload phenomenom [START_REF] Russell L Ackoff | Management misinformation systems[END_REF][START_REF] Keller | Effects of quality and quantity of information on decision effectiveness[END_REF][START_REF] Poursabzi-Sangdeh | Manipulating and measuring model interpretability[END_REF]. Thus, in our second proposal, we try to meet these challenges by computing a cluster-specific ranking of expansions of sparse regularized linear models ordered by increasing complexity thanks to symbolic regression.

Interpretable hierarchical symbolic regression

Throughout the previous sections, we described how a relevant hierarchical structure identified by analyzing the results of a post hoc explanation tool can enhance the predictive performance and interpretability of GLM by means of a multilevel formulation. In this section, we further exploit the discovery of this hierarchical structure by using symbolic regression to capture more complex but still interpretable relationships between the explanatory variables and the dependent variable.

First, we detail the related work of symbolic regression (SR) in section 3.4.1. Then, in section 3.4.2, we introduce SR with a focus on the technique we use to conduct the search in the space of mathematical expressions, viz., a multi-objective optimization extension of the simulated annealing algorithm. Moreover, we present in section 3.4.3 how we successfully manage to build, through a partial pooling approach, a hierarchical symbolic regression to account for the structure previously discovered. Finally, in section 3.4.4, we explain how uncertainty knowledge is provided thanks to Bayesian Inference.

Related work

Symbolic regression consists in exploring a large space of functional forms to discover a predictive model with a good trade-off between accuracy and simplicity. Each element of this space is a parametric regression or classification model whose performance is measured (e.g., with cross-validation) on a given dataset after fitting its parameters. Both the parameters and the functional form of a predictive model are learned based on available data. A wide variety of approaches have been tried to effectively explore the space of functional forms.

Evolutionary algorithms and other optimization techniques

In genetic programming, population of a mathematical expressions evolves through selection, crossover and mutation to improve a fitness function [START_REF] Mckay | Using a tree structured genetic algorithm to perform symbolic regression[END_REF][START_REF] Adriano | Symbolic regression via genetic programming[END_REF][START_REF] Schmidt | Distilling free-form natural laws from experimental data[END_REF][START_REF] Amir Haeri | Statistical genetic programming for symbolic regression[END_REF][START_REF] William | Contemporary symbolic regression methods and their relative performance[END_REF]. Other SR are based on the metaheuristic algorithm of Pareto simulated annealing to discover a set of models which are optimal in terms of a balance of both accuracy and simplicity metrics [START_REF] Stinstra | Metamodeling by symbolic regression and pareto simulated annealing[END_REF]. Thanks to use of Meijer G-functions, SR can also be approached by algorithms based on gradient descent [START_REF] Ahmed | Demystifying black-box models with symbolic metamodels[END_REF]. Bayesian processes, with algorithms based on the Markov Chain Monte Carlo strategy have been used as well to solve the SR problem [START_REF] Jin | Bayesian symbolic regression[END_REF].

Symbolic regression based on Machine Learning

Recent studies apply deep learning methods to symbolic regression in order to discover physical laws from experimental data. In [START_REF] Silviu | Ai feynman: A physics-inspired method for symbolic regression[END_REF], Udrescu et al. create a framework for symbolic regression based on neural network to discover hidden simplicity in the data (e.g., symmetry, separability) in order to decompose complex problems into simpler sub-problems. They apply their framework on 120 selected equations from classical mechanics, electromagnetism, quantum mechanics and obtain promising results on simulated data as they manage to discover more than 90% of the equations. Lately, this approach has been improved with the integration of an information complexity metric by means of Pareto optimization [START_REF] Silviu-Marian | Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity[END_REF]. However, on a recent benchmark realised by La Cava et al. [START_REF] William | Contemporary symbolic regression methods and their relative performance[END_REF], this framework obtain poor results on real data and is even dominated by simple linear models.

Recently, Petersen et al. [START_REF] Brenden K Petersen | Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients[END_REF] use a hybrid approach that combines genetic algorithms and a recurrent neural network (RNN) trained by reinforcement learning to generate better symbolic models at each iteration. Finally, Valipour et al. [START_REF] Valipour | Symbolicgpt: A generative transformer model for symbolic regression[END_REF] consider the problem as a sub task of language modelling and train a generative RNN model with reinforcement learning to produce symbolic equation skeletons whose constants are further adjusted by the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [START_REF] Fletcher | Practical methods of optimization[END_REF].

Areas of application

SR finds applications in numerous domains such as physics [START_REF] Schmidt | Distilling free-form natural laws from experimental data[END_REF], finance [START_REF] Chen | Genetic algorithms and genetic programming in computational finance[END_REF], climate modeling [START_REF] Stanislawska | Modeling global temperature changes with genetic programming[END_REF] or renewable energies [START_REF] William | Automatic identification of wind turbine models using evolutionary multiobjective optimization[END_REF]. So far, only few studies applied symbolic regression to safety analysis. Meier et al. [START_REF] Meier | Symbolic regression for precrash accident severity prediction[END_REF] use prioritized grammar enumeration, a dynamic programming version of symbolic regression, to predict crash severity a few milliseconds before collision. Patelli et al. [START_REF] Patelli | Traffic modelling and prediction via symbolic regression on road sensor data[END_REF] design a GP-based symbolic regression to predict the traffic flow. To the best of our knowledge, symbolic regression has not been applied to long term crash predictions.

Symbolic regression with Pareto simulated annealing

General description

Most versions of symbolic regression (SR) discover an expansion of a linear model with the addition of non-linear effects by searching a space of functional forms. Section 3.4.1 gave an overview of the various methods that have been used to perform this search. Among them, we select simulated annealing, an effective metaheuristic known for its robustness in optimization problems involving a large search space [START_REF] Eren | Küçükdemiral, and İlker Üstoğlu. Chapter 2 -introduction to optimization[END_REF][START_REF] Daniel Delahaye | Simulated annealing: From basics to applications[END_REF]. Thus, we represent the problem as a local search. Moreover, we adopt a multi-objective extension of the simulated annealing algorithm to perform the search while optimizing both the complexity and the accuracy of the models [START_REF] Stinstra | Metamodeling by symbolic regression and pareto simulated annealing[END_REF]. The search ends on a set of mutually non-dominated predictive models, the Pareto front.

Definition of a solution

The functional form of a model is extracted from a set of expression trees. Each expression tree is perfect, binary and consists of internal operator nodes and leaves. Leaves are either represented by a constant or an explanatory variable. Operator nodes can be unary (e.g., cos, sin, tan, exp, ln, left, right) or binary (e.g., +, ×, /) and have two children. For unary operators, we indicate with the subscripts " l " (for "left") and " r " (for "right") to which child the operation is applied. For instance, if the operator is ln l , then the logarithm is applied to the left child. The left and right operators apply the identity function to the left and right child, respectively.

We extract the symbolic expression by a breadth-first traversal of the expression tree. In practice, as operators, constants and input variables are defined with Sympy, an open-source Python library for symbolic computation [START_REF] Meurer | Sympy: symbolic computing in python[END_REF], the traversal returns a Sympy expression. Finally, the functional form S of a solution is obtained by the combination and algebraic simplification of the Sympy symbolic expressions of a set of expression trees (see Fig. 3.15).

Initialization of a first solution

Function initialize of algorithm 1 generates a first solution represented by a set M cur of random expression trees, with S cur the associated functional form. To this end, this function first creates balanced binary trees, each of the same depth. Then, for each tree, an inorder traversal associates an index to each node. Odd indices refer to internal nodes while even indices refer to the leaves (see Fig. 3.16a). In this way, we have an efficient means to search for a node in a tree and to know directly what type of node it is (see section Neighbourhood of a solution). At the same time, internal nodes are initialized with an operator chosen with equiprobability from the set of predefined operators introduced in section 3.4.2. Each leaf has a 50% probability of being initialized either to a constant or to one of the explanatory variables. In the latter case, each explanatory variable is equiprobable.

Neighbourhood of a solution

Function generate of algorithm 2 generates a new solution S new in the neighbourhood of the current solution S cur . It randomly selects an expression tree from M cur and a node index from {0, ..., 2 T -2}, T being the tree depth. Then, a recursive search finds the node with the selected index. When the node is an operator (viz., its index is odd), it is replaced by a randomly selected operator. Likewise, when a leaf is selected (viz., its index is even), it is replaced by a randomly selected constant or explanatory variable. Fig. 3.16 illustrates this process.

If unchecked, function generate could lead to ill-defined operators. For instance, a logarithm could be applied to a potentially negative domain. Therefore, function integrityCheck infers recursively the domain of each operator node and, based on rules from interval arithmetic, checks its validity (table 3.1 introduces some of these rules). With interval arithmetic, we have an efficient way to ensure that the functional form generated from the random process does not contain any undefined values [START_REF] Keijzer | Improving symbolic regression with interval arithmetic and linear scaling[END_REF]. For more details on the integrity check, we refer to [START_REF] Stinstra | Metamodeling by symbolic regression and pareto simulated annealing[END_REF].

Cost of a solution

The cost of a solution is measured in terms of both the prediction error (see function measurePerformance of algorithm 2) and the complexity of the functional form (see function measureComplexity of algorithm 2).

Performance To obtain a robust estimate of the prediction error of S new , we compute the average RMSE, for a regression task, or f 1-score, for a classification task, on

Algorithm 1 Symbolic regression with Pareto simulated annealing

Require: m: number of expression trees; T min = 0.0001: initial temperature (heating phase) and minimum temperature (cooling phase); λ h = 1.15, λ c = 0.85: ratios between two adjacent temperatures in the heating phase and cooling phase, respectively; s h = s c = 300: number of iterations between two updates of temperature; γ c = 1.15: ratio that controls the growth of s c ; max: maximum number of iterations during the cooling phase. while α ≤ 0.9 do ▷ Heating phase

T ← T × λ h , α ← acc/(acc + rej) 12:
acc ← 0, rej ← 0

13: i = 0, ζ = ∅ 14:
M curr , S curr ← initialize (m)

15:

while T > T min and i < max do ▷ Cooling phase M curr ← M new 15:

return ζ (a) x 1 x 2 + x 2 (b) x 1 + 2x 2 (c) x 1 + 2x 2 (d) x 1 + x 2 + x 0 x 2
i ← i + 1 16:
return i, ζ, M curr , S curr , acc, rej the validation subsets of a 5-fold cross-validation process. On each training subset, the coefficients β i of S new are learned by solving an l2-regularized linear regression, for regression tasks, or an l2-regularized logistic regression, for classification tasks.

In either case, the regularization parameter is determined on each training subset of the aforementioned 5-fold cross-validation either by an efficient generalized crossvalidation [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF] for ridge regression or by a nested cross-validation process for logistic regression. Once the estimate of the prediction error is obtained, the coefficients β i are fitted one last time on the whole training dataset.

Complexity We improve the strategy introduced in [START_REF] Stinstra | Metamodeling by symbolic regression and pareto simulated annealing[END_REF] to propose a new measure of the complexity of a solution. We penalize both the collinearities and the number of terms present in the symbolic expression of the functional form. The complexity of a solution S composed of m terms is defined as:

Complexity(S ) = m i=1 1 + max |r ij |; j ∈ {1, 2, ..., m} \ i C i (3.3)
where r ij is the Pearson's correlation coefficient, computed on the training dataset, between terms i and j, and C i is the complexity of the term i. We use algebraic rules to compute the complexity of each term, some of which are presented in table 3.2.

The complexity of a unary operator (e.g., the natural logarithm) is determined by approximating the operator, on its inferred domain, by a polynomial of increasing degree (at most 10) until the score of the fit, as measured on a validation set, is below a predefined threshold. The complexity of the unary operator is then defined as the degree of the best polynomial approximation. It should also be noted that, according to equation 3.3, the more terms a solution has, the more complex it is. We were able to confirm experimentally that the measured complexity represents well the complexity perceived by the safety experts.

Comparison of two solutions

The search ends with a set of Pareto optimal solutions that belong to the boundary beyond which neither the prediction error nor the complexity can be improved without deteriorating the other objective. This can be formally defined in terms of a dominance relation. Let U 1 be the prediction error and U 2 be the complexity metric.

S a dominates S b ≡ ∀i ∈ {1, 2} : U i (S a ) ≤ U i (S b ) and ∃j ∈ {1, 2} s.t. U j (S a ) < U j (S b )
Thus, the search returns a set of non-dominated solutions called the Pareto front.

Term i

Complexity C i Example Computed C i const 0 2 0 x 1 x 4 1 f (x) n n × C(f (x)) a x 2 2 2 f (x) × g(y) C(f (x)) + C(g(y)) x 1 x 2 2 3 f (g(y)) C(f (x)) × C(g(y)) ln(3x 2
2 ) C unary (ln) × 2 b a C(.) is the complexity of the inner function b C unary (.) is the complexity of the unary operator Table 3.2: Algebraic rules used to compute the complexity of each term, adapted from [119, p.320] Exploration by Pareto simulated annealing Simulated annealing (SA) is an iterative local search process used to solve optimization problems for which a simple hill-climbing approach would most often converge on a poor local optimum. At each iteration, SA generates randomly a solution S new in the neighborhood of the current solution S cur . The probability P of accepting S new as the new current solution is a function of both a temperature parameter T and the difference in cost ∆E between the two solutions.

P = e -∆E/T (3.4)
Annealing temperature T SA mimics the physical process of annealing in metallurgy where a material is first heated before being gradually cooled in order to reach an equilibrium state with increased ductility and hardness. SA follows a similar two-steps process.

The heating phase aims at discovering an initial temperature T 0 that favors exploration over exploitation in the beginning of the search. The heating process starts from a low temperature at which a deteriorating neighbour of the current solution is rarely accepted. Then, every s h iterations, the temperature is increased according to a geometric series of ratio λ h > 1. The process ends at a temperature T 0 at which at least 90% of the randomly generated neighbours are accepted.

During the cooling phase, the annealing temperature is progressively decreased, every s c iterations, according to a geometric series of ratio λ c < 1. High temperatures favor the exploration of the space of functional forms by preventing the process from converging too early on a local optimum (see Fig. 3.17). On the contrary, the more the temperature decreases, the less likely it is for a deteriorating neighbour to replace the current solution (see Fig. 3.17). The value of λ c controls the speed at which the annealing temperature decreases. If λ c is too small, the optimization may stay stuck too early in the neighborhood of a poor local optimum. Whereas, if λ c is too close to 1, the optimization may take too long to reach a good optimum. Moreover, parameter s c increases according to a geometric series of ratio γ c > 1. Thus, more iterations are allocated to lower temperatures to favor the exploitation of promising functional forms. Finally, the search ends when either the temperature falls below a threshold or the number of iterations reaches a predefined maximum.

∆E and the acceptance of a new solution For a single-objective optimization problem, ∆E is simply the difference of the objective function evaluated at two neighbouring solutions. For our multi-objective optimization problem, we use the dominance-based performance metric introduced above. When a new solution S new dominates, or is as good as, S cur , it is accepted as the new solution (see function accept of algorithm 3). When S new is less effective than S cur , it has a probability P defined by eq. 3.4 to be accepted. In that case, ∆E is defined as:

∆E(S cur , S new ) = 1 | ζ| | ζSnew | -| ζScur | ( 3

Partial pooling approach for symbolic regression

To handle correlations among groups of observations in the dataset, we developed a two-stage partial pooling approach to build hierarchical symbolic regressions. First, we start by running a symbolic regression on the whole dataset to automatically discover global models. Then, for each cluster in the hierarchical structure discovered in section 3.2, we run a new symbolic regression to extract cluster-specific functional forms that are merged afterwards to the functional form of a global model previously selected.

Automatic discovery of global models

In section 4.6, where we illustrate the dynamic interpretative process made possible by our framework, we emphasize the interest of being able to let the user choose a predictive model on the Pareto front. In that way, the end-user can precisely balance between the predictive performance and the simplicity of the model. However, in our proposed methodology, we also need a principled way to automatically select a model on the Pareto front. To do this, first, we consider the point Ω in the Pareto plan that, (i) on the performance axis, is at the level of the most efficient model encountered and, (ii) on the complexity axis, is at the level of the simplest model encountered. Euclidean distance. This model, located in the elbow of the Pareto front, is likely to offer a good trade-off between predictive performance and complexity. In the next stage of our approach, it is used as a starting point to build cluster-specific models.

Cluster-specific models

To discover cluster-specific phenomena, for each cluster discovered by the approach introduced in section 3.2, a modified version of the symbolic regression search is conducted. It consists in merging the functional form built from the expression trees with the fixed functional form of S glob (see Fig. 3.19): common terms are grouped together and new terms are added to the formula. Hence, the marginal effects already represented by S glob can be reduced or amplified and new cluster-specific effects can be discovered. It corresponds to a partial pooling approach where cluster-specific models can benefit from the effects already discovered by the global model.

In classification tasks, the dependent variable can be highly imbalanced on some clusters only. For imbalanced clusters, we apply jointly the synthetic minority oversampling technique (SMOTE) [START_REF] Nitesh V Chawla | Smote: synthetic minority over-sampling technique[END_REF] and edited nearest neighbor (ENN) [START_REF] Dennis L Wilson | Asymptotic properties of nearest neighbor rules using edited data[END_REF]. The former generates samples from the minority class with interpolation. The later applies under-sampling to clean the noisy samples generated with SMOTE. Also, it should be noted that in the case of a homogeneous cluster, a cluster-specific model would not make sense and we propose to stay with the global model S glob .

Uncertainty estimation

Our approach results in global and cluster-specific expansions of linear models. Therefore, the marginal effects of the terms composing the models are readily interpretable. However, since the training set has been used to estimate the l2-regularization hyper-parameters, there is no simple linear relationship between uncertainty in the parameters and uncertainty in the target.

Bootstrap and asymptotic statistics

Bootstrap techniques could estimate the uncertainty in the parameters. Still, a standard bootstrap approach is not appropriate since the bias introduced by the penalty term would not be correctly estimated. Double bootstrap techniques have been proposed [START_REF] Hrishikesh | Double bootstrap for shrinkage estimators[END_REF][START_REF] Mccullough | Implementing the double bootstrap[END_REF] to take into account an estimation of the bias. Nonetheless, they are computationally expensive (O(n 3 ) where n is the number of samples). Also, asymptotic statistics have been derived to measure the uncertainty in the parameters under a fixed setting of the regularization parameter [START_REF] Firinguetti | Asymptotic confidence intervals in ridge regression based on the edgeworth expansion[END_REF]. They wouldn't be appropriate in our case since we estimate the regularization parameter by leave-one-out crossvalidation.

Bayesian interpretation of Ridge regression

In our work, we make use of a well-known equivalence [START_REF] Mehta | A high-bias, low-variance introduction to machine learning for physicists[END_REF] between the ridge regression regularization parameter and the parameters of a Gaussian prior for the Bayesian formulation of linear regression. It can be shown that the variance of the zero-centered Gaussian prior τ 2 must be defined as:

τ 2 ≡
σ 2 λ where λ is the ridge regularization parameter and σ 2 is the variance of the likelihood that can be estimated by measuring the variance of the target on the training dataset.

More specifically, in the frequentist approach, the coefficient estimate βRidge of Ridge regressions is solved from:

βRidge = arg min[(Y -Xβ) ⊤ (Y -Xβ) + λβ ⊤ β] (3.6)
In the Bayesian paradigm, no single coefficient estimate is found but a posterior distribution of β is inferred from the data. From the Bayes theorem, this posterior distribution can be written as:

p(β|Y, X) ∝ p(β) • p(Y |X, β) (3.7)
where Y | X, β ∼ N (Xβ, σ 2 I) with σ > 0 is the Bayesian formulation of linear regression. If we suppose zero-centered gaussian priors with variance τ 2 , then Eq. 3.7 becomes:

p(β|Y, X) ∝ exp - 1 2 (β -0) ⊤ 1 τ 2 I(β -0) • exp - 1 2 (Y -Xβ) ⊤ 1 σ 2 (Y -Xβ) = exp - 1 2σ 2 (Y -Xβ) ⊤ (Y -Xβ) - 1 2τ 2 β ⊤ β 60
From this expression, we can compute the Maximum A Posteriori estimate βMAP :

βMAP = arg max exp - 1 2σ 2 (Y -Xβ) ⊤ (Y -Xβ) - 1 2τ 2 β ⊤ β = arg min 1 σ 2 (Y -Xβ) ⊤ (Y -Xβ) + 1 τ 2 β ⊤ β = arg min (Y -Xβ) ⊤ (Y -Xβ) + σ 2 τ 2 β ⊤ β
which is equivalent to the Ridge regression estimate found in Eq. 3.6 when λ ≡ σ 2 τ 2 . Thus, for each Pareto optimal solution, we start from the discovered functional form and the value of the ridge regularization hyper-parameter λ to infer again the coefficients, but this time, using Bayesian inference with the above prior. The resulting posterior distributions give an estimate of the parameters' uncertainty.

Implementation details

We implement our model in Python. In order to converge towards interpretable models, we restrict the operators available to the symbolic regression to {left, right, ln} for the unary ones, and {×, +, -} for the binary ones. For the algebraic simplification of the expression trees by, e.g., grouping common terms together (see function simplify in algorithm 2), we use a module6 from the Sympy library ( [START_REF] Meurer | Sympy: symbolic computing in python[END_REF]).

To fit the coefficients of a newly discovered functional form, we use the scikitlearn 7 [96] implementations of ridge regression (in the case of a regression task) and l2-regularized logistic regression (in the case of a classification task). The optimal coefficients of the linear models are computed with regularized least square algorithms. Indeed, with the introduction of a weight decay, better generalization performances can usually be achieved and the models are less prone to the negative effects of multicollinearities. We optimize the l2-regularization parameter by either an efficient form of leave-one-out cross-validation (viz., generalized cross validation) for regression tasks or by five-folds cross-validation for classification tasks. However, for even moderately large classification datasets (e.g., Adult in table 4.2), performing a k-fold cross-validation at each iteration of the symbolic regression can take a long time. In that case, we perform the search using the scikit-learn implementation of a ridge classifier to benefit from the efficient computation of the generalized cross validation. Only when a model is added to the Pareto front, do we fit its parameters by l2-regularized logistic regression with k-fold cross-validation to optimize the regularization parameter. Indeed, we observed experimentally that an l2-regularized logistic regression model tends to perform slightly better than a ridge classifier. By adopting this strategy, we reduce the computation time by a factor of 6 on the Adult dataset, while the final Pareto front remains nearly identical.

For imbalanced classification tasks, we use the imbalanced-learn implementation8 of the SMOTE-ENN resampling technique [START_REF] Nitesh V Chawla | Smote: synthetic minority over-sampling technique[END_REF][START_REF] Dennis L Wilson | Asymptotic properties of nearest neighbor rules using edited data[END_REF].

As explained in section 3.4.4, in order to endow our final models with uncertainty estimates, we use Bayesian inference to compute the posteriors for the coefficients of each functional form on the Pareto front. We apply gaussian priors corresponding to the already known optimal value of the regularization hyper-parameter (see section 3.4.4). We rely on the pymc39 library with the No U-Turn Sampler ( [START_REF] Matthew D Hoffman | The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo[END_REF]) to run simultaneously two Markov chains for 3000 iterations, with a burn-in period of 1000 iterations.

Chapter 4 Experiments

In the thesis, we design intrinsically interpretable predictive models to improve highstakes decision-making processes. In our first contribution, we propose a methodology that combines Bayesian learning of hierarchical GLM with automatic detection of a hierarchical structure and interactions through methods borrowed from the field of explainable artificial intelligence (XAI). In our second contribution, we develop a novel methodology to exploit even better the hierarchical structure by combining symbolic regression and multi-objective optimization to obtain models that are sparser, more efficient while capturing more relevant interactions.

To validate that our approaches improve the predictive capacities of interpretable models and get close to the black-box models, our experiments are carried out on the highway dataset and on more than ten public datasets covering different tasks from various domains.

In this chapter, we begin by presenting in section 4.1 the different datasets on which the experiments are carried out and we introduce the different data processing in section 4.2. Then, we introduce in section 4.3 the evaluation metrics for regressions and classifications. In section 4.4, we describe all the parametric and non-parametric models that are used for comparison purposes. Then, we share the results in section 4.5. Finally, in section 4.6, we illustrate the benefits of our approach by introducing, on a realistic case study, an application we designed for highway safety experts 63

Description of the datasets

French highway dataset

To address the problem of crash prediction, we generate a tabular dataset from the different available data sources. Most of the data come from the APRR's relational databases, and the remaining data, in particular topographical surveys of the network infrastructure, was sent to us by field experts in the form of structured flat files. To generate the dataset, we first divide the highway network into segments of equal length, then for each segment, we calculate the different variables on a given time period. In this section, we describe in more detail the process that leads to the generation of this dataset, further referred to as the French Highway dataset.

Highway network meshing

In the first place, the network must be divided into segments in order to be able to identify afterwards more or less crash-prone configurations on the network and to understand the associated influential explanatory variables. For this, we rely on spatial surveys of reference points of each highway belonging to the network. Each highway is divided into one-way segments of equal length. In this study, we select a length equal to 10 kilometers, in order to be able to capture the variations of the explanatory variables on the vast network while having a crash count rarely equal to zero. Of course, at the ends of the network, some segments will have a length less than the predefined one. Thus, we decide not to consider them in the analysis.

Calculation of variables

Considering the data available to us, we chose to restrict the analysis to the APRR network only. The study therefore covers 1,894 km of network, since the 429 km of the AREA subsidiary's highway network (see Fig. 1.1) have been withdrawn. The data covers eleven years, from January 1, 2008 to December 31, 2018. In our work, we focus on predicting the annual crash count. The available temporal explanatory variables are therefore also aggregated to this time scale. The spatial variables remain fixed over all the years of the study. In this section, we describe how we calculate the different variables, starting with the accident count. Descriptive statistics of the computed variables are given in table 4.1.

Crash count

The accident data comes from police reports or descriptive sheets filled out by APRR staff present at the crash site. The information transmitted are very detailed: location, date and time, number of individuals involved, type of vehicle, etc. About 20% of accidents have no time indication because the protagonists did not wait for the arrival of help, police or group staff before leaving. However, they are identified by the patrols on the same day and are therefore also taken into account.

Then, accident data are stored in the company's large databases. It is from these databases that we are able to retrieve all crash-related information needed for the study. For each highway segment, the accidents are aggregated over each year of the study. The number of crashes per segment varies greatly (see Figure 4.1 for the distribution of observed crash counts). In table 4.1, we observe that the variance of the crash count exceeds its mean. This phenomenon, called over-dispersion, is often witnessed with crash data [START_REF] Lord | The statistical analysis of crashfrequency data: A review and assessment of methodological alternatives[END_REF].

Traffic related variables Two traffic-related variables are considered in the study: the annual average daily traffic (AADT) and the percentage of heavy vehicles. Two data sources are available to estimate these variables. One, coming from the 797 counting loops integrated into the road pavement, represents the number of vehicles and the percentage of heavy vehicles having passed through the sensors during a fixed period of time (6 min). The other comes from transactions at toll gates, taking into account the class of vehicle (e.g., light vehicle, two-axled truck, etc.). While the former often has noisy data (e.g., outliers, long absence of data), mainly due to sensor failures, the latter provides a very reliable estimate of the AADT and percentage of heavy vehicles on highway segments. We therefore select toll data as a way to estimate the yearly traffic related variables on the network. We use the weighted mean to evaluate the traffic on our segments and to be able to capture traffic variations inside the segments:

x = n i=1 w i x i n i=1 w i (4.1)
with w i being the distance for which the variable has the value x i . 4.1. For the variables representing specific elements such as interchanges, tunnels, etc., we received from field experts a spatial statement of their coordinates. Variables that have very small variation on the network (viz., ramps, tunnels, toll barrier, structures) are binarized to indicate whether an element is present on a roadway segment or not. The others (viz., interchanges, resting places) are discrete variables referring to as a count of the element of interest on the roadway segment. The remaining variables (viz., right shoulder width, altitude and speed limit), representing continuous elements on the network, are calculated from the weighted mean defined in eq. 4.1.

Name

Other datasets

As said before, we extend the study to 13 other datasets covering regression and classification tasks from various domains. The size of these datasets varies greatly, both in terms of volumes (from 546 to 116640 instances) and dimensionalities (from 7 to 117 features). A general description of the datasets is given in table

Data preparation

For all datasets, we do not apply any process of dimensionality reduction such as feature selection or principal component analysis. The original data are transformed in different ways depending on the ML model. Models that introduce regularization terms have standardized data as input. Standardization makes the values of each continuous explanatory variable have zero-mean and a unit-variance:

x ′ = x -µ σ
where µ is the mean value of the variable, and σ its standard deviation. For neural networks, continuous explanatory variables are re-scaled in the [0, 1] interval:

x ′ =
x -min(x) max(x) -min(x) Moreover, categorical explanatory variables are one-hot encoded. More precisely, if n is the number of different categorical values, then one-hot encoding transforms the variable into n -1 binary explanatory variables describing the values of the original explanatory variable.

Predicted class Positive Negative

Actual class

Positive True positive (TP) False negative (FN) Negative False positive (FP) True negative (TN)

Regression

To measure the performance of predictive models on regression tasks, we compute the root mean square error and the mean absolute deviation. Given n the number of observations, y i the target and ŷi the predicted value:

Root mean square error is a measure of how spread out are the prediction errors, in the context of numerical predictions. This metric is defined as:

RM SE = 1 n n i=1 (ŷ i -y i ) 2
Mean Absolute Deviation is a measure of variability that indicates the average distance between the predictions and the average value of the observed variable:

M AD = 1 n n i=1 |ŷ i -ȳ|
with ȳ being the mean value of y.

Classification

To evaluate our models on binary classification tasks, we use metrics derived from the confusion matrix presented in table 4.3. Each row of the matrix represents the instances in an actual class while each column represents the instances in a predicted class. However, many machine learning and statistical models predict a probability.

To know the associated predicted class, a classification threshold is used, most of the time equal to 0.5. Probability values greater than this threshold are mapped to one class, and the remaining are mapped to another.

From the confusion matrix in table 4.3, we can compute:

accuracy is the proportion of correct results that a classifier achieved: acc = T P + T N T P + T N + F P + F N This very intuitive metric loses interest when data are imbalanced, as it can give an over-optimistic measure of the model's predictive performance. For instance, on the Breastcancer dataset where the percentage of negative cases (viz., the tumor is benign) is 0.94 (see table 4.2), a model that has learned to detect only benign tumors will obtain a high accuracy. However, measuring the model's ability to predict malignancies can be equally useful, depending on the context. Thus, other metrics must be considered.

precision is the number of correctly predicted positive observations over the total number of predicted positive observations: pre = T P T P + F P recall is the number of correctly predicted positive observations over the total number of observations in the actual class: rec = T P T P + F N f 1-score is the harmonic mean of the precision and recall:

f 1 = 2 *
pre * rec pre + rec Both the false positives and false negatives are considered which makes this metric very useful, especially when the data are imbalanced.

Receiving Operator Characteristic (ROC) curve is designed by plotting the true positive rate (TPR) against the false positive rate (FPR) at various threshold values (see Fig. 4.2a). TPR indicates how many positive predicted outcomes occur among all positive samples: T P R = T P T P + F N FPR defines how many incorrect positive outcomes occur among all negative samples: F P R = F P F P + T N ROC curve is appropriate when the observations are well balanced between each class. However, it can give an overly optimistic view of the model's performance if datasets are very imbalanced. Precision-recall curve illustrates the trade-off between precision and recall for various threshold values (see Fig. 4.2b). This curve is tailored for problems where classes are very imbalanced, but also when the positive class is more interesting than the negative one.

The two aforementioned parametric curves are not performance measures by themselves but representations of how well classifiers can discriminate between the two classes for different threshold settings. Thus, to obtain a measure of performance across all possible classification thresholds, we compute the area under the curve. These metrics are further referred to as AUROC and AUPR for the ROC curve and the precision-recall curve, respectively.

Models used for comparison

Parametric interpretable models

Among the most simple and interpretable models, we select the scikit-learn implementations of ordinary least square regression (OLS ), logistic regression (LR) and decision trees of depth no more than 5 (to preserve interpretability). We also consider a standard Bayesian inferred GLM (B-GLM ) implemented with the Pymc3 library. This model does not include the prior knowledge of a hierarchical structure nor interactions between variables. A variant with interactions discovered by the polynomial network presented in section 3.3.3 is also considered and is further referred to as B-GLM-int. On the French Highway dataset, we also train a generalized linear mixed model (GLMM ) implemented with the statsmodels1 library.

Moreover, we include two variants of Generalized Additive Models (GAM). For the first one, GAM-splines, based on a spline basis, we use the PyGAM2 implementation. For the second one, explainable boosting machine (EBM ), based on gradient boosting with bagging, we use the implementation provided by the InterpretML framework [START_REF] Nori | Interpretml: A unified framework for machine learning interpretability[END_REF].

We also compare our approach to genetic programming based symbolic regressions with the reference implementations of the gplearn3 package (SR-GP ) and GP-GOMEA 4 (SR-Gomea) [START_REF] Virgolin | Improving model-based genetic programming for symbolic regression of small expressions[END_REF], the latter being known to perform well on many real world datasets [START_REF] William | Contemporary symbolic regression methods and their relative performance[END_REF]. For both implementations, the set of operators is restricted to the one we use in our approach (viz., {+, -, ×, ln}). We also consider SR-Gomea-op, the same model with a less restricted set of operators (viz., {+, -, ×, ln, cos, sin, √ }), the same as the one used by [START_REF] William | Contemporary symbolic regression methods and their relative performance[END_REF] in their recent survey. Note that GP-GOMEA does not currently have an implementation for classification problems.

For all the aforementioned interpretable models, we apply a no pooling approach that accounts for the clusters discovered by our Hierarchical structure module (see section 3.2). This approach fits a separate model for each cluster and considers that no similarities exist between them.

Non-parametric black-box models

We select three highly flexible black-box models: (i) the scikit-learn implementation of Support Vector Machines (SVM) and (ii) Multilayer Perceptrons (MLP), and (iii) the XGBoost ( [START_REF] Chen | Xgboost: A scalable tree boosting system[END_REF]) gradient tree boosting library5

Models from our proposal

We consider the first proposal of ours namely, the Bayesian hierarchical GLM (BH-GLM-int). Recall that the latter is based on the bayesian inference of a linear hierarchical model with data-driven discovery of objective priors in the form of i) a hierarchical structure (chapter 3.2) and ii) strong first-order interactions obtained through the analysis of the structure of a trained self-adaptive polynomial network (chapter 3.3). We also include the variant without interactions between variables (BH-GLM ). For each dataset, the associated parameters and priors are given in table 4.4. 4.4: Parameters and priors for the Bayesian hierarchical models Moreover, we consider several variants of our second proposal. SR-trad and SRmax use only the global model of section 3.4.3 while HSR-trad and HSR-max use the cluster-specific models of section 3.4.3 (prefix "H" stands for hierarchical ). We also train our hierarchical symbolic regression on clusters discovered with a hierarchical agglomerating clustering applied to the training data, including the dependent variable. These models are further referred to as HSR-naive-trad and HSR-naive-max. Finally, SR-NP-trad and SR-NP-max are cluster-specific models learned with a no pooling approach, meaning that they do not include the knowledge of the global model. The suffixes trad and max are used to distinguish models selected near the elbow of the Pareto front, that should have a good trade-off (whence trad ) between complexity and predictive performance, from models of maximum complexity (whence max ).

Dataset Link a #clusters 1st level 2nd level 3rd level Regression FrenchHighway Log 4 Y ik ∼ NB(µ ik , α) b , α ∼ G(a = 0.5, b = 2.5) c β j ∼ N (µ, σ) µ ∼ N (0, 10), σ ∼ H(5) Insurance Id d 2 Y ik ∼ N (µ, 1000) β j ∼ N (µ, σ) µ ∼ N (0, 100), σ ∼ H(5) Airbnb Id 4 Y ik ∼ N (µ, 5) β j ∼ N (µ, σ) µ ∼ N (0, 10), σ ∼ H(2.5) House Id 4 Y ik ∼ N (µ, 1000) β j ∼ N (µ, σ) µ ∼ N (0, 100), σ ∼ H(5) Puma Id 3 Y ik ∼ N (µ, 10) β j ∼ N (µ, σ) µ ∼ N (0, 10), σ ∼ H(2.5) Satellite Id 3 Y ik ∼ N (µ, 20) β j ∼ N (µ, σ) µ ∼ N (0, 50), σ ∼ H(5) Wind Id 3 Y ik ∼ N (µ, 10) β j ∼ N (µ, σ) µ ∼ N (0, 10), σ ∼ H(2.5) Breast tumor Id 2 Y ik ∼ N (µ, 20) β j ∼ N (µ, σ) µ ∼ N (0, 50), σ ∼ H(5) Music Id 3 Y ik ∼ N (µ, 10) β j ∼ N (µ, σ) µ ∼ N (0, 10), σ ∼ H(2.5) Wine Id 3 Y ik ∼ N (µ, 10) β j ∼ N (µ, σ) µ ∼ N (0, 10), σ ∼ H(2.5) Toxicity Id 2 Y ik ∼ N (µ, 5) β j ∼ N (µ, σ) µ ∼ N (0, 10), σ ∼ H(2.5) Gas Id 2 Y ik ∼ N (µ, 20) β j ∼ N (µ, σ) µ ∼ N (0, 50), σ ∼ H(5) Classification Breastcancer Logit 2 Y ik ∼ B(p = exp(µik) 1+exp(µik) ) e β j ∼ N (µ, σ) µ ∼ N (0, 100), σ ∼ H(5) Adult Logit 4 Y ik ∼ B(p = exp(µik) 1+exp(µik) ) β j ∼ N (µ, σ) µ ∼ N (0,

Hyper-parameters tuning

For fair comparisons, the hyper-parameters of the models presented in section 4.4.1 and section 4.4.2 are optimized by cross-validation with grid-search. For each model, the grid of hyper-parameters' values are given in table 4.5.

Moreover, we conducted a grid-search for the SR's hyper-parameters (viz., number of expression trees, depth of an expression tree, parameters controlling the annealing temperature in simulated annealing) on two datasets, Insurance and Adult. We obtained the following results:

Expression tree depth Deeper trees lead to more complex models, mainly due to the possibility of deep compositions of functions. Motivated by finding a good com- promise between the complexity of the final models and their predictive performance, we restrict the tree depth to 4. In this way, each expression tree, being a perfect binary tree, has 8 leaves.

Number of expression trees

There is no noticeable improvement in predictive performance when the number of expression trees exceeds two-thirds of the number of features after data preparation (viz., one hot encoding).

Parameters of the simulated annealing The grid-search is also conducted on the parameters of the simulated annealing. For the heating phase, we find that λ h = 1.15 is a good value for reaching quickly a suitable starting temperature T 0 . For the cooling phase, the ratio λ c = 0.85 is a good balance between algorithmic efficiency and not falling into local optimums too quickly. The parameter γ c , which controls the growth of s c (viz., the number of iterations between two updates of temperature) is set to 1.15.

Results

We apply a 5-fold cross-validation to measure the predictive performances of the models. Results obtained on regression datasets are reported in table 4.7 and table 4.8, the latter for cluster-specific interpretable models trained with a no pooling approach. Results for classification datasets are shown in table 4.9.

Hierarchical structure module

For each dataset, the optimal number of clusters computed in the Hierarchical structure module is given in table 4.6. Moreover, to validate the ability of this module to associate an unknown sample to a cluster, a train-test split approach is applied on each training subset of the 5-fold cross-validation. For each training subset, the decision tree classifier is trained, on 80% of the data, to predict, based on the explanatory variables, the cluster to which a new observation belongs. A f 1-score is computed on the remaining 20% of each training subset. The decision tree classifier is highly accurate on all datasets (see table 4.6).

Regression datasets

BH-GLM and BH-GLM-int obtain better RMSE and MAD than the fully interpretable models (viz., OLS, decision tree, B-GLM, B-GLM-int) on most datasets. Moreover, on the French Highway dataset, BH-GLM obtains a 6.64% decrease (6.77% As expected, we observe the superior predictive performance of modern tree-based algorithms (viz., XGBoost and EBM ). However, the black-box models, SVM and MLP, do not appear to have the expected high predictive performances as BH-GLM and BH-GLM-int obtain similar if not better RMSE and MAD on most datasets. SR-GP, the symbolic regression based on genetic programming, obtains poor results and is even dominated by the fully interpretable models on all datasets. The more recent approach SR-Gomea obtains better predictive performance than SR-GP but is still dominated by HSR-trad and HSR-max. SR-Gomea-op does not highlight significant predictive gains compared to SR-Gomea on most datasets. This validates that restricting the operators makes it possible to obtain interpretable functional forms with more than satisfactory predictive performance on real world datasets.

HSR-trad, the model that, according to our second proposal, should offer a good trade-off between performance and complexity, obtains better RMSE and MAD than BH-GLM and BH-GLM-int on most datasets. As expected, HSR-max, the most complex model resulting from our approach, performs better than HSR-trad, except for the Insurance dataset where they obtain similar predictive performance. Moreover, the fact that these models have performance metrics with low standard deviations testifies to their robustness. Indeed, they are likely to discover similar solutions on similar datasets.

HSR-trad and HSR-max, the cluster-specific models, often show a clear improvement when compared to the global models SR-trad and SR-max. The partial pooling approach has a clear interest given that HSR-trad and HSR-max outperform SR-NPtrad and SR-NP-max, their no pooling variants. Also, our approach to discover a hierarchical structure is robust, efficient and obtain better results on all datasets compared to the approach that considers more naive clusters. Indeed, on all datasets, HSRtrad and HSR-max are substantially better than HSR-naive-trad and HSR-naivemax. Moreover, incorporating the data-driven discovery of a hierarchical structure not only provides better predictive performances, it also offers better interpretability by capturing cluster-specific phenomena (see section 4.6).

HSR-max and GAM-splines have similar performances on all datasets but the Insurance dataset where HSR-max discovers a significant interaction between the body mass index and being a smoker. However, the no pooling variant of GAMsplines is slightly better than HSR-max on the Insurance dataset. Finally, EBM performs well on all datasets. It obtains the best performances on the Insurance and Gas datasets and is similar, if not better, to XGBoost on the French Highway, Airbnb, and Wine datasets.

Classification datasets

For the breast cancer dataset, the search for a hierarchical structure (see section 3.2) identifies two homogeneous clusters and there is no need to learn cluster-specific models. Thus, for this dataset, we consider only the global models SR-trad and SRmax. Still for this same dataset, all models perform very well. Nonetheless, SR-max is as, or even more, efficient than LR and black-box models. Also, we observe that BH-GLM and BH-GLM-int, the Bayesian models that account for the hierarchical structure in the data, perform better than the standard Bayesian models B-GLM.

Embedding major first-order interaction also improves the predictive performance.

XGBoost and EBM achieve the best results on the adult dataset. HSR-trad obtains similar results to LR, thus suggesting that sparse models can have similar predictive performances than models that consider all explanatory variables, the latter being prone to the information overload phenomenon. HSR-max obtains similar, if not better, results than black-box models. Its performance is close to the one of GAM-splines. As observed previously for regression datasets, HSR-trad and HSRmax still perform better than BH-GLM and BH-GLM-int. Also, they are better than the global models (viz., SR-trad and SR-max ), thus confirming the benefits of incorporating a hierarchical structure in case of classification problems.

Discussion

Confirming previous studies [START_REF] Lou | Intelligible models for classification and regression[END_REF][START_REF] Caruana | Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission[END_REF], we observe that EBM, as a variant of GAM, is very efficient on both regression and classification tasks. Moreover, this model meets many of the expected criteria for interpretability enumerated in [START_REF] Barredo Arrieta | Explainable artificial intelligence (xai): Concepts, taxonomies, opportunities and challenges toward responsible ai[END_REF]. However, it also has limitations that can make it unsuitable for road safety analysis. First, different optimization strategies adopted to learn an EBM model, can lead to different interpretations of its predictions [START_REF] Chang | How interpretable and trustworthy are gams[END_REF]. However, for road safety analysis, trust in the identification of the main risk factors is required by experts when they elaborate remedial actions. Moreover, for satisfactory interpretability, it helps if a GAM has a small number of components and if each component function is relatively smooth. However, EBM, due to their reliance on boosted trees, can hardly maintain these constraints [START_REF] Rudin | Interpretable machine learning: Fundamental principles and 10 grand challenges[END_REF]. With our approach, field experts are more likely to be confident in models with cluster-specific behaviors and stable functional forms that highlight a selection of relevant factors and their interactions.

Furthermore, for the French Highway dataset, as already observed in [START_REF] Veran | Crash prediction for a french highway network with an xai-informed bayesian hierarchical model[END_REF], the best known strategy to estimate the number of crash counts is to average, for each highway network segment, the number of accidents that occurred in previous years (c.f., the Local model in table 4.7). We also observe that the predictive performances of GLMM are equivalent to those of the Local model. When looking at the model in more detail, we understand that the introduced random effects adjust the predictions of fixed effects such that they match the averages of the historical crashes observed on the segments. However, these models (viz., GLMM and Local ) do not offer much insight about the associations between crash counts and risk factors. We observed that flexible models, such as EBM, are able to approach in performance these local models by discovering quasi-identifiers of road segments. For example, an EBM discovers a complex nonlinear relationship between the altitude and the number of accidents, see Fig. 4.3. Accidents appear more likely for the lowest altitudes. However, this phenomenon should not be interpreted as a potential risk factor linked to the altitude.

In fact, the model is using the altitude as a proxy variable to identify a group of nearby road segments. Therefore, in that particular context, EBM, despite its good predictive performance, does not always provide relevant information to field experts. It can even, at times, mislead them.

Although plots like the one of Fig. 4.3 make it possible to identify these potentially misleading models' behaviors, the EBM model does not provide alternative associations between the explanatory variables and the target. With our approach, the risk of misinterpretation is reduced thanks to the successive models on the Pareto front: from less complex, which capture only overall effects, to most complex, which are flexible enough to focus on hazardous configurations specific to a few roadway segments. Through such a dynamic interpretative process, field experts can use the model best suited to meet their needs. In the next section, we illustrate this process with an experimental study applied to the French Highway dataset.

Dynamic interpretative process

Introduction

A glass-box interpretable crash prediction model able to highlight the marginal effects of the potential risk factors is a valuable tool to help design effective highway safety policies. Indeed, in such a context, the predictive model must limit as much as possible the risks of misinterpretation. In this section, we show how, on a realistic use case validated by field experts, the framework we propose responds to this challenge by making possible an efficient dynamic interpretative process that leads to selecting a predictive model well suited to the task at hand.

Each functional block of our framework (see figure 3.2) plays an important role in the methodology we propose. The hierarchical structure module brings out a relevant partition of the highway network. This partitioning of the measured crash counts, with their associated explanatory variables, not only improves the performance of the models (see section 4.5) but also refines the analysis of marginal effects by capturing cluster-specific phenomena. Moreover, based on the multi-objective optimization leading to global and cluster-specific models, safety experts can explore a list of optimal models that offers a variety of alternatives along the performance/complexity trade-off axis. They will often start with the simpler ones that capture only the global effects, and move towards the more complex ones that can represent more localized phenomena. To support this methodology, we developed, in dialogue with field experts, a graphical user interface (see Fig 4.4).

In the following description of a case study, we suppose that the partitioning of the highway network and the training of the global and cluster-specific models have already been done on ten years of data, from January 1st, 2008 to December 31th, 2017. Data from 2018 is used to validate that, based on out-of-sample predictions, the framework provides useful information to safety experts.

Use case

Partitioning the highway network

From the hierarchical structure module of section 3.2, safety experts identified four relevant clusters (see Fig. 4.5). In table 4.10, we share some descriptive statistics of explanatory variables and crash count for each cluster. We observe that clusters are heterogeneous on most of the explanatory variables and on the crash count. For instance, cluster 1 is representative of hazardous segments with high traffic, whereas cluster 3 is made of less accidental secondary segments that connect the main highways. where the additional variables x 0 , x 1 and x 8 are, respectively, the speed limit, the number of interchanges and the presence of tunnels. Out-of-sample predictions from models 1 and 2 differ locally (see Fig 4.7a and Fig 4.7b). In particular, in mountainous areas, segments considered as moderately hazardous by the first model, are now associated with a high risk of accidents due to the discovery of a first-order interaction between the presence of tunnel and the speed limit. Safety experts, by combining prior knowledge of the network with the observed transition from model 1 to model 2, are confident that this interaction is one of the main reasons why a large number of accidents have occurred on these segments during the ten years covered by the training data. This discovery may support a proposal for reducing the authorized speed limit on these specific segments.

Also, more refined models are sometimes able to better identify less hazardous 86 configurations. For instance, a segment close to Belfort was identified as highly hazardous by model 1 (between 20 and 25 predicted crashes) while model 2 considers it as moderately hazardous (between 15 and 20 predicted crashes). As remedial or preventive actions may involve significant costs, predicting less hazardous configurations is therefore beneficial because actions will be implemented at first on more prominent hazardous segments. So, thanks to this more complex model, field experts optimise the order of safety policies and funding shall be better allocated. Nonetheless, one of the major difficulties for the interpretation of more complex models is related to the introduction of collinearities and interactions between continuous variables. To illustrate this, consider model 3 from Fig. 4.6: Model 3: ŷ = -33.3 -0.00489x 10 + 1.57x 2 -9.37 • 10 -6 x 3 x 6 + 0.00286x 3 + 0.15x 6 Model 3 is characterized by an interaction between the averaged altitude x 6 and the traffic x 3 . By focusing on this interaction, model 3 better fits training data than model 2 but its effects plots (see Fig 4.11a) are arguably more difficult to interpret. However, our framework only produces differentiable closed-form expressions for which it is always possible to compute the partial analytical derivatives (PD) w.r.t. variables of interest, to quantify explicitly their partial effects (i.e., a measure of the conditional effect of a variable on the target) [START_REF] Seidyo | Measuring feature importance of symbolic regression models using partial effects[END_REF]. In this sense, we can understand how a unit change in an explanatory variable affects the crash count when other variables are held constant. For instance, the partial derivatives for the trafic x Histograms of the pointwise partial derivatives can be useful interpretative tools (see Fig. 4.10). Although the partial effects of the traffic x 3 are mostly positive, a few are negative for segments of high altitudes and above average traffic (see Table 4.11). Thus, these variations in partial derivatives emphasize that the relation between the crash count and x 3 is more complex than the linear dependency proposed by model 1 and model 2. By introducing this novel interaction, model 3 manages to capture more variability in the dependent variable than the previous models. Finally, model 4 is much more complex:

Model 4: ŷ = -34.7 + 8.5 • 10 -6 x 2 1 -0.000671x 1 x 3 -8.5 • 10 -6 x 1 x 6 + 13.5x 1 + 0.0113x 10 + 1.92x 2 -8.5 • 10 -6 x 3 x 6 + 0.000679x 3 x 8 + 0.00269x 3 -0.0276x 4 + 1.61x 5 + 0.133x 6 + 0.12x 7 -0.276x 8 + 0.0472x 9

As we can see from Fig. 4.11b, the introduction of new correlated terms improves slightly the fit. To capture extra variability in the dependent variable, model 4 where many accurate-but-different models exist to describe the same data [START_REF] Semenova | A study in rashomon curves and volumes: A new perspective on generalization and model simplicity in machine learning[END_REF]. As discussed by [START_REF] Rudin | Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead[END_REF], we argue that the availability of multiple efficient predictive models is useful since field experts may have more flexibility in choosing a model that they find interpretable. Moreover, we help them in this process as our definition of the complexity warns them when models are likely to be difficult to understand.

Finally, the dynamic interpretative process can be a useful tool to construct new handmade predictive models, based on the knowledge learnt by analyzing the Pareto optimal models. For instance, we have seen that the interactions introduced in model 2 and model 3 are both valuable. The user could consider building a new model with both of them.

Towards causality

A central question remains: among these different models, how can be distinguished the trustworthy ones from the ones based on spurious associations due to inductive bias? As illustrated above, one way is to rely on the diligence of the user equipped with expert knowledge and effective tools. This could also be partially automated when prior knowledge of the conditional independences between variables is formalized, e.g., as a causal graph [START_REF] Pearl | Causality[END_REF]. Such approaches are beyond the scope of our current work. However, our framework fosters a dynamic interpretative process that, combined with a clear quantification of uncertainty, is a useful tool to identify variables of interest and to understand how they interact. Therefore, we can surmise that our framework could facilitate the development of causal models.

Chapter 5 Conclusion and perspectives

Conclusion

Predictive models are being used increasingly to make high stake decisions. For many applications, there is a need for both accuracy and interpretability. For instance, in highway safety analysis, we argue that a preference should be given to predictive models that are both accurate and glass-box interpretable in order to increase the confidence of safety experts in the identification of hazardous segments.

Motivated by these requirements, in our first contribution, we introduced a framework to build efficient and interpretable Bayesian hierarchical models for regression or classification tasks. We proposed a data-driven discovery of objective priors in the form of a hierarchical structure and strong first-order interactions between explanatory variables. We start with a trained and usually efficient, albeit opaque, ML algorithm in order to compute for each observation the Shapley values of the explanatory variables. Then, a partition of the instances, related to how the ML algorithm predicts the target, emerges from the hierarchical agglomerative clustering of the observations described by the Shapley values. Furthermore, we analyze the structure of a trained self-adaptive polynomial network to discover important firstorder interactions. This prior knowledge is then embedded into the definition of a Bayesian hierarchical GLM with nonlinear functional form.

In our second contribution, we proposed to exploit even better the discovery of the hierarchical structure by using symbolic regression to discover sparse interpretable models with rich interactions. We combined a multi-objective simulated-annealingbased symbolic regression and a partial pooling approach to discover models that capture global effects and cluster-specific effects. More specifically, we start by computing a Pareto front of global predictive models. We select among these models the one offering a good trade-off between its predictive performance and its complexity.

Afterwards, for each cluster of the previously discovered hierarchical structure, the global model is used as the starting seed for a new multi-objective symbolic regression. Finally, the best models, i.e. the ones appearing on the Pareto fronts, are re-estimated through Bayesian inference in order to associate uncertainty estimates to their coefficients.

On fourteen datasets, covering both regression and classification tasks, our two proposals outperform most interpretable models. On some datasets, we achieve performance comparable to that of non-parametric black-box models. Furthermore, we presented a case study based on the highway network dataset to validate the new dynamic interpretative process made possible by our second proposal. As our approach discovers transparent and parsimonious symbolic models, safety experts can be more confident in their understanding of the relations between the explanatory variables and the dependent variable. Moreover, thanks to Bayesian inference, the risk factors are associated with measures of uncertainty. In addition, the use of Pareto optimization allows field experts to build a multi-scale view of the risk factors, from the most general to the most specific.

Perspectives

First, we plan to improve the dynamic interpretative process by generating humanreadable explanations of the relationships between Pareto optimal models, such as: model X is an extension of model Y and improves the predictive performance on training data by 3% thanks to the introduction of a new interaction between variables x A and x B . In this way, field experts will have a direct understanding of the Pareto front.

Then, we are going to test several divisions of the road network into segments. Currently, we generate the dataset for crash prediction, we first divide the network with segments of equal length and then compute the variables for each segment (see section 4.1.1). With this procedure, we can capture the heterogeneity of the explanatory and dependent variables along the network. However, other types of meshes can be used. For instance, we can use the segments identified by the SURE approach (see section 1.1). Also, Deublein et al. [START_REF] Deublein | A bayesian network model to predict accidents on swiss highways[END_REF] divide the network into homogeneous segments where explanatory variables are constant, assuming that the dependent variable is weighted by an exposure term that depends on the length of the section and traffic. Thus, future work will analyze the effects of different types of meshes. Also, the user interface will allow to choose the type of mesh. Thus, field experts will be more aware of the effects of this important prior.

Moreover, our framework relies on a specific approach to discover a hierarchical structure. Even though we validated its robustness on numerous datasets, promis-ing next steps involve analyzing other methods that are compatible with ours. For instance, in [START_REF] Rengasamy | Towards a more reliable interpretation of machine learning outputs for safety-critical systems using feature importance fusion[END_REF][START_REF] Rengasamy | Mechanistic interpretation of machine learning inference: A fuzzy feature importance fusion approach[END_REF], the authors propose an efficient ensemble feature importance method where multiple feature importance approaches are applied to a set of ML models and their crisp importance values are combined to produce a final importance for each feature. Thus, we will constitute a benchmark of feature importance methods [START_REF] Barredo Arrieta | Explainable artificial intelligence (xai): Concepts, taxonomies, opportunities and challenges toward responsible ai[END_REF] and evaluate them based on their efficiency, scalability, and on the quality of computed clusters.

Lastly, we will also extend the framework for near real-time crash risk assessment. New spatio-temporal explanatory variables, such as weather-related variables (e.g. rainfall, snow, wind, etc.) or the pavement quality, will be introduced in the analysis. In this context, since remedial actions will probably affect humans' lives even more directly, having both efficient and interpretable models will be all the more important to assist safety experts in their work. Les modèles paramétriques simples ont tendance à être plus interprétables mais moins performants que les modèles non paramétriques très flexibles qui fonctionnent comme des boîtes noires. Lorsqu'ils réfléchissent à des décisions à fort enjeu, comme dans le contexte de la sécurité routière, les experts métier s'attendent à ce que les modèles prédictifs soient à la fois performants et interprétables. En effet, les modèles doivent les aider à concevoir et à déployer des actions de sécurité préventives ou correctives. Dans ces travaux, nous contribuons à améliorer les performances prédictives des modèles paramétriques tout en conservant leur interprétabilité. En premier lieu, une structure hiérarchique bien choisie peut gérer les corrélations entre groupes d'observations et améliorer significativement la qualité des prédictions des modèles et leur interprétation. Nous proposons de l'apprendre en exploitant le résultat d'une méthode d'interprétabilité post-hoc (viz., SHAP) appliquée à un modèle boîte noire flexible (viz., XGBoost). Puis, nous introduisons deux approches algorithmiques (viz., un réseau de neurones polynomial, et une extension de la régression symbolique multi-objectif) pour découvrir des transformations non linéaires des variables explicatives. Ces dernières, tout en restant simple (e.g. interactions de premier ordre, logarithme), permettent aux modèles de capturer une plus grande partie de la variabilité dans la variable dépendante. De plus, avec la régression symbolique multiobjectif, nous calculons un classement spécifique à chaque cluster des expansions de modèles linéaires régularisés ordonnés par complexité croissante. Cela facilite un processus d'interprétation dynamique qui permet de découvrir des modèles prédictifs efficaces, efficients et interprétables. Des expériences ont été menées sur un jeu de données de sécurité routière et sur plus de dix jeux de données publics couvrant des problèmes de classification et de régression variés. Les résultats obtenus sont prometteurs étant donné que nos deux contributions surpassent les modèles interprétables traditionnels et se rapprochent des meilleurs modèles non paramétriques boîtes noires. Enfin, nous illustrons les bénéfices de notre approche en présentant, sur une étude réelle de cas, une application que nous avons conçue pour les experts de la sécurité routière.

MOTS-CLÉS

4 Experiments 4 . 1

 441 Description of the datasets . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 French highway dataset . . . . . . . . . . . . . . . . . . . . . . 4.1.2 Other datasets . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Data preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Performance metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Models used for comparison . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Parametric interpretable models . . . . . . . . . . . . . . . . . 4.4.2 Non-parametric black-box models . . . . . . . . . . . . . . . . 4.4.3 Models from our proposal . . . . . . . . . . . . . . . . . . . . 4.4.4 Hyper-parameters tuning . . . . . . . . . . . . . . . . . . . . . 4.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.1 Hierarchical structure module . . . . . . . . . . . . . . . . . . 4.5.2 Regression datasets . . . . . . . . . . . . . . . . . . . . . . . . 4.5.3 Classification datasets . . . . . . . . . . . . . . . . . . . . . . 4.5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Dynamic interpretative process . . . . . . . . . . . . . . . . . . . . . 4.6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6.2 Use case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi List of Figures

4 . 4

 44 Parameters and priors for the Bayesian hierarchical models . . . . . . 4.5 Hyper-parameters tuning . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Performance of the prediction to associate a new observation to its cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Results on the regression datasets . . . . . . . . . . . . . . . . . . . . 4.8 Results obtained by cluster-specific interpretable models . . . . . . . 4.9 Results on the classification datasets . . . . . . . . . . . . . . . . . . 4.10 Descriptive statistics of variables for each cluster . . . . . . . . . . . . 4.11 Description of explanatory variables for the overall cluster-specific data and for samples where partial derivatives w.r.t. the traffic variable are the lowest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: APRR's highway network[START_REF] Aprr | Cartography of the highway network[END_REF]

Figure 1 . 2 :

 12 Figure 1.2: Bias and variance contributing to total error[START_REF] Fortmann-Roe | Understanding the bias-variance tradeoff[END_REF]

Figure 2 . 1 :

 21 Figure 2.1: Example of a shallow decision tree trained on the French highway dataset

Figure 2 . 2 :

 22 Figure 2.2: Partial dependence of the annual average daily traffic for a gradient boosting tree model on the highway dataset

Figure 2 . 3 :

 23 Figure 2.3: Toy example of a linear SVM, with the red line being the optimal hyperplane.

Figure 2 . 4 :

 24 Figure 2.4: Multilayer perceptron (a) simple architecture with one single hidden layer (b) toy example of the computation of an output

  [START_REF] Huang | Predicting crash frequency using an optimised radial basis function neural network model[END_REF] compare a radial basis functions neural network (RBFNN) with BPNN and NB regression for crash frequencies prediction. RBFNN obtains the best results. Xie et al.[START_REF] Xie | Predicting motor vehicle collisions using bayesian neural network models: An empirical analysis[END_REF] use a bayesian neural network for crash-count prediction. The neural network model outperforms BPNN and NB regression.

Figure 2 .

 2 Figure 2.5: Sensitivity analysis performed on a Bayesian neural network for the traffic variable, from Xie et al. [133, p.930]

Figure 3 . 1 :

 31 Figure 3.1: Description of a Bayesian hierarchical generalized linear model

Figure 3 . 2 :

 32 Figure 3.2: Proposed framework for an interpretable hierarchical symbolic regression

Figure 3 . 3 :

 33 Figure 3.3: Force plot of the explanatory variables' contributions to the estimated crash count for a specific observation. On this example, the traffic has the biggest positive contribution and explains most of the crash count shift from its overall expected value

Figure 3 . 4 :

 34 Figure 3.4: Observations grouped by similarity of their forceplots

Figure 3 .

 3 Figure 3.5: Clusters automatically discovered by the hierarchical structure module, for 2018

  3.8b) between Lyon and Dijon and in the vicinity of Paris and Belfort.

  (a) Observed crash count (b) Annual average daily traffic

Figure 3 . 6 :Figure 3 . 7 :

 3637 Figure 3.6: Variables' distributions for 4 identified clusters

  (a) Observed crash count (b) Annual average daily traffic

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: Variables' distributions for 2 identified clusters

Figure 3 . 10 :

 310 Figure 3.10: Variables' distributions for 6 identified clusters

Figure 3 . 11 :

 311 Figure 3.11: Two examples of posterior distributions

Figure 3 . 12 :

 312 Figure 3.12: Structure of the GMDH model

  (a) Without interaction effects (b) With interaction effects

Figure 3 . 13 :

 313 Figure 3.13: Triptych plots of predicted crash counts vs. annual average daily traffic. Note that explanatory variables are standardized (see section 4.2).

Figure 3 . 14 :

 314 Figure 3.14: Posterior Predictive Checks (PPC) on two clusters. Left column: observed data (cluster 0: 1491 samples; cluster 2: 965 samples). Right column: replicated data (2000 simulations)

Figure 3 .

 3 Figure 3.15: A functional form associated with a set of expression trees.

1 :M

 1 function Simulated Annealing(T min , λ h , s h , λ c , s c , γ c , max) curr , S curr ← initialize (m) 8:

  , M curr , S curr , acc, rej ← explore(T, i, ζ, M curr , S curr , acc, rej) 17: if i mod s c = 0 then 18: T ← T × λ c , s c ← s c × γ c 19:

Figure 3 .

 3 Figure 3.16: A sequence of transformations applied to an expression tree. (a) An initial expression tree. Blue integers refer to node indices. (b) A transformation is applied to operator node 1, thus modifying the underlying functional form. (c) The transformation applied to leaf node 6 is muted due to its left operator parent. (d) Later in the process, node 6 can be reactivated when its parent is transformed.

Figure 3 . 17 :

 317 Figure 3.17: Effect of the annealing temperature on the acceptance probability.

  .5) with ζ the set of solutions that approximate the Pareto front, | ζ| the cardinality of ζ ∪ {S cur , S new }, and | ζS | the number of solutions in | ζ| that dominate S (see Fig. 3.18). Moreover, to smooth the estimated acceptance probability distribution, new artificial points are added to the attainment surface to get an evenly spread attainment surface over the two dimensions of the Pareto front [115]. Updates of the Pareto front Finally, when S new is accepted, the Pareto front ζ is updated (see function update of algorithm 3) by removing the solutions dominated by S new and then adding S new to ζ when it is not dominated by any other solution

Figure 3 . 18 :

 318 Figure 3.18: Example of an approximated Pareto front and its attainment surface, adapted from [119, p. 322]. From Eq. 3.5, ∆E(S cur , S new ) = (2 -4)/9 = -2/9

Figure 3 . 19 :

 319 Figure 3.19: Extraction of a cluster-specific functional form from a set of expression trees and the fixed functional form of the global model

Figure 4 . 1 :

 41 Figure 4.1: Distribution of the observed crash counts

Figure 4 . 2 :

 42 Figure 4.2: Examples of (a) ROC curve and (b) precision-recall curve obtained on the FrenchHighway dataset

  100), σ ∼ H(5) a Link function; b Negative Binomial distribution; c Gamma distribution; d Identity function; e Bernoulli distribution Table

Figure 4 . 3 :

 43 Figure 4.3: Global explanation plot provided by the InterpretML framework for an EBM model on the altitude variable on the French Highway dataset

Figure 4 . 4 :

 44 Figure 4.4: Screenshot of the web application

10 :Figure 4 . 6 : 85 (

 104685 Figure 4.5: Clusters identified for 2018

Figure 4 . 7 :

 47 Figure 4.7: Crash count predictions for 2018

Figure 4 . 8 :

 48 Figure 4.8: Effects plots and posteriors of the most influencing variables for model 1. Effects plots are obtained by computing for all observations the effect of a variable j on the crash count, defined by effect(i) j = β j x (i) j, where β j is the coefficient estimate of the j-th variable of the model and x (i) j is the value of variable j on the i-th observation[START_REF] Molnar | Interpretable Machine Learning[END_REF].

Figure 4 . 9 :

 49 Figure 4.9: Effects plots and posteriors plots of the most influencing variables for model 2

3 and altitude x 6

 6 are respectively: PD(x 3 ) = δ ŷ δx 3 = 0.00286 -9.37 • 10 -6 x 6 , PD(x 6 ) = δ ŷ δx 6 = 0.15 -9.37 • 10 -6 x 3

Figure 4 .

 4 Figure 4.11: Effects plots

  'INSA LYON, MEMBRE DE L'UNIVERSITE DE LYON NOM : VERAN DATE de SOUTENANCE : 04/11/2022 Prénoms : Thomas, Marcel, Simon, Nicolas TITRE : Efficient and interpretable crash prediction models : an application to a French highway network NATURE : Doctorat Numéro d'ordre : 2022ISAL0096 Ecole doctorale : INFOMATHS Spécialité : Informatique RESUME : Dans le monde entier, les accidents de la route ont des impacts sociaux et financiers importants. Pour réduire leur fréquence et leur gravité, les modèles de prédiction d'accidents (CPM) sont utilisés pour identifier les segments de route dangereux et fournir des indices exploitables sur les facteurs de risque associés. Les CPM sont soit des modèles statistiques paramétriques, en particulier des modèles linéaires généralisés (GLM), soit des modèles d'apprentissage automatique avec un nombre important de paramètres sans estimation d'incertitude associée (e.g., ensemble d'arbres de décision, machine à vecteurs de support…).

  : Machine Learning, model interpretability, SHAP, Bayesian hierarchical modeling, symbolic regression, multiobjective optimization Laboratoire (s) de recherche : LIRIS Directeur de thèse : Prof. Dr. Jean-Marc PETIT Président de jury : Prof. Dr. Florence SEDES Composition du jury : Prof. Dr. Pierre GANCARSKI, Prof. Dr. Gabriele GIANINI, Prof. Dr. Florence SEDES, Prof. Dr. Julien JACQUES, Prof. Dr. Jean-Marc PETIT, Dr. Pierre-Edouard PORTIER

  3.2.2 Training of a non-parametric ML model . . . . . . . . . . . . 3.2.3 Local post hoc explanations . . . . . . . . . . . . . . . . . . . 3.2.4 Discovery of a hierarchical structure . . . . . . . . . . . . . . . 3.2.5 Representation of the hierarchical structure in highway safety 3.3 Bayesian hierarchical generalized linear model . . . . . . . . . . . . . 3.3.1 Model description . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Bayesian inference . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Supervised learning of interactions with a polynomial network 3.3.4 Model evaluation and validation . . . . . . . . . . . . . . . . . 3.3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Interpretable hierarchical symbolic regression . . . . . . . . . . . . . . 3.4.1 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Symbolic regression with Pareto simulated annealing . . . . . 3.4.3 Partial pooling approach for symbolic regression . . . . . . . .

3.4.4 Uncertainty estimation . . . . . . . . . . . . . . . . . . . . . . 3.4.5 Implementation details . . . . . . . . . . . . . . . . . . . . . .

  1.1 APRR's highway network[START_REF] Aprr | Cartography of the highway network[END_REF] . . . . . . . . . . . . . . . . . . . . . . . 1.2 Bias and variance contributing to total error[START_REF] Fortmann-Roe | Understanding the bias-variance tradeoff[END_REF] . . . . . . . . . . . . 2.1 Example of a shallow decision tree . . . . . . . . . . . . . . . . . . . . 2.2 Partial dependence on the French highway dataset . . . . . . . . . . . 2.3 Toy example of a linear SVM . . . . . . . . . . . . . . . . . . . . . . 2.4 Multilayer perceptron and the computation of a neuron output . . . . 2.5 Sensitivity analysis performed on a Bayesian neural network . . . . . .15 A functional form associated with a set of expression trees. . . . . . . 3.16 A sequence of transformations applied to an expression tree . . . . . 3.17 Effect of the annealing temperature on the acceptance probability. . . 3.18 Example of an approximated Pareto front and its attainment surface 3.19 Extraction of a cluster-specific functional form . . . . . . . . . . . . . 4.1 Distribution of the observed crash counts . . . . . . . . . . . . . . . . 4.2 Examples of ROC curve and precision-recall curve . . . . . . . . . . . xiii 4.3 Global explanation plot provided by the InterpretML framework . . . 4.4 Screenshot of the web application . . . . . . . . . . . . . . . . . . . . 4.5 Clusters identified for 2018 . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Pareto front for cluster 0 specific models . . . . . . . . . . . . . . . . 4.7 Crash count predictions for 2018 . . . . . . . . . . . . . . . . . . . . . 4.8 Effects plots and posteriors of the most influencing variables for model 1 4.9 Effects plots and posteriors of the most influencing variables for model 2 4.10 Histograms of partial derivatives . . . . . . . . . . . . . . . . . . . . . 4.11 Effects plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .1 Illustrative description of the french highway dataset . . . . . . . . . 3.1 Rules for interval arithmetic . . . . . . . . . . . . . . . . . . . . . . . 3.2 Algebraic rules used to compute the complexity of each term . . . . .

3.1 Description of a Bayesian hierarchical generalized linear model . . . . 3.2 Proposed framework for an interpretable hierarchical symbolic regression 3.3 Example of a forceplot for a specific observation . . . . . . . . . . . . 3.4 Observations grouped by similarity of their forceplots . . . . . . . . . 3.5 Clusters automatically discovered by the hierarchical structure module 3.6 Variables' distributions for 4 identified clusters . . . . . . . . . . . . . 3.7 Clusters when the number of clusters is reduced to 2 . . . . . . . . . 3.8 Variables' distributions for 2 identified clusters . . . . . . . . . . . . . 3.9 Clusters when the number of clusters is increased to 6 . . . . . . . . . 3.10 Variables' distributions for 6 identified clusters . . . . . . . . . . . . . 3.11 Two examples of posterior distributions . . . . . . . . . . . . . . . . . 3.12 Structure of the polynomial network . . . . . . . . . . . . . . . . . . 3.13 Triptych plots of predicted crash counts vs. annual average daily traffic 3.14 Posterior Predictive Checks on two clusters . . . . . . . . . . . . . . . 3xiv List of Tables

14.1 Description of the French Highway dataset . . . . . . . . . . . . . . . 4.2 Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Confusion matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1 .

 1 1: Illustrative description of the french highway dataset

	Segment	Year Crash count AADT a %HV b #Interchanges . . . Presence of tollgates
	A6-2-354-364 c 2008 20	26456	19.59 0	. . . 1
	A6-2-354-364 2009 20	26343	16.9	0	. . . 1
	A6-2-354-364 2010 27 . . . . . . . . .	26656 . . .	17.31 0 . . . . . .	. . . 1 . . . . . .
	A40-1-102-112 2017 10	12701	14	0	. . . 0
	A40-1-102-112 2018 7	13015	14	0	. . . 0
	a annual average daily traffic; b percentage of heavy vehicles
	c highway name -direction -reference point begin -reference point end

Table 3 .

 3 1: Rules for interval arithmetic, from [119, p.318]. We suppose that an operator node has two children. The left one is defined on [a, b] and the right one on [c, d]. if accept(perf new , compl new , perf curr , compl curr , ζ, T ) then M curr , S curr , perf curr , compl curr ← update( ζ, M new , S new , perf new , compl new )

	Algorithm 2 Pseudo code of function explore
	3:	if integrityCheck(M new ) then
	4:	S new ← simplify(M new )
	5:	if S new ̸ = S curr then
	6:	perf new ← measurePerformance(S new )
	7:	compl new ← measureComplexity (S new )
	8:
	9: ζ, 10: acc ← acc + 1
	11:	else
	12:	rej ← rej + 1
	13:	else
	14:

1: function explore(T, i, ζ, M curr , S curr , acc, rej) 2: M new ← generate(M curr )

:

  Then, we select the model S glob on the Pareto front closest to Ω in the sense of the 57 Algorithm 3 Pseudocode of accept and update functions 1: function accept(perf new , compl new , perf curr , compl curr , ζ, T ) if perf new ≤ perf curr and compl new ≤ compl curr then function update(ζ, M new , S new , perf new , compl new ) M new , S new , perf new , compl new

	2:	is_accepted ← False
	3:
	4:	is_accepted ← True	▷ new solution dominates, or is as good as, the
		current one
	5:	else
	6:	compute P according to Eq. 3.4
	7:	draw randomly j in [0, 1]
	8:	if P ≥ j then
	9:	is_accepted ← True
	10:	return is_accepted
	1112:	is_dominated = False
	13:	for solution S in the Pareto front ζ do
	14:	if S dominates S new then
	15:	is_dominated = True
	16:	if is_dominated = False then
	17:	remove solutions in ζ dominated by S new
	18:	add S new to ζ
	19:	return ζ,

Table 4 . 1 :

 41 Description of the French Highway datasetSpatial variables A list of spatial variables is given in table

		type	Min Max Mean Std.
		Dependent variable
	Crash count	discrete	0	78	11.43	8.06
		Explanatory variables
	Temporal
	AADT	continuous 659	37644 12828 6791
	Percentage of heavy vehicles continuous 1	35.1	17.39	5.80
	Spatial
	Speed limit	continuous 104	130	128.37 4.91
	Number of interchanges	discrete	0	7	0.2	0.68
	Number of resting places	discrete	0	3	0.63	0.57
	Right shoulder width	continuous 0.6	3	2.93	0.25
	Altitude	continuous 71.97 615.8 250.83 105.71
	Presence of ramps	binary	0	1	0.24	0.43
	Presence of tunnels	binary	0	1	0.02	0.13
	Presence of tollgates	binary	0	1	0.43	0.5
	Presence of bridges	binary	0	1	0.18	0.39
	Number of instances in the dataset: 4152

Table 4 .

 4 4.2. Datasets taken from https://epistasislab.github.io/pmlb/index.html † Datasets taken from https://archive.ics.uci.edu/ml/index.php

	66

* 2: Datasets

Table 4 .

 4

	3: Confusion matrix

Table 4

 4

	.5: Hyper-parameters tuning

Table 4 .

 4 6: For each dataset: number of clusters selected and performance of the prediction to associate a new observation to its cluster for BH-GLM-int) in RMSE when compared to B-GLM, the standard Bayesian model.

	Dataset	#clusters f1 (std)
	Regression
	French Highway 4	0.995 (0.003)
	Insurance	2	1.0 (0.0)
	House	4	0.908 (0.026)
	Puma	3	0.975 (0.013)
	Satellite	3	0.964 (0.003)
	Wind	3	0.944 (0.022)
	Breast tumor	2	0.995 (0.004)
	Music	3	0.96 (0.029)
	Wine	3	0.957 (0.012)
	Toxicity	2	0.95 (0.039)
	Gas	2	0.99 (0.003)
	Airbnb	4	0.985 (0.005)
	Classification
	Breastcancer	2	0.961 (0.028)
	Adult	4	0.995 (0.002)

Table 4 .

 4 number of trees for SR-* and HSR-* models b averages and standard deviations of performance metrics obtained on 5-fold cross-validation 7: Results obtained on 12 regression datasets

		French Highway (8 a )	Insurance (6)	Airbnb (12)	Puma (6)
		RMSE (std) b	MAD (std)	RMSE (std)	MAD (std)	RMSE (std)	MAD (std)	RMSE (std)	MAD (std)
	Local	5.052 (0.250)	3.677 (0.120)	-	-	-	-	-	-
	GLMM	5.044 (0.287) 3.676 (0.140) -	-	-	-	-	-
	OLS	6.213 (0.496)	4.458 (0.216)	6077 (287)	4203 (144)	0.50 (0.011)	0.360 (0.004)	4.471 (0.068)	3.649 (0.026)
	Decision tree	6.134 (0.422)	4.405 (0.196)	4739 (324)	2738 (125)	0.490 (0.011)	0.352 (0.004)	3.688 (0.05)	2.881 (0.011)
	GAM-splines	5.80 (0.344)	4.245 (0.183)	6021 (299)	4229 (158)	0.464 (0.011)	0.330 (0.004)	4.236 (0.065)	3.492 (0.03)
	EBM	5.363 (0.233)	3.968 (0.134)	4533 (339)	2528 (131)	0.450 (0.013)	0.321 (0.005)	3.283 (0.049)	2.553 (0.037)
	SR-GP	8.06 (0.730)	5.602 (0.408)	5168 (360)	2611 (95)	0.563 (0.026)	0.411 (0.014)	4.504 (0.086)	3.593 (0.045)
	SR-Gomea	6.309 (0.303)	4.539 (0.128)	4885 (251)	2931 (113)	0.502 (0.010)	0.363 (0.004)	3.362 (0.033)	2.626 (0.057)
	SR-Gomea-op	6.297 (0.267)	4.543 (0.155)	4815 (250)	2852 (160)	0.514 (0.01)	0.371 (0.004)	3.238 (0.059)	2.488 (0.044)
	B-GLM	6.468 (0.531)	4.58 (0.241)	6080 (266)	4228 (124)	0.498 (0.011)	0.360 (0.004)	4.464 (0.068)	3.648 (0.026)
	B-GLM-int	6.356 (0.513)	4.474 (0.235)	5151 (293)	2950 (140)	0.497 (0.01)	0.360 (0.004)	4.282 (0.065)	3.537 (0.026)
	XGBoost	5.571 (0.377)	4.084 (0.175)	4667 (346)	2673 (193)	0.440 (0.011) 0.312 (0.003) 3.257 (0.056)	2.545 (0.028)
	SVM	6.310 (0.645)	4.382 (0.268)	4953 (272)	2968 (147)	0.503 (0.011)	0.356 (0.004)	4.493 (0.089)	3.612 (0.037)
	MLP	6.140 (0.462)	4.392 (0.212)	4867 (347)	2916 (205)	0.463 (0.013)	0.330 (0.008)	3.170 (0.052) 2.445 (0.029)
	BH-GLM	6.102 (0.441)	4.380 (0.187)	4926 (300)	2930 (130)	0.476 (0.019)	0.354 (0.01)	3.873 (0.133)	3.057 (0.091)
	BH-GLM-int	6.004 (0.46)	4.315 (0.191)	4925 (301)	2929 (131)	0.474 (0.02)	0.352 (0.011)	3.871 (0.13)	3.056 (0.088)
	SR-trad	6.258 (0.509)	4.488 (0.208)	5219 (330)	3176 (209)	0.510 (0.018)	0.373 (0.011)	3.961 (0.031)	3.162 (0.026)
	SR-max	6.186 (0.469)	4.44 (0.224)	4889 (293)	3095 (366)	0.507 (0.011)	0.369 (0.008)	3.528 (0.042)	2.791 (0.041)
	HSR-naive-trad 6.472 (0.462)	4.644 (0.261)	6429 (307)	3304 (237)	0.584 (0.067)	0.410 (0.031)	4.156 (0.184)	3.064 (0.162)
	HSR-naive-max 6.386 (0.473)	4.582 (0.252)	6328 (307)	3170 (243)	0.545 (0.018)	0.401 (0.026)	4.062 (0.078)	3.002 (0.09)
	HSR-trad	5.921 (0.54)	4.250 (0.226)	4840 (308)	2930 (153)	0.475 (0.012)	0.343 (0.011)	3.30 (0.084)	2.547 (0.084)
	HSR-max	5.80 (0.507)	4.210 (0.282)	4844 (304)	2933 (149)	0.470 (0.011)	0.340 (0.011)	3.277 (0.082)	2.532 (0.087)
		Satellite (24)	Wind (9)	Breast tumor (6)	Music (78)
	OLS	1.213 (0.009)	1.02 (0.011)	3.289 (0.104)	2.521 (0.077)	10.023 (0.036) 7.88 (0.027)	0.465 (0.039) 0.35 (0.03)
	Decision tree	1.061 (0.022)	0.621 (0.017)	3.839 (0.112)	2.980 (0.103)	9.844 (0.039)	7.632 (0.03)	0.705 (0.066)	0.497 (0.053)
	GAM-splines	0.90 (0.009)	0.624 (0.007)	3.082 (0.084)	2.366 (0.057)	9.663 (0.047)	7.537 (0.034)	0.898 (0.101)	0.678 (0.069)
	EBM	0.851 (0.035)	0.575 (0.021)	3.140 (0.089)	2.398 (0.061)	9.519 (0.049)	7.401 (0.04)	0.60 (0.036)	0.441 (0.035)
	SR-GP	1.628 (0.438)	1.055 (0.095)	3.858 (0.216)	2.979 (0.191)	10.441 (0.277) 8.151 (0.23)	0.71 (0.154)	0.51 (0.102)
	SR-Gomea	1.164 (0.028)	0.897 (0.036)	3.306 (0.128)	2.545 (0.099)	9.988 (0.056)	7.795 (0.044)	0.523 (0.081)	0.379 (0.056)
	SR-Gomea-op	1.102 (0.032)	0.826 (0.029)	3.296 (0.101)	2.534 (0.077)	9.973 (0.049)	7.788 (0.025)	0.499 (0.036)	0.369 (0.036)
	B-GLM	1.213 (0.008)	1.019 (0.012)	3.308 (0.098)	2.527 (0.074)	10.018 (0.033) 7.776 (0.03)	0.473 (0.033)	0.353 (0.03)
	B-GLM-int	1.117 (0.044)	0.905 (0.069)	3.295 (0.098)	2.517 (0.073)	9.995 (0.034)	7.791 (0.026)	0.469 (0.045)	0.358 (0.031)
	XGBoost	0.667 (0.033) 0.35 (0.016)	3.084 (0.075)	2.365 (0.050) 9.435 (0.048) 7.288 (0.041) 0.507 (0.046)	0.367 (0.037)
	SVM	1.261 (0.028)	1.003 (0.022)	3.307 (0.102)	2.521 (0.073)	10.045 (0.036) 7.86 (0.029)	0.472 (0.042)	0.348 (0.03)
	MLP	0.789 (0.047)	0.448 (0.03)	3.076 (0.087) 2.367 (0.062)	9.67 (0.04)	7.519 (0.023)	0.498 (0.042)	0.36 (0.04)
	BH-GLM	0.986 (0.036)	0.598 (0.034)	3.297 (0.101)	2.538 (0.079)	9.76 (0.035)	7.611 (0.031)	0.472 (0.035)	0.356 (0.027)
	BH-GLM-int	0.952 (0.029)	0.553 (0.028)	3.291 (0.106)	2.531 (0.084)	9.751 (0.036)	7.598 (0.033)	0.467 (0.037)	0.353 (0.027)
	SR-trad	1.175 (0.05)	0.948 (0.059)	3.342 (0.104)	2.568 (0.080)	10.096 (0.064) 7.917 (0.051)	0.543 (0.081)	0.401 (0.049)
	SR-max	1.018 (0.04)	0.784 (0.035)	3.176 (0.081)	2.445 (0.059)	10.03 (0.079)	7.864 (0.06)	0.476 (0.045)	0.359 (0.037)
	HSR-naive-trad 1.012 (0.064)	0.645 (0.068)	3.547 (0.041)	2.731 (0.050)	11.932 (0.16)	9.289 (0.241)	0.524 (0.045)	0.372 (0.02)
	HSR-naive-max 0.991 (0.034)	0.631 (0.031)	3.562 (0.068)	2.711 (0.046)	11.915 (0.108) 9.267 (0.185)	0.507 (0.056)	0.372 (0.039)
	HSR-trad	0.95 (0.034)	0.554 (0.032)	3.205 (0.074)	2.457 (0.054)	9.727 (0.056)	7.595 (0.06)	0.497 (0.05)	0.366 (0.031)
	HSR-max	0.934 (0.056)	0.529 (0.062)	3.198 (0.074)	2.455 (0.054)	9.662 (0.061)	7.531 (0.049)	0.471 (0.053)	0.351 (0.027)
		House (5)	Wine (8)	Toxicity (6)	Gas (7)
	OLS	41563 (1270)	24354 (202)	0.754 (0.02)	0.586 (0.012)	1.256 (0.097)	0.949 (0.072)	8.112 (0.133)	5.796 (0.081)
	Decision tree	35752 (1199)	20035 (401)	0.753 (0.015)	0.596 (0.013)	1.394 (0.12)	1.059 (0.086)	7.705 (0.127)	5.638 (0.068)
	GAM-splines	33460 (1405)	18634 (212)	0.728 (0.029)	0.565 (0.015)	1.245 (0.106)	0.92 (0.08)	5.993 (0.10)	4.187 (0.038)
	EBM	31062 (1192)	16921 (133)	0.689 (0.019)	0.537 (0.01)	1.20 (0.095)	0.899 (0.06)	5.476 (0.072) 3.763 (0.038)
	SR-GP	56615 (21299) 24687 (2204)	0.857 (0.068)	0.664 (0.055)	1.457 (0.264)	1.087 (0.152)	10.75 (1.107)	8.165 (0.775)
	SR-Gomea	36750 (1693)	20768 (976)	0.742 (0.022)	0.582 (0.015)	1.343 (0.189)	0.978 (0.088)	8.703 (0.144)	6.583 (0.176)
	SR-Gomea-op	36865 (1434)	20916 (788)	0.739 (0.021)	0.579 (0.015)	1.267 (0.116)	0.956 (0.086)	8.742 (0.278)	6.764 (0.227)
	B-GLM	42104 (1406)	22632 (220)	0.753 (0.02)	0.586 (0.012)	1.237 (0.099)	0.934 (0.073)	8.112 (0.134)	5.797 (0.081)
	B-GLM-int	39811 (1375)	20872 (335)	0.751 (0.019)	0.583 (0.012)	1.246 (0.089)	0.941 (0.056)	8.112 (0.137)	5.797 (0.081)
	XGBoost	29630 (1237) 15733 (120)	0.68 (0.014)	0.531 (0.006) 1.157 (0.129) 0.872 (0.092) 5.705 (0.190)	4.001 (0.126)
	SVM	44879 (1676)	21201 (228)	0.748 (0.018)	0.582 (0.01)	1.286 (0.195)	0.959 (0.132)	6.954 (0.372)	4.875 (0.453)
	MLP	36004 (851)	20174 (396)	0.757 (0.071)	0.587 (0.051)	1.280 (0.153)	0.973 (0.112)	6.023 (0.239)	4.223 (0.134)
	BH-GLM	35212 (1661)	19444 (708)	0.741 (0.017)	0.576 (0.013)	1.235 (0.095)	0.933 (0.069)	6.87 (0.282)	4.768 (0.159)
	BH-GLM-int	35039 (1820)	19079 (443)	0.737 (0.017)	0.574 (0.014)	1.237 (0.098)	0.934 (0.073)	6.87 (0.282)	4.769 (0.159)
	SR-trad	37412 (2150)	21067 (524)	0.744 (0.024)	0.584 (0.015)	1.253 (0.089)	0.953 (0.078)	8.153 (0.686)	6.131 (0.515)
	SR-max	34716 (2049)	18976 (491)	0.731 (0.019)	0.571 (0.016)	1.233 (0.106)	0.935 (0.071)	7.395 (0.553)	5.467 (0.565)
	HSR-naive-trad 37518 (1563)	19592 (433)	0.735 (0.021)	0.575 (0.011)	1.363 (0.112)	0.989 (0.08)	7.46 (0.708)	5.557 (0.607)
	HSR-naive-max 37640 (1741)	19294 (472)	0.725 (0.019)	0.566 (0.009)	1.266 (0.071)	0.944 (0.057)	6.621 (0.15)	4.762 (0.152)
	HSR-trad	33542 (1127)	17964 (258)	0.724 (0.019)	0.566 (0.013)	1.243 (0.097)	0.933 (0.047)	7.181 (0.618)	5.32 (0.776)
	HSR-max	33102 (833)	17403 (155)	0.712 (0.014)	0.559 (0.011)	1.214 (0.065)	0.921 (0.038)	6.421 (0.150)	4.662 (0.153)

a

Table 4 .

 4 8: Results obtained by cluster-specific interpretable models

Table 4 . 9

 49

: Results on the classification datasets

  42 7685 14658 26606.63 3513.14 20949 37113 15920.91 1721.15 13382 20945 5627.28 1225.75 659

	8438	1 35.06 16.28 4.67 1 34.13	1.03 3 2.97 0.07 2.6 3	84.03 521.02 188.75 100.25 79.76 549.63
		5.35	0.34	96.58
		5.2 31.4 16.84	2.68 3 2.89	71.97 427.78 286.83
		4.54	0.06	70.51
		6.73 2.08 30.3 17.13	0.27 1.15 3 2.97	104.72 79.76 615.78 210.17	0.46 0 1 0.17
		%age heavy vehicles 18.22	Right shoulder width 2.91	Altitude 274.53	Presence of ramps 0.30

Cette thèse est accessible à l'adresse : https://theses.insa-lyon.fr/publication/2022ISAL0096/these.pdf © [T. Véran], [2022], INSA Lyon, tous droits réservés

Laboratoire d'InfoRmatique en Image et Systèmes d'information: https://liris.cnrs.fr/

https://www.data-newroad.com/

A model is said to be Pareto optimal if there is no alternative model that can improve one of its objective function without deteriorating the others.

Note that, in[START_REF] Scott M Lundberg | From local explanations to global understanding with explainable ai for trees[END_REF], Lundberg et al. call this process a supervised clustering. However, this usually refers to clustering with access to a "teacher" that knows the right cluster for some of the instances[START_REF] Awasthi | Supervised clustering[END_REF]. In our case, our approach for clustering is unsupervised.

https://xgboost.readthedocs.io/en/stable/

The test dataset materializes a production context as our model did not use these data for training.

https://github.com/kvoyager/GmdhPy

https://docs.sympy.org/latest/modules/simplify/simplify.html

https://scikit-learn.org/stable/

https://imbalanced-learn.org/stable/references/generated/ imblearn.combine.SMOTEENN.html

https://docs.pymc.io/en/v3/

https://www.statsmodels.org/stable/mixed_glm.html

https://pygam.readthedocs.io/en/latest/

https://gplearn.readthedocs.io/en/stable/

https://github.com/marcovirgolin/GP-GOMEA

https://xgboost.readthedocs.io/en/latest/python/python_intro.html

Acknowledgments

The doctoral student was part of the research groups for distributed systems, information retrieval and mobility (DRIM) and data base (BDThis work has received financial support from the French Association for Research and Technology (ANRT).

x 3

x 6 overall PD(x 3 ) < -0.0015 overall PD(x introduces highly correlated combinations of terms. From model 3 to model 4, a sharp increase in complexity for a small gain in performance should alert the user to the risk of no longer understanding the inner workings of model 4: time must be spent at studying the various partial effects before deciding if the model can still be trusted.

Specificities and benefits of the ranking-by-complexity approach

Thanks to our complexity metric (see Eq. 3.3), model 3 does not dominate model 2 even though they have the same number of terms and their terms have equal complexities. If we had not penalize collinearities, then model 2, which is of high interest to field experts, would not have been included in the Pareto front. In this sense, the ranking-by-complexity favors a progressive analysis of numerous instructive models. Among Fig. 4.6 models, some can attain similar predictive performances while bringing out different effects of the explanatory variables. This can be understood from the point of view of the Rashomon effect [START_REF] Breiman | Statistical modeling: The two cultures (with comments and a rejoinder by the author)[END_REF] which characterizes problems