Perrine Lacroix
email: perrine.lacroix@universite-paris-saclay.fr

Mélina Gallopin
email: melina.gallopin@universite-paris-saclay.fr

Marie-Laure Martin
email: marie_laure.martin@agroparistech.fr

AN OVERVIEW OF VARIABLE SELECTION PROCEDURES USING REGULARIZATION PATHS IN HIGH-DIMENSIONAL GAUSSIAN LINEAR REGRESSION

Keywords: Variable selection, Gaussian linear regression, High-dimension, Regularization path, Comparison study

Current high-throughput technologies provide a large amount of variables to describe a phenomenon. Only a few variables are generally sufficient to answer the question. Identify them in a highdimensional Gaussian linear regression model is the one of the most-used statistical methods. In this article, we describe step-by-step the variable selection procedures built upon regularization paths. Regularization paths are obtained by combining a regularization function and an algorithm. Then, they are combined either with a model selection procedure using penalty functions or with a sampling strategy to obtain the final selected variables. We perform a comparison study by considering three simulation settings with various dependency structures on variables. In all the settings, we evaluate (i) the ability to discriminate between the active variables and the non-active variables along the regularization path (pROC-AUC), (ii) the prediction performance of the selected variable subset (MSE) and (iii) the relevance of the selected variables (recall, specificity, FDR). From the results, we provide recommendations on strategies to be favored depending on the characteristics of the problem at hand. We obtain that the regularization function Elastic-net provides most of the time better results than the ℓ 1 one and the lars algorithm has to be privileged as the GD one. ESCV provides the best prediction performances. Bolasso and the knockoffs method are judicious choices to limit the selection of non-active variables while ensuring selection of enough active variables. Conversely, the data-driven penalties considered in this review are not to be favored. As for Tigress and LinSelect, they are conservative methods.

Introduction

Recent scientific advances allow us to have access to large-scale data: the size of the data sets is exploding, as well as the complexity of each of them. For instance, in genomics, to describe the molecular activities, microarrays and RNA-sequencing technologies provide a quantification of the expression of all the genes simultaneously, making it easier to study their interactions. Studies on genetic associations are diversifying, considering a wide range of phenotypes, including gene expression or proteomic and metabolomic data. Other fields are concerned with such massive data sets, such as medicine with the development of imaging, treatment or disease monitoring, or market strategy. In a statistical point of view, the number of parameters to estimate explodes and reduction of dimension is required to select only relevant variables and summarize the redundant information for a given model. In this review, we focus on the variable selection procedures in high-dimensional linear Gaussian regression models. The considered dataset with a number of variables p close to or slightly higher than the number of observations n is a real challenge since it hampers the use of the traditional estimation methods. A regularization of the cost function is required so that only a subset of variables is selected to explain the response variable.

In most reviews on variable selection in high-dimensional Gaussian linear regression, a focus is done on the construction of the regularization path. It is based on the minimization of a cost function penalized by a regularization function and provides an order on variables. [START_REF] Wainwright | Sharp thresholds for noisy and high-dimensional recovery of sparsity using l 1 -constrained quadratic programming[END_REF] provides a meticulous theoretical analysis of the ℓ 1 regularization function. In particular, for a given number of active variables, the author discusses the choice of the number of observations to ensure asymptotic properties to recover these active variables. [START_REF] Bühlmann | High-dimensional variable screening and bias in subsequent inference, with an empirical comparison[END_REF] compared several regularization functions in a simulation study by using semi-real datasets. In their simulation design, they considered several numbers of observations, of variables and of active variables. They also modified the signal-to-noise ratio and considered two scenarios of variable correlations. The results of the different regularization functions are inspected with ROC curves and partial ROC curves when 0.5 × n and 0.9 × n variables are selected. [START_REF] Wang | High-dimensional regression in practice: an empirical study of finite-sample prediction, variable selection and ranking[END_REF] compared a large set of regularization functions with a simulation design similar to [START_REF] Bühlmann | High-dimensional variable screening and bias in subsequent inference, with an empirical comparison[END_REF]. They evaluated both prediction and variable identification but the main difference with our investigations is that the only considered model selection procedure applied from the regularization path is the cross-validation one. Finally some reviews considered different contexts. [START_REF] Wu | A selective review of robust variable selection with applications in bioinformatics[END_REF] were interested in robust variable selection strategies when heavy-tailed errors and outliers in response variables exist. They discussed the different steps from the modification of the least squares function to the choice of the parameters for the model selection through a presentation of algorithms accounting for outliers. [START_REF] Vinga | Structured sparsity regularization for analyzing high-dimensional omics data[END_REF] considered a variety of models from survival models to generalized linear models, frequently used in biomedical research. [START_REF] David | A review on variable selection in regression analysis[END_REF] considered a wide range of model structures (linear, grouped, additive, partially linear and non-parametric) and discussed three main categories of algorithms for the variable selection.

Our review distinguishes itself from the previous ones since we propose an evaluation of both construction of regularization paths and choice of the final selected variables. This leads to 33 combinations. Moreover, for model selection procedure, we add in this review non-asymptotic methods which are generally not considered. To construct the regularization path, we test two regularization functions (Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and Elastic-Net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]) combined with two algorithms (LARS [START_REF] Efron | Least angle regression[END_REF] and gradient descent algorithm [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]). Each regularization path provides a collection of variable subsets. To choose one of them, we compare model selection and variable identification approaches. On the one hand, the model selection uses penalization criteria of the least squares (eBIC [START_REF] Chen | Extended bayesian information criteria for model selection with large model spaces[END_REF], data-driven calibration strategies [START_REF] Birgé | Minimal penalties for gaussian model selection[END_REF][START_REF] Lebarbier | Detecting multiple change-points in the mean of gaussian process by model selection[END_REF][START_REF] Baudry | Slope heuristics: overview and implementation[END_REF][START_REF] Arlot | Minimal penalties and the slope heuristics: a survey[END_REF] and LinSelect [START_REF] Baraud | Gaussian model selection with an unknown variance[END_REF][START_REF] Giraud | High-dimensional regression with unknown variance[END_REF]). On the other hand, the variable identification methods (ESCV [START_REF] Lim | Estimation stability with cross-validation (escv)[END_REF], Bolasso [START_REF] Bach | Bolasso: model consistent lasso estimation through the bootstrap[END_REF], Stability Selection [START_REF] Meinshausen | Stability selection[END_REF], Tigress [START_REF] Haury | Tigress: trustful inference of gene regulation using stability selection[END_REF] and the knockoffs method [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF]) use sampling strategies to stabilize the selected variable subset while limiting the selection of non-active variables. Methods based on multiple testing procedures [START_REF] Kos | On the asymptotic properties of slope[END_REF] and Bayesian approaches are not included in this review. We refer the readers to [START_REF] Celeux | Regularization in regression: comparing bayesian and frequentist methods in a poorly informative situation[END_REF] for an empirical comparison of frequentist and Bayesian points of view. Lastly, we assume no prior knowledge between interactions, spatial localization and chronological information and refer to [START_REF] Bondell | Consistent high-dimensional bayesian variable selection via penalized credible regions[END_REF][START_REF] Razaghi-Moghadam | Supervised learning of gene-regulatory networks based on graph distance profiles of transcriptomics data[END_REF] for such approaches.

After a description of the methods and some of their theoretical properties, we compare them in a simulation study by considering three settings. In the first one, the variables are independent and are drawn from a Gaussian distribution. It allows a comparison in the theoretical framework used to develop them. In the second setting, two structures of the correlation between variables are considered to evaluate how dependencies usually observed affect the methods. Observations are generated according to a Gaussian linear model, the most favorable case where assumptions broadly hold. Finally, the third setting mimics the biological complexity of transcription factor regulations. Observations are generated using the FRANK algorithm [START_REF] Carré | Reverse engineering highlights potential principles of large gene regulatory network design and learning[END_REF].

In these three settings, performances of the methods are evaluated for their prediction performance and for their ability to identify the active variables. To discriminate active variables to the others, we use the pROC-AUC metric. We use the mean squared errors (MSE) to measure the prediction performance, and the recall, specificity and false discovery rate (FDR) metrics to assess the quality of the selected variables in terms of active variables. As [START_REF] Bühlmann | High-dimensional variable screening and bias in subsequent inference, with an empirical comparison[END_REF][START_REF] Wang | High-dimensional regression in practice: an empirical study of finite-sample prediction, variable selection and ranking[END_REF], we notice that there is no unambiguous winner among all the studied approaches. Our goal is to provide recommendations for a judicious choice of a method according to the application. In particular, Elastic-Net should be preferred to Lasso for the regularization function, as well as the lars algorithm. Moreover, to ensure a prediction ability, ESCV and the knockoffs seem to be the most judicious choices. If the goal is to recover the active variables, ESCV and eBIC are preferable whereas Bolasso, the knockoffs and LinSelect should be privileged to limit the non-active variables in the selected subset.

The rest of the paper is organized as follows. Section 2 describes the statistical framework and all the methods we compare in this review. Section 3 presents the simulation settings, technical aspects about the implementation of methods and the evaluation metrics we consider. Section 4 is devoted to all the results. Section 5 offers practical recommendations to choice the best methods. Lastly, a discussion is provided in Section 6.

Methods

Statistical framework

For the sequel, the norms |.| 0 , |.| 1 and ||.|| are respectively the usual norms 0, 1 and 2 on R q , for any q ∈ N * . It means that for a vector β ∈ R q ,

|β| 0 = q j=1 1 {βj ̸ =0} ; |β| 1 = q j=1 |β j |; ||β|| = q j=1 β 2 j .
We consider a Gaussian linear regression model where the response variable Y is explained by a linear combination of p variables X = (X 1 , . . . , X p ):

Y = Xβ * + ε.
The parameter β * ∈ R p is the vector of unknown parameters and ε follows a centered Gaussian distribution with an unknown variance denoted σ 2 .

To estimate the parameters β * and σ 2 , we observe y i ∈ R and the variables (x i1 , ..., x ip ) ∈ R p for i ∈ {1, ..., n}.

We suppose that observations are independent. We consider the high-dimensional framework with p ∼ n or p ≫ n, preventing the traditional least squares estimation. In this context, we assume that only a small number of variables among the p ones explains the response variable. These variables are associated to a non-zero coefficient in β * and are named active variables. Under this sparsity assumption, the target is for t > 0:

min β∈R p :|β|0≤t ||Y -Xβ|| 2 ,
where its associated Lagrangian form is for λ > 0:

min β∈R p ||Y -Xβ|| 2 + λ|β| 0 . (1)
The proof of the equivalence and the link between t and λ are provided in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

Determining the hyperparameter λ is one of the major issues and the challenge lies in its calibration to adjust a trade-off between sparsity and a good linear adjustment. A large value of λ provides a small subset of variables (assumption of sparsity satisfied) but it might correspond to a linear adjustment far from the response variable. A small value of λ provides a large subset of variables (assumption of sparsity not satisfied) but it might correspond to a linear adjustment close to the response variable. Moreover, the criterion (1) being non-convex, the existence and the uniqueness of the solution are not guaranteed. So, as presented in [START_REF] Vinga | Structured sparsity regularization for analyzing high-dimensional omics data[END_REF], Equation ( 1) is replaced with the optimization problem:

min β∈R p ||Y -Xβ|| 2 + λF (β) , (2)
where F is a continuous and convex regularization function satisfying the existence of a minimum for any λ.

Regularization functions

Several regularization functions exist and we present the most used. The first one is the Lasso regularization [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] with

F (β) = |β| 1 .
If λ is well chosen, it provides a consistent estimator of β * . This procedure achieves the best trade-off between regularity (convexity, reasonable computationally solving) and sparsity for independent variables. However, when some variables are correlated, Lasso tends to select randomly only one of them rather than selecting none or all of them. A solution to this problem is the Adaptive Lasso procedure where each variable is weighted with respect to an initial estimator [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF].

Another well-known regularization function is the Ridge regularization [START_REF] Hoerl | Ridge regression[END_REF]:

F (β) = ||β|| 2 .
In addition to taking variable dependencies into account, Ridge provides a strictly convex and derivable optimization problem with an explicit estimator of β * . However, this estimator is not sparse. To combine sparsity and correlated variables, the Elastic-Net regularization [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] has been proposed:

F (β) = (1 -α)|β| 1 + α||β|| 2 ,
where the parameter α controls the trade-off.

When prior knowledge on variable dependencies is available, there exist other regularization functions, not considered here: the Group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], Overlap Group Lasso [START_REF] Jacob | Group lasso with overlap and graph lasso[END_REF], Hierarchical Group Lasso [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF], double sparse Lasso [START_REF] Li | Sharp minimax optimality of lasso and slope under double sparsity assumption[END_REF] and fused Lasso [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF].

Regularization path construction for Lasso and Elastic-Net

The optimization problem (2) has generally no explicit solution and requires a computational approach. A first algorithm is LARS [START_REF] Efron | Least angle regression[END_REF] where one variable is added at each step and each step corresponds to a value of λ. A grid Λ of λ is then naturally provided. Briefly speaking, the first subset contains the variable X j which has the largest absolute correlation with Y . The second subset contains exactly two variables: X j and the variable which is the most correlated with the residuals of the regression of Y on X j . LARS provides an exact solution of the optimization problem with nested subsets. It is an important property for theoretical considerations. A second algorithm is based on the gradient descent method [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]. This algorithm constructs a regular grid Λ of a given size by starting with the largest λ corresponding to the first nonempty variable subset. Then, a variable subset is obtained for each λ of this grid by solving (2) with the cyclic coordinate descent method. In contrast to LARS, the gradient descent method provides a proxy of the optimization problem with independent solutions along the grid. Consequently, the gradient descent method may provide a richer collection with several models of the same dimension.

Whatever the choice of both regularization function and algorithm, a collection of variable subsets (m λ ) λ∈Λ is obtained at the end of the regularization path construction. Each m λ is associated with an estimator of β * , however this estimator being known to be biased [START_REF] Meinshausen | Relaxed lasso[END_REF], it is commonly replaced with the least-squares estimator calculated on the variable subset m λ and is denoted βλ [START_REF] Connault | Calibration d'algorithmes de type Lasso et analyse statistique de données métallurgiques en[END_REF]. The number of non-zero coefficients of βλ or in an equivalent manner, the number of variables in m λ is denoted D λ . The next step consists in selecting a variable subset among the collection either with a model selection criterion or with a variable identification procedure.

Model selection

To select the best subset m λ, model selection approaches consist of minimizing a penalized loss function in λ ∈ Λ:

γ(m λ ) + pen(n, p, D λ ). (3)
The loss function γ(m λ ), quantifying the quality of the model fit, is either the least-squares function ||Y -X βλ || 2 or the deviance -2 log(L(Y, X; βλ , σ2 λ )), where L is the likelihood function calculated with βλ and σ2 λ , the empirical estimators associated to the m λ . The penalty function pen(n, p, D λ ) accounts for the model complexity and the characteristics of the sample: higher the penalty values, smaller the number of selected variables and farther the linear combination X βλ to the response variable Y .

Asymptotic criteria. The first criteria are asymptotic: their properties are verified when the sample size n tends to infinity. In this review, we focus on the more recent asymptotic criterion, called eBIC [START_REF] Chen | Extended bayesian information criteria for model selection with large model spaces[END_REF], used to get a consistent estimator by penalizing the deviance by:

pen eBIC (n, p, D λ ) = D λ log(n) + 2δ log( p D λ ), (4)
where δ is a value in [0, 1].

Non-asymptotic criteria. In a practical consideration, having guarantees for n going to infinity has no sense [START_REF] Giraud | High-dimensional regression with unknown variance[END_REF] and applying criteria with properties confirmed for any fixed sample n size is more relevant. Introduced by Birgé and Massart [START_REF] Birgé | Gaussian model selection[END_REF], the goal of non-asymptotic criteria is to achieve the risk oracle:

inf λ∈Λ E[ ||Xβ * -X βλ || 2 ],
and instead of getting asymptotic equality of the kind

P   lim n→+∞ E[ ||Xβ * -X βλ || 2 ] inf λ∈Λ E[ ||Xβ * -X βλ || 2 ] = 1   = 1,
they get an inequality holding for any value of n:

E[ ||Xβ * -X βλ || 2 ] ≤ C n inf λ∈Λ {E[ ||Xβ * -X βλ || 2 ] } + R n , (5)
where C n ≈ 1 at least for n large and R n is small comparable to the risk oracle. The selected model is the minimizer of Equation ( 3) where the loss function is the least-squares function and two penalty functions, which do not require the knowledge of the variance, are available.

The first penalty is a data-driven penalty [START_REF] Birgé | Minimal penalties for gaussian model selection[END_REF]:

pen Data-driven (n, p, D λ ) = 2κD λ 2.5 + log p D λ , (6)
where the constant 2.5 has been proposed in a context of changepoint detection in a signal [START_REF] Lebarbier | Detecting multiple change-points in the mean of gaussian process by model selection[END_REF]. The constant κ is calibrated from the sample. For that, two strategies are proposed. The first one is the slope heuristics: assuming that the least-squares function is linear in D λ 2.5 + log( p D λ ) as soon as D λ is large enough (see Figure 2 of [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF]), the constant κ is then equal to the estimated slope. The second strategy is the dimension jump: assuming the existence of κ * such that for all the values smaller than κ * , the associated model has a very high dimension, whereas for all the values greater than κ * , the associated model has a reasonable dimension (see Figure 1 of [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF]), the constant κ is then equal to the estimated κ * . For more practical and theoretical details, we refer the reader to [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF] and to the survey [START_REF] Arlot | Minimal penalties and the slope heuristics: a survey[END_REF].

The second penalty function is LinSelect proposed in [START_REF] Baraud | Gaussian model selection with an unknown variance[END_REF] and generalized for a high dimensional context in [START_REF] Giraud | High-dimensional regression with unknown variance[END_REF]. It is built from the empirical estimator of the variance onto each m λ :

pen LinSelect (n, p, D λ ) = 1.1 × n -D λ n -D λ -1 Ψ D λ + 1, n -D λ -1, e -L λ , (7)
where the L λ are weights satisfying some properties and the function Ψ[D, N, q] is the unique solution of the equation:

ϕ [D, N, Ψ(D, N, q)] = q,
where ϕ[D, N, x] is defined for x ≥ 0:

ϕ[D, N, x] = 1 D E max 0, χ 2 D -x χ 2 N N ,
for χ 2 D and χ 2 N two independent χ 2 random variables with degrees of freedom D and N respectively.

Variable identification

The high-dimensional framework usually leads to unstable results: addition, suppression, or modification of some observations could radically change the selected variable subset. For prediction, different sets of variables can give the same prediction performances. However, when the objective is the identification of the active variables, this instability is a drawback.

To circumvent this sampling uncertainty, the idea is to work with perturbed datasets generated from the original sample. Cross-validation [START_REF] Allen | The relationship between variable selection and data augmentation and slow feature analysis[END_REF][START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF] is commonly proposed. It consists in splitting K times the original sample into a training set and a testing set. The training set is used to calculate an estimator βk λ and the testing set is used to evaluate the mean squared error. The selected model is the one which minimizes the mean squared error in λ. Applying cross-validation in a high-dimensional context is computationally expensive and known to be unstable. An alternative is ESCV [START_REF] Lim | Estimation stability with cross-validation (escv)[END_REF] which estimates the instability along the regularization path with the K perturbed datasets and selects the variable subset which minimizes the instability in λ:

1 K K k=1 ||X βk λ -1 K K ℓ=1 X βℓ λ || 2 || 1 K K ℓ=1 X βℓ λ || 2 .
Sampling strategy is also a solution. Two widely used approaches are Bolasso [START_REF] Bach | Bolasso: model consistent lasso estimation through the bootstrap[END_REF] and Stability Selection [START_REF] Meinshausen | Stability selection[END_REF]. They mainly differ in their sampling strategy: Bolasso generates datasets of n data uniformly chosen with replacement among the original sample, whereas Stability Selection generates datasets of n 2 distinct data randomly chosen and includes also the complement of each generated dataset in the sampling strategy to limit the subsampling effects. In Stability Selection, they also propose a random perturbation in the Lasso regularization:

F (β) = p j=1 |β j | w j ,
where w j ∼ U([θ, 1]) with θ > 0. Sampling strategies get an occurrence frequency of each variable and those having the highest occurrence frequencies are retained to constitute the final variable subset. Tigress method [START_REF] Haury | Tigress: trustful inference of gene regulation using stability selection[END_REF] modifies the calculation of the occurrence frequency by averaging on the grid.

The last type of variable identification method is the knockoffs method [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF], which controls the False Discovery Rate (FDR). This method starts with the construction of a matrix X such that X and X have the same covariance structure with Xj the least correlated to X j . It is done through linear algebra tools [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF][START_REF] Jean Candès | Panning for gold: Model-free knockoffs for high-dimensional controlled variable selection[END_REF]. Then, a regularization path is constructed on the augmented matrix X X of size n × 2p where the active variables are expected to be selected very earlier than their copy. Let denote

W j = max Z j , Zj × sign Z j -Zj , (8)
where Z j and Zj correspond to the largest λ for which X j and Xj are selected respectively. A positive value of W j states that X j is selected before its copy Xj and a large positive value indicates that X j is selected rapidly. Let q be the target FDR, the final variable subset is composed by the X j such that W j ≥ T with:

T = min t ∈ {|W j |, j = 1, ..., p} \ {0}, 1 + #{j : W j ≤ -t} min (1, #{j : W j ≥ t}) ≤ q .
3 Comparison study

Three simulation settings

The design of the simulation study is composed of three simulation settings to study the behavior of the methods with respect to the dependency structures between the variables.

Simulation under independent design. This is the simplest setting where the high-dimensional framework is the single handicap [START_REF] Fan | Adaptive robust variable selection[END_REF][START_REF] Wang | High-dimensional regression in practice: an empirical study of finite-sample prediction, variable selection and ranking[END_REF]. The matrix X is simulated by concatenation of p independent standard Gaussian vectors of size n.

The number of non-zero coefficients of the vector β * is generated from a uniform variable on integers between 10 and 15. Theirs values are generated from a uniform distribution between 0.5 and 2 and the response variable Y is defined by Y = Xβ * + ε, where ε ∼ N (0, I n ). The independent design is a benchmark in this comparison study since the statistical properties of each method must be verified on the independent structure.

Simulation under a Gaussian graphical model. An equivalence exists between the network inference by Gaussian Graphical model and support estimation in Gaussian linear regressions [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF]. An edge between the nodes i and j in the network states that X i and X j are dependent and either X j is on the support of the regression when X i is the response variable; or X i is on the support of the regression when X j is the response variable. In this direction, datasets with dependency structure are simulated using a Gaussian Graphical model. More precisely, a sample of size n is generated from a (p + 1) multivariate centered Gaussian distribution with covariance matrix Σ, where the dependency structure is encoded in the precision matrix Σ -1 [START_REF] Lauritzen | Graphical Models[END_REF][START_REF] Córdoba | On generating random gaussian graphical models[END_REF]. The response variable Y is chosen as a column of the (p + 1) multivariate centered Gaussian and the remaining columns constitute the matrix X of size n × p. It differs from the previous papers [START_REF] Zou | Regularized simultaneous model selection in multiple quantiles regression[END_REF][START_REF] Wang | Quantile regression for analyzing heterogeneity in ultra-high dimension[END_REF][START_REF] Peng | An iterative coordinate descent algorithm for high-dimensional nonconvex penalized quantile regression[END_REF][START_REF] Wang | High-dimensional regression in practice: an empirical study of finite-sample prediction, variable selection and ranking[END_REF], where the response variable is simulated once the matrix X is fixed, but our simulation choice has been motivated by applications such that regulatory network inference where transcription factors can be response variables and also active variables.

We consider two graph patterns:

• Cluster: the precision matrix is simulated as a block diagonal matrix with B blocks of equal size. The response variable Y is defined as the first variable. • Scale-free: a few variables have a lot of neighbors in the network while all the others have a few ones. We consider two response variables corresponding to the variables having the highest and the smallest number of neighbors. These simulation designs are called scale-free-max and scale-free-min respectively.

Simulation under a dynamical process. It is the most realistic setting, based on the algorithm FRANK [START_REF] Carré | Reverse engineering highlights potential principles of large gene regulatory network design and learning[END_REF] which simulates large networks with characteristics of gene regulatory networks. In this algorithm, variables are categorized into a set of transcription factors that activate or inhibit a set of target genes and the FRANK data are generated from a dynamic process and deviate from the statistical model assumptions, especially the Gaussian distribution. We use FRANK with only transcription factor variables in order to compare the results with those from the other settings. We consider as response variables the variables having the highest and the smallest number of neighbors. These simulation designs are called FRANK-max and FRANK-min respectively.

Simulation parameters. For all the settings, we set n = 150 and p = 199. We generate 100 samples of size 2n to create a training set of size n for the estimation and a validation set of size n to evaluate the methods. For each sample, before being used, the variables are centered and scaled. To generate data from Gaussian graphical model, we use the function huge.generator from the R package huge (version 1.3.4.1). For the cluster design, the block number B equals 5 and the probability of connection within a component is set to the default value 0.3. For the FRANK algorithm, we use the online version available on the website https://m2sb.org/?page=FRANK with p + 1 transcription factors and 2n observations. The number of eigenvalues of the matrix on the unit circle is fixed to 2 and the minimum and maximum of sparsity are set to 1 and 50. Other parameters are set to default values. In subsection 4.8, we study the impact of the high-dimension by increasing the initial dataset from the independent design and from the Gaussian graphical model design to sizes n = 300, 600, 1200.

Investigated methods and their parameters

A total of 16 methods of model selection are defined by the combination of a regularization function (Lasso or Elastic-net) with an algorithm (LARS or the gradient descent method) and a penalty function (eBIC, LinSelect or the 2 data-driven penalties). A total of 17 methods of variable identification are defined: when the sampling strategy is performed before the definition of the grids, 8 methods are defined by the combination of the sampling strategy (Bolasso or Stability Selection) with a regularization function (Lasso or Elastic-net) and an algorithm (LARS or the gradient descent method). When the grid Λ is fixed, the sampling strategy is performed for each λ of the grid, which implies using the gradient descent algorithm, hence 4 methods are defined. Furthermore, we include Tigress, the knockoffs method and ESCV. For the last two, a gradient descent algorithm is used with either Lasso or Elastic-net.

For the LARS algorithm, we use the function enet of the R package elasticnet (version 1.1.1). The maximal number of steps to define the grid size is the default value 50 × min(p, n -1). For the gradient descent method, we use the function glmnet of the R package glmnet (version 3.0) and set the grid size at 1000. Both functions propose the Lasso and elastic-net regularization functions. We set α = 0.5 for Elastic-Net regularization.

To perform model selection, eBIC is implemented with δ = 1. LinSelect is implemented in the function tuneLasso of the R package LINselect (version 1.1.3). The data-driven penalties are calculated by using the function capushe of the R package capushe (version 1.1.1). The parameters are set to the default values except for the minimum percentage of points for the plateau selection set to 0.1.

To perform variable identification, the function escv.glmnet of the R package HDCI (version 1.0.2) is used for the ESCV strategy, with a number of groups K = 10. Bolasso and Stability Selection are implemented with 100 samples. We do not investigate the impact of the presence of a random perturbation in the LASSO regularization function for Stability Selection and set θ = 1. A variable is selected when its occurrence frequency is higher than 0.8. For Tigress, we use the function tigress of the R package tigress (version 0.1.0) with 50 steps for the LARS algorithm. For the knockoffs method, we use the function knockoff.filter with option create.second_order of the R package knockoff (version 0.3.2), we calculate the W j with the function stat.lasso_lambdasmax and set the FDR to 0.1.

Evaluation metrics

First, we evaluate the performance of regularization path constructions. We use the partial area under the receiver operating characteristic curve (pROC-AUC) where the x-axis is the proportion of selected non-active variables among the non-active variables and the y-axis is the proportion of selected active variables among the active variables. A high value of pROC-AUC states that the regularization path is able to discriminate the active variables from the others. Since the lengths of the regularization paths differ according to the chosen regularization and algorithm, the pROC-AUC are calculated by truncating the values of the x-axis at the largest value common to all the regularization paths, to compare them fairly.

Second, we evaluate the prediction performance of the methods by calculating the mean squared errors (MSE) on each validation set ( Ỹ , X):

1 n n i=1 Ỹi -( X βλ ) i 2 , ( 9
)
where βλ is the estimator of β * calculated on the associated training set. As data are centered and scaled, a MSE value lower than 1 means that the method has a prediction ability: the selected variables predict Y better than the empty set.

Finally, we evaluate the variable identification by using three metrics: the recall (the proportion of the selected active variables among the active variables), the specificity (the proportion of the non-active variables not selected among the non-active variables) and the false discovery proportion (the proportion of selected non-active variables among the selected variables). By averaging on the 100 simulated samples, the false discovery proportion becomes the False Discovery Rate (FDR). While the recall and the specificity respectively control the number of active variables which are selected and the number of non-active variables which are non selected, the False Discovery Rate evaluates the quality of the selected variables subset. So, considering each of these metrics when assessing method performance yields different information. Since the objective is to limit the selection of non-active variables while selecting as many active variables as possible, recall and specificity are expected to be close to 1 while the FDR is expected to be low or slightly smaller than the threshold fixed by the knockoffs method.

Results

The first six subsections are dedicated to the results obtained for the independent design and when samples are generated from a Gaussian graphical model when n = 150. In our simulation settings, the sizes of the estimated supports of β * do not take extreme values: the sparsity hypothesis is respected (we refer to Table 1 summarizing the active variable number per simulation setting). The mean of the active variable numbers for the independent design is close to the one for the cluster design (around 12, relatively small). As the independent design is a benchmark in this study, comparison of results from independent and cluster designs allows evaluating the impact of the presence of a dependency structure on variables. The scale-free design allows us to investigate the method behavior with respect to the number of active variables. In scale-free-min, the support size always equals 1. In contrast, for scale-free-max, the support size is 31.41 on average with a high standard deviation (9.70) to fully investigate the methods. Subsection 4.7 is devoted to the FRANK data. As these datasets deviate from the statistical model assumptions, we analyses performances of the methods in an independent way. Lastly, we discuss the behavior of the methods for different values of n to investigate the impact of the high-dimension.

For the sequel, the notation ind is the diminutive of the independent design. The notations GD, E-Net and ℓ 1 denote respectively the gradient descent algorithm, the Elastic-Net regularization function and the Lasso regularization function.

The slope heuristics method and the dimension jump are named slope and jump respectively. Bolasso and Stability Selection are named bol and ss respectively. Lastly, grid and sub denote respectively the strategy when grids are first generated and the strategy when samples are first generated. Throughout this section, we mainly discuss the median of the evaluation metrics obtained from the 100 simulated samples of each scenario.

Size of the selected variable subsets

Tables 2 and3 summarize the number of selected variables by the different methods.

The best combination to get closer to the number of active variables at the step of the regularization path construction is the lars algorithm with the E-Net regularization and when samples are first generated. For the independent setting, we observe that eBIC provides the closest number to the size of the support among the model selection methods. Linselect generally selects the empty set. Concerning the data-driven penalties, the number of selected variables is always significantly higher than the number of active variables. Compared to the slope heuristics strategy, the dimension jump strategy selects less variables and seems to be less sensitive to the choice of the algorithm and the regularization function. Of note, the slope heuristic based on lars and E-Net shows similar results to dimension jump. Among the variable identification methods, ESCV and the knockoffs method provide a number of selected variable very close to the number of active variables, whereas all the other methods select less variables than expected. Tigress is very conservative with less than 5 selected variables in average. Bolasso and Stability Selection with lars and bolasso always select the empty set.

When a dependency structure is introduced, for the model selection methods, we observe that eBIC becomes as conservative as Linselect. The dimension jump still has to be preferred to the slope heuristic. All the variable identification methods show difficulties to manage the dependency between the variables. They select very small subsets of variables in the cluster design and favor the empty subset for the scale-free-min. For the scale-free-max design, ESCV, the knockoff method, bolasso and Stability Selection based on lars and E-Net when samples are first generated (sub) outperform.

4.2 Area under the ROC curve (pROC-AUC)

Figure 1 summarizes the pROC-AUC values of the different regularization paths. The target median value is as high as possible.

For the independent setting, the highest median values are obtained with the lars algorithm: almost 0.99 for the E-Net regularization but the median value decreases around 0.73 for the ℓ 1 regularization. By contrast, the GD algorithm discriminates less well the active and non active variables since the median values are around 0.52 for both regularization functions.

For correlation structure settings, the highest median values are obtained with the lars algorithm combined with the E-Net regularization: 0.84 and 0.98 for cluster and scale-free-max designs respectively and 1 for scale-free-min design.

Lars combined with the ℓ 1 regularization does not stand out from the GD algorithm: 0.55, 0.64 and 0.73 for cluster, scale-free-max and scale-free-min designs, whereas the values of the GD algorithm are around 0.55, 0.55 and 0.72.

In conclusion, whatever the settings, lars algorithm with E-Net regularization is the best combination to discriminate the active variables from the non-active ones.

Mean squared errors (MSE)

Figure 2 and 3 summarize the MSE values of the different methods. The target median value is 0.

For the independent setting, all the model selection strategies are predictive but differences can be observed: the median value is around 0.07 for eBIC and the data-driven penalties and increases to around 0.96 for LinSelect. For the variable identification methods, the smallest median value is around 0.07 for ESCV and the knockoffs method. In contrast, Tigress has a median value ten times higher. Concerning the Bolasso and Stability Selection, the smallest median values are respectively 0.27 and 0.17 and are obtained with the lars algorithm combined with the E-Net regularization when samples are first generated (sub).

For correlation structure settings, among the model selection methods, the data-driven penalties are only predictive for the scale-free-max design with median values between 0.30 and 0.45. For the cluster and scale-free-min designs, the median values of eBIC and LinSelect are elevated, higher than 0.92 but lower than 1, whatever the algorithm and the regularization function. These methods provide median values between 0.6 and 1 for the scale-free-max design. Among the variable identification methods, the median values are also larger than 0.92 but lower than 1 for the cluster and the scale-free-min designs, indicating that having prediction performances is also difficult for these methods. In constrast, for the scale-free-max design, the median values are around 0.33 for ESCV and the knockoffs method; 0.5 and 0.44 for Bolasso and Stability Selection when samples are first generated (sub) and with lars algorithm combined with E-Net regularization. The other methods provide median values between 0.6 and 1.

To conclude, when the variables are independent, the best methods are eBIC, ESCV and the knockoffs method, whatever the algorithm and the regularization function used. When a dependency structure exists, the predictive performances are deteriorated. For the cluster and scale-free-min designs, the methods are not really predictive since the smallest values of MSE are higher than 0.92. For scale-free-max design, the best methods are ESCV, the knockoffs method and the data-driven penalty with the dimension jump strategy when lars algorithm and E-Net penalty are combined. For the independent setting, all the model selection methods except LinSelect, have a median value equals 1, meaning that all the active variables are selected. Of note, eBIC combined with the lars algorithm and the E-Net regularization shows more variability than the other methods. Concerning LinSelect, the median value is always equal to 0.07. Among the variable identification methods, ESCV and the knockoffs method median values equal 1 whatever the algorithm and the regularization used. In contrast, Bolasso and Stability Selection results depend on the algorithm and the regularization used. The highest median values are respectively 0.57 and 0.71 and are obtained with the lars algorithm combined with the E-Net regularization function and when samples are first generated. If E-Net is replaced by the ℓ 1 regularization, the median values go down to 0. Finally Tigress obtains a median value around 0.18.

Recall

When a dependency structure exists, among the model selection methods, the data-driven penalties have a better recall than eBIC and Linselect for the three designs. However, their results depend on the combination of the algorithm and the regularization function used: the slope heuristics strategy combined with the GD algorithm and E-Net is preferable to the dimension jump strategy for the cluster design, but for the scale-free-max the dimension jump combined with lars and E-Net is sligthly better than the slope heuristics strategy combined with lars and E-Net. Of note, the values obtained by the slope heuristics have higher variability than those obtained by dimension jump. Among the variable identification methods, the results depend strongly on the design. For the cluster design, the best median values are obtained with ESCV and Stability Selection combined both with the lars algorithm and the E-Net regularization and when the samples are generated first (sub). Nevertheless, we can observe that the values remain low, around 0.2 and there exists a large variability across the 100 samples. For the scale-free-max design, the knockoffs method and ESCV provide the best results with median values both close to 0.75, whatever the algorithm and the regularization function.

The other methods are clearly worst. For the scale-free-min design, only ESCV has a median value of 1, the others have a median value of 0.

To conclude, when variables are independent, all the model selection methods except LinSelect, ESCV and the knockoffs method recover the active variables. When a dependency structure exists, the results are deteriorated. For the model selection method, the results of the data-driven penalties are better for the three designs but the choice of the algorithm and regularization function depend on the design. An important instability is observed for Bolasso and Stability Selection with respect to the choices of regularization function, the algorithm and the sampling strategy.

Specificity

Figures 6 and7 summary the specificity values of the different methods. The target median value is 1.

For the independent setting, only the data-driven penalties show a specificity different from 1. The results obtained by the different algorithms and regularization functions suggest that the reason is the construction of the regularization path since lars combined with E-Net provides the best specificity. But even with this combination, for both data-driven penalties, some non active variables are selected and it might explained by the size of the estimated support which is larger than the size of the true support. All the variable identification methods have a specificity very close to 1.

We draw the same conclusions when a correlation structure exists. Of note, with slope heuristics strategies, values decrease with more variability.

To conclude, the dependency structure between variables and the choice of both algorithm and regularization function do not impact the specificity. All the methods except the data-driven penalties have a specificity very close to 1. For the data-driven penalties, the estimated support being larger than the true support (see subsection 4.1): it suggests that the penalty term ( 6) is too small.

FDP and FDR

Figures 8 and9 summary the FDP values of the different methods. The target median value is low but not equal 0 or slightly smaller than the threshold fixed by the knockoffs method.

For the independent setting, the median values are generally smaller for the ℓ 1 than E-Net penalty, whatever the choice of the algorithm used. Among the model selection methods, LinSelect has a FDP equals 0 because the method generally selects no variable. The data-driven penalties have a high FDP. As already observed for the specificity, it could be due to a too small penalty term (6) which selects too much variables. Only eBIC has a reasonable estimated FDR at 0.06 and 0.17 with the ℓ 1 and E-Net regularization respectively. For all the variable identification methods, the median value of the FDP is equal to 0 except the knockoffs method where the median value equals 0.08 for an expected control at 0.1.

Unsurprisingly, the presence of a dependency structure impacts negatively the FDP values, especially for the data-driven penalties and ESCV. The ℓ 1 penalty generally provide smaller median values than E-Net penalty. Among the model selection methods, for eBIC and Linselect, the median value is always equal to 0. If we calculate an estimated FDR from only the samples where the selected subset is non empty, we get an estimated FDR higher than 0.8 for both methods and for the cluster and scale-free-min designs. Only for the scale-free-max design, eBIC achieves an estimated FDR around 0.11 whereas Linselect achieves 0.34. The data-driven penalty strategies provide high median values. Among the variable identification methods, for the cluster design, only ESCV and the Stability Selection when the lars algorithm is combined with E-Net and when the samples are generated first provide a median value of FDP different from 0. However there exists a large variability and it leads to an estimated FDR higher than 0.24. For scale-free-max, only the knockoffs method has a median value different from 0, leading to an estimated FDR at 0.08. For the scale-free-min, the methods generally select the empty subset and if it is not the case, the selected subset contains non-active variables leading to a high estimated FDR.

To conclude, in the independent setting, eBIC has a reasonable estimated FDR and the knockoffs method gets an estimated FDR close to the threshold provided as an input of the method. When a dependency structure exists, eBIC, Bolasso, Tigress en the knockoffs method should be preferred. Conversely, the data-driven penalty strategies always provide high median values. Because of the variability of the results, the ℓ 1 regularization is preferred over the E-Net one.

Results from the FRANK designs

This subsection is devoted to analyze performances of the methods from the FRANK data. These datasets are interesting for an application point of view since the dependency structure of the data is close to a gene regulatory network. A complementary study shows that the distribution of the FRANK data is far from the Gaussian distribution. Thus, the FRANK data deviate from the statistical model assumptions and this subsection investigates the impact of the non-Gaussian assumption. We recall that we consider the variables having either the highest (setting denoted by FRANKmax) or the smallest (setting denoted by FRANK-min) number of neighbors as response variables. We only present the conclusions. Graphical representations are available at https://sites.google.com/view/placroix/research.

Area under the ROC curve (pROC-AUC):

Similarly with scenarios from independent and Gaussian models, the combination of the E-Net regularization with the lars algorithm achieves the highest value of pROC-AUC. However, these values equal 0.5 for FRANK-max and 0.57 for FRANK-min, significantly smaller values compared to other scenarios. Median values from other combinations are smaller than 0.3.

To conclude, the quality of the obtained regularization paths are clearly deteriorated on FRANK data.

MSE:

For FRANK-max design, all the model selection methods lead to a median value larger than 1, meaning that they are not predictive. As for the variable identification methods, the median values are strictly smaller than 1 but strictly larger than 0.99 except for ESCV with values between 0.97 and 0.99.

For FRANK-min design, among the model selection methods, eBIC and LinSelect are the best methods, with values between 0.97 and 0.99. Those from the data-driven penalties are always larger than 1. ESCV is the best method among the variable identification ones, with median value between 0.97 and 0.99. The others provide a median value between 0.99 and 1.

To conclude, the prediction performances are deteriorated on the FRANK designs compared to the independent and Gaussian models. The model selection methods lose even their prediction performances.

Recall:

For FRANK-max design, all the median values are null expect for ESCV, the dimension jump and slope heuristics with values respectively smaller than 0.06, 0.08 and 0.24.

For FRANK-min design, values are null for all the methods.

To conclude, sensibility is drastically deteriorated on FRANK data meaning that all the methods fail to select the active variables.

Specificity:

For FRANK-max design, concerning the model selection methods, median values for eBIC and LinSelect are larger than 0.99 while those for the data-driven penalties decrease until 0.93 and 0.7 for respectively the dimension jump and the slope heuristics strategies. As for the identification methods, the median values equal 1 except for ESCV where values are larger than 0.96 and Bolasso with lars, E-Net and when samples are first generated where values are larger than 0.98.

For FRANK-min design, the same conclusions are observed except for the slope heuristics where median values decrease until 0.84.

To conclude, the specificity values remain high on FRANK data for all methods and the data-driven penalty strategies are the worst ones.

FDP and FDR:

For FRANK-max design, concerning the model selection methods, the median values equal 1 for eBIC and LinSelect, meaning that the selected variables are all non-active. Median values equal 0.9 and 0.88 for the slope heuristics and the dimension jump respectively. Concerning the identification methods, median values equal 0.8 for ESCV and Stability Selection with lars, E-Net and when samples are first generated, and 0.6 for Bolasso with lars, E-Net and when samples are first generated. All the other combinations provide a null median value.

For FRANK-min design, all the model selection methods achieve 1. As for the variable identification methods, all the median values are null expect for ESCV, Bolasso and Stability Selection with lars, E-Net and when samples are first generated with values equal 1, 0.7 and 1 respectively.

To conclude, the FDR values are drastically deteriorated on FRANK data with values close to 1 in most cases.

Impact of the high-dimension

This subsection looks at the impact of the high-dimension. If the results of a method remain similar whereas the sample size n increases, then it ensures that the high-dimension is not the reason of poor performances. For a neutral comparison, results of this subsection are obtained from the same datasets used from Subsections 4.1 to 4.6 that we complete to get samples of size n = 300, 600 and 1200. We only present the conclusions. Graphical representations are available at https://sites.google.com/view/placroix/research.

Area under the ROC curve (pROC-AUC)

In the independent setting, the pROC-AUC equals 1 for lars when n increases. For the GD algorithm, the result is less expected since the pROC-AUC increases from 0.5 to 0.9 when n goes from 150 to 300 and then decreases up to 0.76 when n equals 600 or 1200. The difference between the two regularization functions with lars algorithm that we observed at n = 150 does not exist when n increases.

When a dependency structure exists, lars combined with E-Net remains the best but the differences between the four combinations between an algorithm and a regularization function become negligible.

MSE

For the independent design, prediction performances of LinSelect improve drastically since values are below 0.1 when n = 600. When n = 1200, all the model selection methods provide median values smaller than 0.3. For the variable identification methods, the median values decrease with n except for Bolasso and Stability Selection when the grids are generated first. For n = 300, the knockoffs method has a median value close to 1 and then values decrease with n. As soon as n = 600, all the methods are similar with median values lower than 0.2 except for Tigress which still remains less predictive. We observe that the combination between lars and E-Net is recommended for Bolasso and Stability Selection.

When a dependency structure exists, approximately no change is observed for cluster and scale-free-min designs. Slight decreases in values are observed but median values remain close to 1. For scale-free-max design, increasing the value of n changes the best method which become eBIC instead of the data-driven penalty strategy with the dimension jump strategy. All model selection methods give values between 0.2 and 0.3 when n = 600. Concerning the variable identification methods, only Bolasso and Stability Selection used with lars and E-Net show a decreasing with n up to a median value indicating predictive performances (at most 0.4). Tigress and the knockoffs remain with high values (respectively 0.8 and around 1).

Recall:

For the independent design, all the model selection methods have a recall equals 1 as soon as n = 600. Among the variable identification methods, it is not verified for Tigress and for Bolasso and Stability Selection where values are between 0.7 and 1 when samples are first generated but lower than 0.5 when grids are first generated. Of note that the knockoffs methods show an instability with respect to n.

When a dependency structure exists, the recall is generally higher for the model selection methods than for the variable identification methods. Linselect does not manage properly the recall in the cluster design. Concerning the variable identification methods, they behave differently in designs reflecting an instability: the knockoffs method remains at 0 until n = 1200 for cluster and scale-free-min but finishes at 1 for the scale-free-max; ESCV values achieve 1 for scale-free-min, 0.9 for scale-free-max with E-Net but 0.1 for the ℓ 1 regularization, and 0.3 for cluster. As for Bolasso and Stability Selection, median values are 0.9 when grids are first generated but only 0.6 when samples are first generated. Lastly, Tigress is not sensitive, with values always lower than 0.2.

Specificity

In the independent design, only the data-driven penalty strategies do not have a high median value when n = 150, but when n increases, 1 is achieved. Of note, the LinSelect performances are deteriorated as soon as n = 600 and are unstable according to the choice of algorithm and regularization function. For the variable identification methods, the median values remain at 1 for all the methods when n increases.

When a dependency structure exists, the GD algorithm is better than the lars one. We observe that the specificity decreases when n increases for the data-driven penalties. One possible explanation is that the size of the selected subsets tend to be larger that the set of the active variables when n increases and consequently some non-active variables are selected. As soon as n = 600, eBIC has a decreasing specificity for the cluster design, even if the median size of the estimated support is smaller than the true support. Linselect has also a decreasing specificity in the scale-free-max design as soon as n = 600 and the explanation is the number of selected variables which is really to high compared to the number of active variables. For the variable identification methods, the specificity is always 1 except for Bolasso and Stability Selection in the cluster and scale-free-min designs and for the knockoffs method when n = 1200 in the scale-free-max design.

FDP and FDR

In the independent design, the median value decreases when n increases for all the methods except for Linselect which selects too much variables. The median values of the variable identification methods tend to 0 when n increases.

When a dependency structure exists, Linselect has surprisingly a median value increasing with n expect for the scalefree-min design. This is also the case for the data-driven penalties. eBIC is the model selection method having the best estimated FDR whatever the designs. The median values of the variable identification methods remain very low when n increases, surprisingly except for Stability Selection in scale-free-min and cluster designs.

Take home message

As preliminary recommendations, Tables 2 and3 suggest that the E-Net regularization, the lars algorithm and the strategy of samples first generated have to be privileged to get a size of selected variables set close to the size of the set of active variables. The considered data-driven penalties provide a too large set of selected variables. Conversely, Tigress and LinSelect are both conservative leading to an almost empty set of selected variables.

Recommendation per method

The structure of dependencies between variables in a dataset is often unknown. Therefore, we summarize in the following table the best combination of the algorithm, the regularization function and eventually the strategy of the procedure for each of the variable selection method, regardless of the dependency between the variables:

method

Discussion

High-dimensional Gaussian linear regression is commonly used to model interactions between entities when the number variables is close to or larger than the number of observations. This framework raises many methodological questions and this review aims at highlighting the performances of some methods according to some different objectives. This work is focused on variables selection methods from regularization path constructions. We propose an evaluation of both the regularization path construction and the choice of the final selected variables. The first step is based on the minimization of the least-squares penalized by a regularization function and the main question is the choice of the regularization function. The second step is based on selection of the final variable subset either by penalized criterion minimization with either asymptotic or non-asymptotic properties (model selection methods), or by data sampling strategies (variable identification methods). To evaluate the different methods in a fair way, we simulated different settings, each one having its own characteristics. The independent setting is irrelevant for most of the applications since having completely independent variables is rare in practice. However, this setting is commonly used to develop the statistical methods and is a benchmark to evaluate performances of any method. The settings based on the Gaussian graphical model generating correlated variables are the most favorable case where assumptions of our model broadly hold. The FRANK setting provides a more realistic framework (it is based on a dynamic process to generate a gene regulatory network) but deviates from the statistical model's assumptions, especially the Gaussian distribution. Lastly, the methods are evaluated for different performances: the ability to discriminate the active variables from the others through the pROC-AUC, the prediction performances through the MSE, the ability to recover the active variables through the recall, the ability to not select the non-active variables through the specificity and the quality of the selected variables subset through the FDR.

The impact of dependency structure between variables is evaluated by comparing results on the independent case with those on the three other designs. Our results show significant degradation of performances for all the metrics on the three last settings in comparison with the independent one. This proves that controlling the variable dependency structure is an important assumption for the statistical procedure. Our observations suggest that when dependence between variables exists, the methods work better when the support is large enough (scale-free-max).

One of the most striking conclusions is that choosing the best method depends on the metric to control. Analyses are performed by studying the medians and the variability of the results metric by metric. The first step of the statistical framework is to order the variables through regularization path constructions. According to the simulation study, the combination of the lars algorithm and the E-Net regularization is the best one to discriminate the active variables and the non-astive ones. In a complementary study, we observe that some non-active variables appear at the beginning of the regularization paths and their number is reduced when n increases. So, the high-dimensional context impacts the variable selection procedure at the step of the regularization paths constructions. To define a final variable subset from the regularization paths, ESCV and the knockoffs are the best methods for prediction performances. This is surprising since model selection procedures are constructed for theoretical guarantees on the predictive risk. Among them, eBIC and LinSelect are competitive and provide slightly higher values of the MSE compared with ESCV and the knockoffs. To recover all active variables, all the model selection methods excepted LinSelect are to be favored: they provide a high recall. ESCV is the best method among the variable identification ones. All variable identification methods are specific and eBIC is the best method among the model selection ones for specificity. So, to control both recall and specificity, ESCV and eBIC are the best choices. Lastly, to ensure that the selected variables are active variables, a small value of FDR is expected. In this direction, Bolasso and the knockoffs get a FDR value smaller than 0.1, the common threshold used in practice. Tigress and LinSelect are conservative providing often an empty set of variables.

As general conclusions, while methods based on asymptotic properties are commonly used in literature, this review emphasizes that the non-asymptotic ones should also be considered. We often observe bad performances with the data-driven penalties. One reason may be that data-driven methods are based on a heuristic whereas LinSelect, the other non-asymptotic model selection procedure, was constructed from an oracle inequality. Moreover, we also want to point out that for the data-driven penalty methods, the penalty shape and the multiplicative constant 2.5 in ( 6) are derived from a variable selection procedure in a context of changepoint detection [START_REF] Lebarbier | Detecting multiple change-points in the mean of gaussian process by model selection[END_REF], so in a different statistical framework from the one we are studying in this review. These quantities may be not adapted in a high-dimensional Gaussian linear regression and more work is required to propose other calibrations adapted in high-dimensional regression. More generally, the choice of parameters (setting to default values) for each method can also be discussed.

A statistical framework always requires some assumptions. For the high-dimensional Gaussian linear regression, data are assumed to be distributed according to a Gaussian distribution, observations are supposed to be independent, variable dependencies are well controlled and the set of active variables is supposed to be small. Experimental datasets generally do not verify all these assumptions. In this work, we just investigate the impact of the high-dimension and the relaxation of the Gaussian assumption but some other relaxations can be studied in future works. We observe that the MSE continue to be high for the data-driven penalties and Tigress, as well as the recall remaining small, even if n increases. So, high-dimension is not the reason to explain these bad performances. Surprisingly, FDR values drastically increase with n for LinSelect, eBIC and Stability Selection and the specificity values decrease for LinSelect with n. Since our study shows deterioration of values for all metrics in the FRANK data, the Gaussian distribution seems to be an important assumption for the considered statistical model. Conclusions could indicate the importance of data transformations and preprocessing steps if the Gaussian assumption is not satisfied (see [START_REF] Liu | High-dimensional semiparametric gaussian copula graphical models[END_REF]). Of note, in a complementary study, the shrinkage transformation (available in the R function huge.generator from the R package huge) is tested on FRANK-max and FRANK-min to replace the classical normalization per variable. The MSE values are slightly improved but all the other metrics are slightly deteriorated (boxplots are provided in supplementary material available in 1 ).

Supplementary data

The scripts are available from 2 . Boxplots of results in Subsection 4.8 and 4.7 are provided in supplementary material available in3 .

8 Appendix: Boxplots for n = 150 per metric In this appendix, table of the number of active variables per simulation setting and table of the number of selected variables per method and per simulation setting are provided when n = 150. Then, boxplots of results from scenarios independent, cluster, scale-free-max and scale-free-min when n = 150 are provided. They are arranged by metric. For boxplots of results when n = 300, 600, 1200 as well as for boxplots of results from the FRANK scenarios, the reader is invited to consult the supplementary material available in 4 .

The glm, enet and lasso denote respectively the gradient descent algorithm, the Elastic-Net regularization function and the Lasso regularization function. The slope heuristics method and the dimension jump are named respectively by slope and jump. Bolasso and Stability Selection are named respectively by b and ss. Lastly, grid and sub denote respectively the strategy when grids are first generated and the strategy when samples are first generated.

Figures 4 and 5

 5 Figures 4 and 5 summary the recall values of the different methods. The target median value is 1.

Figure 2 :

 2 Figure 2: Boxplots of the MSE values for model selection procedures and for n = 150.

Figure 3 :

 3 Figure 3: Boxplots of the MSE values for variable identification procedures and for n = 150.

Figure 5 :

 5 Figure 5: Boxplots of the recall values for variable identification procedures and for n = 150.

Figure 7 :

 7 Figure 7: Boxplots of the specificity values for variable identification procedures and for n = 150.

  More precisely, the data-driven penalties have to be avoided due to bad performances and instabilities in results. Among the variable identification methods, we suggest to use Bolasso instead of Stability Selection since (i) Stability ESCV recall ones, which makes eBIC and ESCV competitive methods for control of both recall and specificity.In the following Table are summarizing the best and the worst methods to provide variable sets with a small number of non-active variables (small value of FDR): Bolasso and the knockoffs have FDP values smaller than 0.1, the threshold usually fixed in practice. LinSelect and Tigress are set apart since they often select no variable.

	FDR:
	Best methods:		Bolasso, the knockoffs
	Best regularization function:	ℓ 1
	Best algorithm:		GD, lars
	methods to avoid:		ESCV, data-driven penalties
	More precisely,
		regularization function	algorithm	strategy
	path:	E-Net	lars
	eBIC:	ℓ 1	lars
	data-driven penalties:	to avoid
	LinSelect:	indifferent	indifferent
	ESCV:	indifferent	indifferent
	The knockoffs:	indifferent	indifferent
	Tigress:	ℓ 1 (per default)	lars (per default)
	Bolasso:	E-Net	lars	samples first generated
	Stability Selection:	to avoid

Table 1 :

 1 Active variable number

	Active variable number	ind	cluster	scale_free_max	scale_free_min
	mean	12.59	11.63	31.41	1
	(sd)	(1.76)	(2.75)	(9.70)	(0)
	Selected variables number	ind	cluster	scale_free_max	scale_free_min
	GD_E-Net_ebic (mean)	14.34	1.27	7.53	1.03
	(sd)	(2.84)	(0.63)	(2.30)	(0.17)
	GD_E-Net_slope (mean)	36.80	62.69	60.63	57.99
	(sd)	(16.64)	(34.98)	(21.47)	(36.28)
	GD_E-Net_jump (mean)	25.83	23.37	42.20	24.38
	(sd)	(12.23)	(17.13)	(20.19)	(19.70)
	GD_E-Net_linselect (mean)	1.07	1.32	1.59	1.04
	(sd)	(0.33)	(1.12)	(1.51)	(0.20)
	GD_ℓ 1 _ebic (mean)	13.50	1.27	7.67	1.03
	(sd)	(2.21)	(0.63)	(2.40)	(0.17)
	GD_ℓ 1 _slope (mean)	35.69	62.94	62.61	59.67
	(sd)	(16.12)	(35.49)	(21.30)	(35.06)
	GD_ℓ 1 _jump (mean)	28.23	23.85	42.49	24.12
	(sd)	(11.93)	(17.61)	(20.16)	(18.67)
	GD_ℓ 1 _linselect (mean)	1.07	1.30	1.53	1.04
	(sd)	(0.33)	(1.02)	(1.49)	(0.20)
	lars_E-Net_ebic (mean)	14.80	1.25	7.43	1.02
	(sd)	(3.76)	(0.61)	(2.43)	(0.14)
	lars_E-Net_slope (mean)	20.81	23.15	45.33	22.09
	(sd)	(10.43)	(15.93)	(17.30)	(14.17)
	lars_E-Net_jump (mean)	23.47	23.34	43.78	22.86
	(sd)	(13.39)	(13.89)	(15.88)	(13.69)
	lars_E-Net_linselect (mean)	1.04	1.23	1.55	1.03
	(sd)	(0.28)	(0.76)	(1.51)	(0.17)
	lars_ℓ 1 _ebic (mean)	13.50	1.26	7.71	1.02
	(sd)	(2.21)	(0.61)	(2.36)	(0.14)
	lars_ℓ 1 _slope (mean)	61.65	54.34	67.72	53.51
	(sd)	(29.42)	(32.72)	(24.92)	(33.65)
	lars_ℓ 1 _jump (mean)	32.30	25.13	44.72	23.84
	(sd)	(15.33)	(19.54)	(22.14)	(18.35)
	lars_ℓ 1 _linselect (mean)	1.04	1.26	1.46	1.03
	(sd)	(0.28)	(0.81)	(1.44)	(0.17)

Table 2 :

 2 Selected variables number for model selection methods

	Selected variables number	ind	cluster	scale_free_max	scale_free_min
	GD_bol_grid_E-Net (mean)	4.23	0.20	1.00	0.01
	(sd)	(1.28)	(0.45)	(0.86)	(0.10)
	GD_bol_grid_ℓ 1 (mean)	4.08	0.10	0.46	0.00
	(sd)	(1.32)	(0.36)	(0.59)	(0.00)
	GD_bol_sub_E-Net (mean)	2.48	0.17	0.42	0.01
	(sd)	(1.02)	(0.43)	(0.57)	(0.10)
	GD_bol_sub_ℓ 1 (mean)	2.34	0.08	0.23	0.00
	(sd)	(0.99)	(0.31)	(0.42)	(0.00)
	GD_escv_E-Net (mean)	11.61	6.81	22.42	4.64
	(sd)	(4.42)	(8.71)	(11.18)	(10.37)
	GD_escv_ℓ 1 (mean)	11.59	6.86	22.42	5.07
	(sd)	(4.41)	(8.62)	(11.13)	(10.46)
	GD_knockoffs_E-Net (mean)	12.11	0.00	25.74	0.00
	(sd)	(4.68)	(0.00)	(8.44)	(0.00)
	GD_knockoffs_ℓ 1 (mean)	12.11	0.00	25.74	0.00
	(sd)	(4.68)	(0.00)	(8.44)	(0.00)
	GD_ss_grid_E-Net (mean)	6.18	0.49	2.84	0.19
	(sd)	(1.74)	(0.70)	(1.27)	(0.42)
	GD_ss_grid_ℓ 1 (mean)	6.16	0.33	1.73	0.09
	(sd)	(1.70)	(0.59)	(1.10)	(0.29)
	GD_ss_sub_E-Net (mean)	3.89	0.34	1.64	0.15
	(sd)	(1.23)	(0.59)	(1.05)	(0.39)
	GD_ss_sub_ℓ 1 (mean)	3.81	0.26	1.04	0.07
	(sd)	(1.20)	(0.50)	(0.90)	(0.26)
	lars_bol_sub_E-Net (mean)	7.15	1.48	9.19	0.43
	(sd)	(1.35)	(1.65)	(2.21)	(0.54)
	lars_bol_sub_ℓ 1 (mean)	0.00	0.00	0.04	0.00
	(sd)	(0.00)	(0.00)	(0.20)	(0.00)
	lars_ss_sub_E-Net (mean)	9.11	2.79	12.64	0.99
	(sd)	(1.35)	(1.97)	(2.25)	(0.80)
	lars_ss_sub_ℓ 1 (mean)	0.00	0.00	0.16	0.00
	(sd)	(0.00)	(0.00)	(0.42)	(0.00)
	lars_tigress_ℓ 1 (mean)	2.14	0.31	0.40	0.14
	(sd)	(0.89)	(0.56)	(0.59)	(0.35)

Table 3 :

 3 Selected variables number for variable identification methods

	independent

https://sites.google.com/view/placroix/research

https://forgemia.inra.fr/GNet/high-dimensional_regression_comparison

https://sites.google.com/view/placroix/research

https://sites.google.com/view/placroix/research

Acknowledgments

We are grateful to the INRAE MIGALE bioinformatics facility (MIGALE, INRAE, 2020. Migale bioinformatics Facility, doi: 10.15454/1.5572390655343293E12) for providing computing resources. IPS2 benefits from the support of the LabEx Saclay Plant Sciences-SPS (ANR-17-EUR-0007). This work was supported by a public grant as part of the Investissement d'avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH.

Selection gives similar or poorer results compared to Bolasso, (ii) Stability Selection results often depend on the variable dependency structure and are unstable, (iii) the Stability Selection method is processed on samples of size n 2 , which is even smaller than p, providing additional estimation problem due to the high-dimension.

Recommendation per metric

In this subsection, we summarize metric per metric the best methods as well as those to avoid, regardless of the dependency between the variables. GD, lars methods to avoid: Tigress, LinSelect, data-driven penalties

In the following

Recall and specificity metrics have to be considered together. Indeed, only controlling the recall would lead to a variable set containing all the active variables but also many non-active ones. Conversely, only controlling the specificity would lead to a variable set with a few non-active variables but many active variables would not be selected. Thus, only controlling specificity leads to a conservative method with a small recall, while controlling recall leads to a sensitive method with a small specificity. Controlling both metrics simultaneously should provide a variable set closer to the active variable set. Table above provides methods with high values for both metrics. More precisely, Tigress and LinSelect are set apart since their recall values are smaller than 0.5. In contrast, the data-driven penalties are set apart since their specificity values are the smallest with respect to all the others. The remaining variable identification methods have the same performances in terms of specificity metric and ESCV provides the highest recall value. In the same way, the remaining model selection methods have the same performances in terms of the recall metric and eBIC provides the highest specificity value. These values are just below those of ESCV but eBIC recall values are just above independent cluster scale-free-max scale-free-min


Igor Khatuntsev
email: khatuntsev@iki.rssi.ru

Marina Patsaeva

Dmitri Titov

Nikolai Ignatiev

Anna Fedorova

Alexander Turin

Ludmila Zasova

Jean-Loup Bertaux

Winds from visible (513 nm) images obtained by the Venus Express Monitoring Camera

Long-term observations of Venus by the Venus monitoring camera (VMC) [4] onboard Venus Express in the visible range (513 nm) allowed to obtain information on circulation in the middle of Hadley cell. Here we present preliminary results of wind tracking of visible images obtained by VMC in the period 2007/07/01-2013/01/29. The mean zonal speed in middle latitudes (30-65ºS) is monotonously changing with latitude -76.5 to -61.5 m/s. In low latitudes (10-20ºS) zonal speed is about -82 m/s. Differences in the behavior of the meridional component in low and middle latitudes indicate that different cloud layers are observed in the visible range (513 nm). The mean horizontal flow demonstrates a dependence on local solar time. Moreover, the nature of diurnal variations is different for low and middle latitudes. Following Bertaux et al. [1] the mean horizontal flow in the middle of Hadley cell show a connection with the underlying topography.

Hadley cell from VMC observations

During ~9 years (from April 2006 through January 2015) Venus Monitoring Camera (VMC) [START_REF] Markiewicz | Venus Monitoring Camera for Venus Express[END_REF] onboard the Venus Express orbiter has observed the upper cloud layer of Venus. The largest set of images was obtained in the UV (365 nm), visible (513 nm) and two infrared channels -965 nm and 1010 nm. The UV dayside images were used to study the atmospheric circulation at the Venus cloud tops [START_REF] Khatuntsev | Cloud level winds from the Venus Express Monitoring Camera imaging[END_REF], [START_REF] Patsaeva | The relationship between mesoscale circulation and cloud morphology at the upper cloud level of Venus from VMC/Venus Express[END_REF]. Cloud features tracking in the IR images (965 nm) provided information about winds in the middle and lower cloud (49-57 km) [START_REF] Khatuntsev | Winds in the middle cloud deck from the near-IR imaging by the Venus Monitoring Camera onboard Venus Express[END_REF]. Meridional speed for UV observations (cloud tops) has a negative sign, at the same time for IR observations the mean meridional speed in the middle cloud deck (55 km) has a positive sign (Figure 1). This result, together with the earlier measurements of the poleward flow at the cloud tops indicates the presence of a closed Hadley cell in the altitude range 55-67 km. In middle latitudes (30-65ºS) the zonal speed is monotonously changing from -76.5 to -61.5 m/s with South latitude increase. The meridional speed is also decreasing from 2.8 to 0.8 m/s and this parameter is systematically less than the meridional speed obtained from IR images. In low latitudes (10-20ºS) the zonal speed is dramatically higher (about -82 m/s) and looks like equatorial jet. The meridional component is also increasing from 0.4 m/s at -10ºS to 6 m/s at 20ºS with a maximum (~9 m/s) near -15ºS.

Winds in the middle of Hadley cell

EPSC Abstracts

Vol. 13, EPSC-DPS2019-309-2, 2019 EPSC-DPS Joint Meeting 2019 c Author(s) 2019. CC Attribution 4.0 license. Thus, VMC is observing an overlying layer in low latitudes relative to the layer of middle latitudes.

In low latitudes the mean horizontal flow demonstrates diurnal variations. Both zonal and meridional components have afternoon maximum near -15ºS latitude. The zonal speed has a weak maximum near 13.5h, the meridional speed shows a strong peak at 14-15h (Figure 2). 2) In low latitudes the mean horizontal flow shows diurnal variations.

Summary and Conclusions

3) Zonal and meridional components of wind speed in the dayside Hadley cell show a connection with the underlying topography.

  Wind tracking of visible images (513 nm) obtained by VMC in the period 2007/07/01-2013/01/29 allowed us to investigate the circulation in the middle of Hadley cell. Mean zonal and meridional profiles of winds are presented at Figure 1 (green line).

Figure 1 :

 1 Figure 1: Mean zonal (top) and meridional (bottom) wind profiles. Blue correspond to UV (365 nm)[START_REF] Khatuntsev | Cloud level winds from the Venus Express Monitoring Camera imaging[END_REF], redto IR (965 nm)[START_REF] Khatuntsev | Winds in the middle cloud deck from the near-IR imaging by the Venus Monitoring Camera onboard Venus Express[END_REF] and greento visible (513 nm). In middle latitudes (30-65ºS) the zonal speed is monotonously changing from -76.5 to -61.5 m/s with South latitude increase. The meridional speed is also decreasing from 2.8 to 0.8 m/s and this parameter is systematically less than the meridional speed obtained from IR images. In low latitudes (10-20ºS) the zonal speed is dramatically higher (about -82 m/s) and looks like equatorial jet. The meridional component is also increasing from 0.4 m/s at -10ºS to 6 m/s at 20ºS with a maximum (~9 m/s) near -15ºS.

Figure 2 :

 2 Figure 2: Diurnal variations of the mean zonal (left panel) and meridional (right panel) wind speed components obtained from visible (513 nm) images.

Figure 3

 3 Figure 3 displays contour maps of mean zonal and meridional components derived from the VMC visible imaging. As in the case of the cloud top winds [1] and the case of middle cloud deck [3], we see similar evidence of the influence of Venus topography on the zonal and meridional speed. Obviously, Aphrodite Terra (60º-200ºE) significantly affects the circulation in the middle of Hadley cell.

Figure 3 :

 3 Figure 3: Venus geographic map with zonal (upper panel) and meridional (lower panel) speed contour maps.

1 )

 1 VMC visible (513 nm) images of Venus dayside were used to derive atmospheric circulation in the middle of Hadley cell at altitudes 60-62 km.

Acknowledgements J.-L. Bertaux, A. Fedorova, I. Khatuntsev, and M. were supported by the Ministry of High Education and Science of Russian Federation grant 14.W03.31.0017. N. Ignatiev, A. Turin and L. Zasova would like to thank the program #АААА-А18-118052890092-7 of the Ministry of High Education and Science of Russian Federation for the support.