Anne Gégout-Petit
email: anne.gegout-petit@univ-lorraine.fr

Aurélie Gueudin-Muller

Clémence Karmann
email: clemence.karmann@inria.fr

The revisited knockoffs method for variable selection in L 1 -penalised regressions

Keywords: Variable selection, knockoffs, regression, Lasso penalisation

We consider the problem of variable selection in regression models. In particular, we are interested in selecting explanatory covariates linked with the response variable and we want to determine which covariates are relevant, that is which covariates are involved in the model. In this framework, we deal with L 1 -penalised regression models. To handle the choice of the penalty parameter to perform variable selection, we develop a new method based on the knockoffs idea. This revisited knockoffs method is general, suitable for a wide range of regressions with various types of response variables. Besides, it also works when the number of observations is smaller than the number of covariates and gives an order of importance of the covariates. Finally, we provide many experimental results to corroborate our method and compare it with other variable selection methods.

Introduction

Regression methods are really helpful to analyse dependencies between a variable, named the response, and one or several explanatory covariates. This is one of the reasons why they are widely used and studied in statistical analysis [START_REF] Hastie | The elements of statistical learning[END_REF]. Many models have been introduced including the well-known linear regression for a continuous response variable or logistic regression for a binary response variable. Indeed, many data sets involve this last situation such as the occurence of a disease in medicine or voting intentions in econometrics. Another type of data is nominal data (that is unordered categorical data) like housing types or food choice of predators. The situation is a bit more complicated when the response is ordered categorical (ordinal), e.g. different stages of cancer, pain scales, place ratings on Google or data collected from surveys (0: never, 10: always). Logistic regression can naturally be extended to the case where the response is nominal. For such data, many authors [START_REF] Agresti | Analysis of ordinal categorical data[END_REF][START_REF] Mccullagh | Regression models for ordinal data[END_REF][START_REF] Liu | The analysis of ordered categorical data: an overview and a survey of recent developments[END_REF][START_REF] Sai Suggala | Ordinal graphical models: A tale of two approaches[END_REF] provided models based on odds ratios such as cumulative link models, adjacent-categories logit models or continuation-ratio logit models. The choice of one of these models depends on the kind of problem. In this paper, we concentrate on a restricted but large spectrum of regression models including all regression models mentioned above.

Although prediction and interpretation provide major challenges in regression motivations, another important issue is to identify the influential explanatory covariates, that is variable selection. Selection problems often arise in many fields including biology [START_REF] Tong Tong Wu | Genome-wide association analysis by lasso penalized logistic regression[END_REF]. For example, in microarray cancer diagnosis [START_REF] Zhu | Classification of gene microarrays by penalized logistic regression[END_REF], a primary goal is to understand which genes are relevant. For cost and time reasons, it can also be convenient for biologists to restrict their studies to a smaller subset of explanatory covariates (genes, bacteria populations...). Accordingly, the sparsity assumption (that is, a few number of relevant explanatory covariates) is frequently suitable and adequate, even crucial for interpretation. Indeed, with a large number of covariates, it is also useful for interpretation to determine a smaller subset of covariates that have the strongest effects. Besides, when the number of covariates is larger than the number of observations or when covariates are highly correlated, standard regression methods become inappropriate.

Lasso penalisation, or L 1 -penalisation, introduced by Tibshirani [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] offers an attractive solution to these issues. That includes a L 1 penalty in the estimation of the regression coefficients in order to perform variable selection by optimising a convex criterion. The regularisation resulting from Lasso penalisation shrinks down to zero the coefficients of explanatory covariates that have the lowest effects and leads to sparse solutions and more interpretable models, making Lasso one of the most popular penalisation [START_REF] Hastie | Statistical learning with sparsity: the lasso and generalizations[END_REF][START_REF] Zhao | On model selection consistency of Lasso[END_REF][START_REF] Young | L1-regularization path algorithm for generalized linear models[END_REF]. Moreover, using Lasso induces the critical choice of the penalty parameter which controls the number of selected covariables. This choice is major because two close values of the penalty parameter can often lead to very different scientific conclusions. Many general techniques have been proposed in the literature but they do not have the same purposes. For instance, K-fold cross validation emphasises prediction, the validation step involving computing the prediction error and aiming at minimising this. Furthermore, cross validation is often quite greedy and tends to overfit the data [START_REF] Wasserman | High dimensional variable selection[END_REF]. Other techniques, like StARS [START_REF] Liu | Stability approach to regularization selection (stars) for high dimensional graphical models[END_REF], can be adapted to a regression framework and aim at 'overselecting', that is selecting a larger set of covariates which contains the relevant ones, allowing false positive detections. Some frameworks such as gene regulatory networks require this choice: indeed, false positive detections can then be eliminated by further biological experiments whereas omitted interactions cannot be recovered after that. On the contrary, one can prefer selecting a set of covariates included in the set of true covariates to avoid false positive detections that is 'underselecting'. This constraint comes from the fact that after selection, the relevant covariates have to be studied by scientists through new experiments. But new experiments are generally expensive or time-consuming and it would be a waste to involve noisy covariates. In this paper, we focus on variable selection in the former case. Compared with our goals, we propose an intuitive and general method for automatic variable selection, inspired from the knockoffs idea of Barber and Candes [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF][START_REF] Barber | A knockoff filter for high-dimensional selective inference[END_REF]. This method uses a matrix of knockoffs of the covariates built by swapping the rows of the matrix of covariates. This knockoffs matrix is thus random and aims at determining if a covariable belongs to the model using a decision rule based on change detection methods. One of the major advantages is that it can be performed in a wide range of regression frameworks including when the number of covariates is much larger than the number of observations. We will see that our method does not lead to a choice of the penalty parameter. Nevertheless, it provides an order of importance on the covariates allowing to select covariates according to the target.

In this paper, we address the problem of variable selection in L 1 -penalised regressions. Our goal is to determine which covariates are relevant and which are noisy. We achieve it by proposing a new method of type knockoffs. The rest of the paper is organised as follows. In Section 2, we first introduce the background and describe the knockoffs method for variable selection. We also describe briefly our R package kosel in which the revisited knockoffs method is implemented. In Section 3, we give many illustrations and results of our method on simulated data. Furthermore, we propose a way to exploit randomness of our procedure.

f k (µ k (Y |X)) = α k + β 1 X 1 + . . . + β p X p , k = 1, . . . , m, (1)
where f k is a deterministic function, µ k (Y |X) parameters of the distribution of Y given X and α k , β 1 , . . . , β p real coefficients. Note that the vector of regression coefficients β := (β 1 , . . . , β p ) does not depend on k.

This framework is quite general and emcompasses many regression models such as generalised linear models [START_REF] Agresti | Categorical data analysis[END_REF] (linear regression, logistic regression, Poisson regression, multinomial regression), ordinal logistic regression models [START_REF] Agresti | Analysis of ordinal categorical data[END_REF] (cumulative logit models with proportional odds [START_REF] Simon | Alternative analyses for the singly-ordered contingency table[END_REF][START_REF] Williams | Analysis of contingency tables having ordered response categories[END_REF][START_REF] Anderson | Regression, discrimination and measurement models for ordered categorical variables[END_REF], adjacent-categories logit models, continuation-ratio logit models) or cumulative link models [START_REF] Agresti | Analysis of ordinal categorical data[END_REF]. Indeed, for generalised linear models, m = 1, µ 1 (Y |X) = E(Y |X) and f 1 is the link function (identity, log, logit...) of the corresponding model. For ordinal logistic regression models,

f k = logit and µ k (Y |X) = P(Y ≤ k|X) (cumulative), µ k (Y |X) = P(Y = k|Y = k or Y = k + 1, X) (adjacent), µ k (Y |X) = P(Y > k|Y ≤ k, X) (continuation).
Notice that these last models only allow identical effects of the covariates, which implies that the regression coefficients β i do not depend on the modality k of the response variable Y . In particular, this framework includes models for many types of response variable such as binary, continuous, ordinal or categorical.

In this framework, covariates X l , l = 1, . . . , p have to be linked to the response variable Y through a linear expression so that the conditional dependence between Y and X l given X 1 , . . . , X l-1 , X l+1 , . . . , X p can be measured through the regression coefficient β l . More precisely, β l = 0 means that X l and Y are independent conditionally on the other covariates X k , k = 1, . . . , l -1, l + 1, . . . p. We are thus interested in the nullity of the regression coefficients β l to select the relevant covariates. Moreover, we make sparsity assumption, that is a relatively small number of covariates play an important role. This implies that only a few covariates are relevant and thus, only a few regression coefficients β l are non-null. This sparsity assumption is convenient for scientists to restrict their studies to a smaller subset of covariates, namely in high-dimensional settings. Instead of checking the nullity of each coefficient β l by performing statistical tests, we add a L 1 -penalisation on the coefficients β in the estimation of the coefficients of the model. Coefficients are usually estimated by solving the optimisation problem:

argmax (α,β) L(α, β, Y, X), (2)
where L(α, β, Y, X) is a function of the coefficients relative to the model (like log-likelihood), depending on the observations Y and X of the response variable Y and the vector of covariates -→ X respectively. Instead of estimating the coefficients by (2), we add a Lasso penalisation on the coefficients vector β which leads to solve the following optimisation problem: argmax

(α,β) L(α, β, Y, X) -λ||β|| 1 , (3)
where λ > 0 is the penalty parameter.

Usually, all penalisation methods require the choice of the (positive) penalty parameter, also referred as tuning or regularisation parameter. We then need to tune the penalty parameter λ (involved in the optimisation problem (3)) which controls the number of selected covariates: the larger λ is, the fewer the selected covariates are. On the contrary, values of λ closed to 0 lead to the full model, that is the model with all the covariates. We remind that our goal is to select only relevant covariates and thus, to avoid false positive detections (the wrongly detected covariates).

With regard to our problems and goals, we propose a new method, inspired from the knockoffs process used by Barber and Candes [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF] in the linear Gaussian regression setting. Actually, this method does not lead to a choice of the penalty parameter λ but it puts the explanatory covariates in order from the most relevant to the least. Furthermore, it suits any regression of the type presented in (1) including when the number n of observations is smaller than the number p of covariates. Obviously, in the linear Gaussian model, it is much more pertinent to use the procedure described in [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF] because of their theoretical guarantees. Even if their procedure initially held for n > p, they subsequently extended it thanks to a preliminary screening step [START_REF] Barber | A knockoff filter for high-dimensional selective inference[END_REF]. In what follows, we present the principle of our revisited knockoffs method.

Principle and generalities

Let X denote the n × p matrix of the n observations of the p-vector -→ X = (X 1 , . . . , X p ) of covariates, called the design matrix. The principle, given in [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF], is to use a matrix X of knockoffs (of the covariates X i ) whose covariance structure is similar to X but independent from Y. The goal is to determine if a covariate X i is relevant by studying if it enters the model before its knockoff Xi , that is if X i enters the model for a larger value of the penalty parameter λ. Indeed, as the knockoff matrix is independent from Y, if a covariate enters the model after its knockoff, we can rightfully suspect that this covariate does not belong to the model. We mainly differ from the method proposed by [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF] in the construction of the knockoffs. In their paper, they propose a sophisticated construction of the knockoff filter using linear algebra tools. This construction allows to control the false discovery rate (FDR) -the expected fraction of false discoveries among all discoveries-in the linear Gaussian model whenever there are at least as many observations as covariates. This difference in the construction of the knockoffs makes our method suitable for the setting n < p and for a larger set of regression models. Nevertheless, theoretical guarantees about the control of the false discovery rate do not hold anymore. We construct our knockoff matrix X by randomly swapping the n rows of the design matrix X. This way, the correlations between the knockoffs remain the same as the original variables but the knockoffs are not linked to the response Y. Note that this construction of the knockoffs matrix also makes the procedure random. Then, in the same way as [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF], we perform the regression of Y on the n × 2p augmented matrix [X, X], i.e. the columnwise concatenation of X and X. Let us note β(λ) the estimated regression coefficients of the λ-penalised regression of Y on the augmented matrix [X, X]:

α(λ), β(λ) := argmax (α,β) L(α, β, Y, [X, X]) -λ||β|| 1 .
For each variable of the augmented design, that is for each covariate and its corresponding knockoff, we consider T i := sup {λ > 0, βi (λ) = 0}, i ∈ {1, . . . , p, p + 1, . . . , 2p}. Statistics T i correspond to the largest value of λ for which the covariate X i if i ∈ {1, . . . , p} or its knockoff Xi-p if i ∈ {p + 1, . . . , 2p} first enters the model. At this stage, we hope that T i is large for the relevant covariate, that is for X i , i ∈ {1, . . . , p} such that β i = 0 and small for the knockoffs variables X i := Xi-p , i ∈ {p + 1, . . . , 2p} or for the noisy covariates X i , i ∈ {1, . . . , p} such that β i = 0. This yields us a 2p-vector (T 1 , ..., T p , T1 , ..., Tp ) where Ti denotes T i+p . Then, we consider, for all i ∈ {1, ..., p},

W i := max(T i , Ti ) × (+1) if T i > Ti (-1) if T i ≤ Ti .
Statistics W i allo to determine if a covariate enters the model before or after its knockoff. A negative value for W i means that the covariate X i enters the model after its knockoff and we eliminate it. On the contrary, a positive value for W i means that the covariate X i enters the model before its knockoff and is more likely to belong to the model. But covariates X i whose statistic W i is positive are not necessarily relevant: we hope that W i is large for most of relevant covariates and small for the other ones. Thus, we are interested in the largest positive values of the p-vector of statistics W which moreover indicates that the corresponding covariate enters the model early, that is for a large value of λ. Statistics W i allows in fact to sort the covariates according to their importance: the larger W i is, the more relevant the associated covariate X i is. This then implies defining a threshold s for W i over which we will keep the corresponding covariates in the estimated model. On the whole, we will choose the estimated model Ŝ such that: Ŝ := {X i : W i ≥ s}.

Choice of the threshold

The second major difference with Barber and Candes [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF] lies in the choice of the threshold s. They provide in fact a data-dependent threshold that shows attractive results relative to the false discovery rate in the Gaussian setting. Unfortunately, these results do not hold in general, out of the linear Gaussian model. In our method, we make the assumption that there is a breakdown in the distributions between the W i corresponding to the covariates X i belonging to the model and the other ones (see Figure 1). Figure 1 illustrates that distributions of W i depend on whether X i is relevant or not. To generate Figure 1, we have simulated a set of data under a linear Gaussian regression model with p = 20 independent Gaussian covariates. Only the five first ones were linked to Y :

Y = β 1 X 1 + • • • + β p X p + ,
where β = (1, 1, 1, 1, 1, 0, . . . , 0) and ∼ N (0, 1). In our knockoffs procedure applicated to this data set, only statistics W 1 , W

β = (1, 1, 1, 1, 1, 0, . . . , 0)).
W i , i = 1, . . . , 20, which implies that X 3 is the covariate the most likely to belong to the model. We can clearly observe a breakdown between the values of the five first covariates and the other ones.

Consequently, we present two automatic ways to define the threshold s by using two change detection methods: the method proposed by Auger and Lawrence [START_REF] Auger | Algorithms for the optimal identification of segment neighborhoods[END_REF][START_REF] Picard | A statistical approach for CGH microarray data analysis[END_REF][START_REF] Picard | A segmentation/clustering model for the analysis of array cgh data[END_REF] and the CUSUM method for mean change detection. Let W (i) , i = 1, . . . , w, denote the sorted w positive statistics W i , i = 1, ..., w, that is 0 < W (1) ≤ W (2) ≤ . . . ≤ W (w) and e j = W (j+1) -W (j) for all j = 1, . . . w -1 the w -1 gaps between these sorted statistics. Remark that w, the number of positive statistics W i , is random (w = 11 on Figure 1). We propose two automatic thresholds defined as:

• the minimum of the two thresholds obtained by applying these two change detection methods directly on the statistics W (i) , i = 1, . . . , w sorted in ascending order,

• the minimum of the two thresholds obtained by applying these two change detection methods on the gaps e j , j = 1, . . . , w -1.

Let us name the first one 'W -threshold' and the second one 'gaps-threshold' for the sake of simplicity.

R package kosel

Our procedures have been implemented in a R package, called kosel (for knockoffs selection), available on the CRAN. Our package includes three functions: ko.glm, ko.ordinal and ko.sel.

The two first functions construct the knockoffs matrix and return the p-vector of statistics W for L 1 -penalised regressions models respectively implemented in the R functions glmnet and ordinalNet from the packages of the same name. glmnet emcompasses generalised linear models whereas ordinalNet includes ordinal regression models such as cumulative link, adjacent or continuation-ratio or stopping-ratio. By default, a seed is used so that the knockoffs matrix remains the same (and thus, the resulting statistics vector W ). But this can be changed with the option random = TRUE to exploit the randomness of the procedure (see Subsection 3.3 for further details).

The third function ko.sel deals with the choice of the threshold. It uses the statistics vector W obtained by one of the two other functions and returns the p-binary vector of estimation and the threshold s. Three choices are proposed: method = 'stats' and method = 'gaps' respectively correspond to the 'W -threshold' and 'gaps-threshold' while method = 'manual' allows the user to choose its own threshold. The option print = TRUE displays the positive statistics W i sorted in ascending order like in Figure 1. For method = 'manual', they are automatically displayed so that the user can choose its threshold. For the 'W -threshold' (method = 'stats') and 'gaps-threshold' (method = 'gaps'), option print = TRUE also displays an horizontal line corresponding to the threshold.

3 Simulation studies

Settings

We now describe experimental results to study the efficiency of our procedure. For that, we have performed different simulations with various regressions: linear Gaussian regression, logistic regression and cumulative logit regression (with proportional odds). Covariates -→ X are simulated as Gaussian such that E(X k ) = 0 and var(X k ) = 1 for all k = 1, . . . , p and such that X i and X j are dependent conditionally on the other covariates X k , k ∈ {1, . . . , p} \ {i, j} with probability 0.2. The design matrix X of covariates has been simulated with the R function huge.generator from the package huge, for a random graph structure. We have then simulated the observations of the response variable Y as:

Y = β 1 X 1 + . . . + β p X p + , (linear regression) logit(P(Y = 1|X)) = α 1 + β 1 X 1 + . . . + β p X p , (logistic regression) or logit(P(Y ≤ k|X)) = α k + β 1 X 1 + . . . + β p X p , k = 1, 2, (cumulative logit regression)
where ∼ N (0, 1) is Gaussian noise, the vector of regression coefficients β is sparse and given below and intercepts α 1 and α 2 are chosen so that the response variable Y takes enough values in each of its modalities ({0, 1} for logistic regression and {0, 1, 2} for cumulative logit regression). These regressions have been respectively performed with the R functions glmnet and ordinalNet of the eponymous R packages.

We present detection rates of each covariates on B = 100 repetitions for different settings. The detection rate of the covariate X l is the number of times among the 100 repetitions where the estimated model included X 1 . First, we have simulated n = 200 observations of p = 50 covariates for pedagogical reasons and next, n = 1000 observations of p = 2000 covariates to illustrate results in a higher-dimensional framework. For p = 2000 covariates, results are presented as boxplots of detection rates according to the regression coefficient β in order to improve readability. In addition, we compare our results with results obtained by cross validation. Cross validation has been performed with the R functions cv.glmnet for linear and logistic regressions and ordinalNetTune with 'logLik' tune method for cumulative logit regression. For p = 50, we also compare our results in the linear Gaussian setting with results obtained by Barber and Candes' procedure. Their procedure is implemented in the function knockoff.filter from the R package knockoff. We do not perform this comparison for p = 2000 because their procedure is not applicable due to the too few observations.

Efficiency and comparisons

p = 50

In the first place, we present results for n = 200 observations of p = 50 covariates. Regression coefficients are set to β = (1, 1, 1, 1, 1, 0, . . . , 0) and β = (2.5, 2, 1.5, 1, 0.5, 0, . . . , 0). Covariates X are the same for each different regression. But they are different according to the regression coefficients β. In other words, for a fixed value of β, the design matrix X is the same for each of the three regression models. But the response variable Y is simulated according to the regression model and is therefore different. The knockoffs matrix is also different. covariates and the rest of them is really clear, regardless of the regression coefficients or the choice of the threshold. For linear regression, these two thresholds give similar results whereas for logistic and cumulative logit regressions, gaps-threshold tends to give slightly lower detection rates than W -threshold, for both relevant and noisy covariates. In comparison, cross validation provides considerably higher detection rates: although the first five covariates, especially X 4 and X 5 , can be more detected, noisy covariates are also much more detected than with our procedure. For example, for logistic regression in Figure 3, noisy covariates are almost all detected less than 20% with our procedures whereas they are detected between 20% and 40% with cross validation. In practice, using cross validation would give more false positive detections than our procedures. Figure 2 also show detection rates obtained by Barber and Candes' knockoffs in the linear Gaussian regression model. To perform their procedure, we have to choose a target false discovery rate. In practice, we want it to be small but too small values lead to an infinite threshold and thus an empty estimated model. By default, the FDR is set to 0.1 but we set it to 0.4 to avoid too many empty estimated models. For the two different configurations of β, we have obtained 4 empty estimated models on 100 repetitions. Because of that, detection rates of the noisy covariates tend to be a bit higher than ours. For the same reason, the five first ones are a bit less detected (close to 96%, which corresponds to the number of repetitions for which the threshold was not infinite). For β = (2.5, 2, 1.5, 1, 0.5, 0, . . . , 0), X 5 is yet better detected. All of these three figures illustrate also that detection rates depend on the regression coefficient β: the higher β is, the more the associated covariate is detected. Indeed, for β = (2.5, 2, 1.5, 1, 0.5, 0, . . . , 0), we can observe that the covariate X 5 tends to be less detected than the four first ones. Furthermore, we can notice that some noisy covariates are more detected. For example, this is the case for the noisy covariates X 18 , X 19 , X 29 , X 36 , X 39 or X 47 for β = (1, 1, 1, 1, 1, 0, . . . , 0) and for all kind of regressions (since the design matrix is the same). This is probably due to the dependence structure of -→ X . In particular, these covariates are dependent to three of the first five covariates conditionally on the other ones. Similar phenomenon occurs for β = (2.5, 2, 1.5, 1, 0.5, 0, . . . , 0) with noisy covariates X 19 , X 41 and X 47 . Finally, our procedure seems to be quite efficient regardless of the regression model. Nevertheless, results are a little bit better for linear regression. In this case, we remind that Barber and Candes' procedures [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF][START_REF] Barber | A knockoff filter for high-dimensional selective inference[END_REF] also provide theoretical guarantees.

p = 2000

We present now results for n = 1000 observations of p = 2000 covariates to illustrate that our procedure is suitable with thousands of covariates. Regression coefficients are set to

β k =                5, if 1 ≤ k ≤ 20, 4, if 21 ≤ k ≤ 40, 3, if 41 ≤ k ≤ 60, 2, if 61 ≤ k ≤ 80, 1, if 81 ≤ k ≤ 100, 0, otherwise.
. In the same way as for p = 50 (Subsection 3.2.1), covariates X are the same for each different regression. But they are different according to the regression coefficients β. In other words, for a fixed value of β, the design matrix X is the same for each of the three regression models. But the response variable Y is simulated according to the regression model and is therefore different. The knockoffs matrix is also different.

Results and comments. Figures 5, 6 and 7 each contain four graphics: three of them are boxplots of detection rate of each of the 6 groups of covariates according to their regression coefficient β. These detection rates are respectively obtained with the revisited knockoffs methods and cross validation. In order to compare our method with cross validation for the noisy covariates, we present detection rates of the noisy covariates (these for which β = 0) obtained with the knockoffs method (with gaps-threshold) in function of detection rates obtained with cross validation in the last graphic.

Results in the linear regression framework are presented in Figure 5. Comparing the three boxplots, we can remark that detection rates with revisited knockoffs method with W -threshold are lower than with gaps-threshold. More specifically, relevant covariates whose β = 1 are detected between 10 and 95% with W -threshold whereas they are detected more than 80% with gaps-threshold. Detection rates with W -threshold are lower and more widespread for covariates whose regression coefficient β = 1 comparing to the other relevant covariates. Cross validation leads to better detection rates for the relevant covariates. However, Figure 5d illustrates that most of the noisy covariates have higher detection rates with cross validation than with our procedure. Thus, cross validation gives more false positive detections.

Results for logistic and cumulative logit regressions are respectively presented in Figures 6 and7. As for p = 50, we can note on boxplots that detection rates depend on the regression coefficient β: for all the three methods, detection rates are decreasing according to β that is, the higher β is, the more the associated covariates are detected. We also observe this on graphics 5a and 5b for linear regression, although it is less pronounced. As for linear regression, even if cross validation gives better detection rates for the relevant covariates, it also gives more fase positive detections for the noisy covariates as illustrated in graphics 6d and 7d. This phenomenon is even stronger for these two regression models where almost all of the noisy covariates are more detected with cross validation than with our procedure. Contrary to linear regression, detection rates obtained by the knockoffs method with W -threshold are higher than with gaps-threshold and they are higher for both relevant and noisy covariates. Although detection rates are better for linear regression, our procedures lead to satisfying results for the three regression models. Even though relevant covariates are not always enough detected, detection rates of noisy covariates are also often very low, especially in comparison with cross validation. On the whole, it seems to be appropriate for sparse models regardless to the regression model and particularly when the goal is to avoid false positive detections. Notice that the threshold (W or gaps-threshold) to be used to avoid false positive detections may vary according to the regression model.

Randomness of the procedure

Note that revisited knockoffs procedure and cross validation are both random (which is not the case for Barber and Candes' procedure). Indeed, the former is random in the construction of the knockoffs matrix whereas randomness of cross validation lies in the choice of the folds. Hence, applying several time one of these methods leads to different results.

To conclude this section, we compare detection rates obtained by these three methods on the same sample of data. This sample includes n = 200 observations of p = 50 covariates and consists in the first sample of the B = 100 samples used in the Subsection 3.2.1. Dependence structure of the vector -→ X is thus also the same as in Subsection 3.2.1.

Results and comments. Figure 8 displays detection rates of each covariate using randomness of the three procedures: revisited knockoffs method with W and gaps-thresholds and cross validation. Detection rates are obtained on 100 repetitions on the same sample for the three methods and for the three regression models: linear, logistic and cumulative logit regressions. Regression coefficients are set to β = (1, 1, 1, 1, 1, 0, . . . , 0). Thus, only the first five covariates belong to the model. We can notice that these first five covariates are almost always 100% detected except with gaps-threshold for cumulative logit regression (for which they are still detected more than 95%). Noisy covariates, that is covariates X 6 , . . . , X 50 , are always less detected by our procedures than by cross validation. However, some of them are wrongly highly detected: X 18 for linear regression, X 18 , X 25 , X 30 and X 39 for logistic regression or X 19 , X 21 , X 36 and X 40 for cumulative logit regression. This is probably due to the sample (in which the dependence structure is more or less pronounced).

It should be recalled that this dependence structure is the same as in Subsubsection 3.2.1 for p = 50. In comparison with Figures 2, 3 and 4 of Subsubsection 3.2.1, we can see that these covariates have also higher detection rates. For all these covariates, cross validation gives higher detection rates than our procedures. Cross validation gives also much higher detection rates for some other noisy covariates. For instance, X 6 , X 11 , X 13 , X 24 , X 40 and X 48 are always detected for linear regression whereas revisited knockoffs with W and gaps threshold never detect them. For cumulative logit model, X 8 , X 41 , X 46 and X 50 are always detected whereas our procedures detect them less than 15%.

In practice, with real data, this randomness opens up to further ways to perform variable selection.

Discussion

In this paper, we proposed a method for variable selection in regression models based on the construction of a matrix of knockoffs of the covariates. This method is quite intuitive and suitable for many types of regressions, including when the number of observations is much smaller than the number of covariates. Two different thresholds can be chosen, leading to two procedures, which have been implementend in the R package kosel. We have seen that these procedures both turn out to be very pertinent and efficient as the many and diverse simulations exemplify. Our two procedures are particularly appropriate when the goal is to avoid false positive detections. Indeed, even if there are false negative detections, there are also a very small rate of false positive detections. Simulations show also that efficiency of our procedures depends on the regression model. In general, we can try the two thresholds and choose results according to its target. Furthermore, randomness of our procedures provides other techniques to perform variable selection. Nonetheless, in the case of linear Gaussian regression, Barber and Candes' procedures also offer theoretical guarantees. In addition, our procedures give better results than cross validation with regard to false positive detections, even when we make use of randomness. However, if we aim at overselecting, it is more appropriate to use other techniques such as cross validation.

  β = (2.5, 2, 1.5, 1, 0.5, 0, . . . , 0).

Figure 2 :

 2 Figure 2: Detection rates of each covariate for the four methods: revisited knockoffs W -threshold and gaps-thresholds, cross validation and Barber and Candes' knockoffs. Linear Gaussian regression model with n = 200 observations of p = 50 covariates with regression coefficients β = (1, 1, 1, 1, 1, 0, . . . , 0) (a) and β = (2.5, 2, 1.5, 1, 0.5, 0, . . . , 0) (b). Covariates are dependent Gaussian with a random structure. The number of i.i.d. repetitions is B = 100.

Figure 3 :

 3 Figure 3: Detection rates of each covariate for the three methods: revisited knockoffs Wthreshold and gaps-thresholds and cross validation. Logistic regression model with n = 200 observations of p = 50 covariates with regression coefficients β = (1, 1, 1, 1, 1, 0, . . . , 0) (a) and β = (2.5, 2, 1.5, 1, 0.5, 0, . . . , 0) (b). Covariates are dependent Gaussian with a random structure. The number of i.i.d. repetitions is B = 100.

  β = (2.5, 2, 1.5, 1, 0.5, 0, . . . , 0).

Figure 4 :

 4 Figure 4: Detection rates of each covariate for the three methods: revisited knockoffs Wthreshold and gaps-thresholds and cross validation. Cumulative logit model with n = 200 observations of p = 50 covariates with regression coefficients β = (1, 1, 1, 1, 1, 0, . . . , 0) (a) and β = (2.5, 2, 1.5, 1, 0.5, 0, . . . , 0) (b). Covariates are dependent Gaussian with a random structure. The number of i.i.d. repetitions is B = 100.

  Revisited knockoffs with W -threshold.

  Revisited knockoffs with gaps-threshold.

  Comparison of detection rates of noisy covariates for the revisited knockoffs method with gapsthreshold and for cross validation.

Figure 5 :

 5 Figure 5: Boxplots of detection rates of each covariate according to their regression coefficient β for the three methods: revisited knockoffs W -threshold (a), gaps-thresholds (b) and cross validation (c). Linear Gaussian regression model with n = 1000 observations of p = 2000 covariates. Covariates are dependent Gaussian with a random structure. The number of i.i.d. repetitions is B = 100.

  Revisited knockoffs with W -threshold.

  Revisited knockoffs with gaps-threshold..

  Comparison of detection rates of noisy covariates for the revisited knockoffs method with gapsthreshold and for cross validation.

Figure 6 :

 6 Figure 6: Boxplots of detection rates of each covariate according to their regression coefficient β for the three methods: revisited knockoffs W -threshold (a), gaps-thresholds (b) and cross validation (c). Logistic regression model with n = 1000 observations of p = 2000 covariates. Covariates are dependent Gaussian with a random structure. The number of i.i.d. repetitions is B = 100.

  Revisited knockoffs with W -threshold.

  Revisited knockoffs with gaps-threshold..

  Comparison of detection rates of noisy covariates for the revisited knockoffs method with gapsthreshold and for cross validation.

Figure 7 :

 7 Figure 7: Boxplots of detection rates of each covariate according to their regression coefficient β for the three methods: revisited knockoffs W -threshold (a), gaps-thresholds (b) and cross validation (c). Cumulative logit model with n = 1000 observations of p = 2000 covariates. Covariates are dependent Gaussian with a random structure. The number of i.i.d. repetitions is B = 100.

  Cumulative logit regression.

Figure 8 :

 8 Figure 8: Detection rates of each covariate for the three methods: revisited knockoffs method with the W -and gaps-thresholds (see Subsection 2.3) and cross validation. Detection rates are obtained on 100 repetitions on the same sample of n = 200 observations of p = 50 covariates with regression coefficients β = (1, 1, 1, 1, 1, 0, . . . , 0). Covariates are dependent Gaussian with a random structure.

  2 , W 3 , W 4 , W 5 , W 6 , W 7 , W 13 , W 14 , W 16 , W 19 associated to the covariates X

							W
	Sorted positive statistics W	0.2 0.4 0.6 0.8 1.0				W 2	W 4	W 1	W 5
		0.0	W 7	W 14 W 19 W 16 W 13	W 6
				2	4	6	8	10

1 , X 2 , . . . , X 14 , X 16 and X 19 have a positive value. For example, the covariate X 1 entered the model for λ = 1.002 (thus, T 1 = 1.002) and entered the model before its knockoff X1 . This means that the knockoff variable X1 entered the model for λ < 1.04 and implies that T1 < 1.002. W 3 takes the largest value among all the statistics 3 Figure 1: Example of positive statistics W i sorted in ascending order. Linear Gaussian regression model with n = 500 observations of p = 20 covariates. Only covariates X 1 , X 2 , X 3 , X 4 and X 5 belong to the model (regression coefficients are set to


Anne Gégout-Petit
email: anne.gegout-petit@univ-lorraine.fr

Aurélie Muller-Gueudin
email: aurelie.gueudin@univ-lorraine.fr

Clémence Karmann
email: clemence.karmann@inria.fr

Aurélie Gueudin-Muller

Graph estimation for Gaussian data zero-inflated by double truncation

Keywords: doubly truncated Gaussian, zero-inflation, graph estimation, precision matrix, Gaussian graphical model, graphical lasso, sparsity. maximum

come

Introduction

Multivariate data analysis often involves describing and explaining the relationships among a set of variables. Undirected graphical models offer a way to address this demand by using a graph to represent a model. A graph is a set of nodes and edges which can be represented as a graphic in order to make it easier to study, visually or computationally. Undirected graphical models are based on the conditional independence: a relation between two variables, represented by an edge in the graph, means that the corresponding variables are conditionally dependent given all the remaining variables. Among undirected graphical models, Gaussian graphical model provides a particularly convenient framework. This model assumes that the observations have a multivariate Gaussian distribution with mean µ and covariance matrix Σ. In this Gaussian setting, a direct relation between two variables corresponds to a non-zero entry in the precision matrix Σ -1 . In other words, if Σ -1 ij is zero, then variables i and j are conditionally independent given the other variables. Thus, graph estimation involves finding the zero pattern in the inverse covariance matrix. The theoretical approaches developed solve a covariance matrix whose theoretical guarantees, including the control of the convergence rate in infinite norm, required for the graphical lasso procedure, are studied.

The rest of the paper is organized as follows. In Section 2, we explicit the model and present the two step estimation procedure. The first step consists in estimating the covariance matrix, by estimating each term by maximising the corresponding bivariate marginal log-likelihood of the truncated vector, which is a non-convex optimisation problem. The second one relies on the graphical lasso to estimate the precision matrix. Section 3 contains two theoretical results about the convergence rate of the covariance and the precision matrix estimators. We first use recent results of Mei, S., Bai, Y. and Montanari, A. (2017) to set out that our covariance matrix estimator concentrates well in infinite norm around the theoretical covariance matrix. These results concern properties of the stationary points of non-convex empirical risk minimisation problems. Next, we use this first result and the consistency properties of graphical lasso, studied by Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B. (2011) in a more general framework, to show consistency and sparsistency of our final estimator of the precision matrix. In Section 3, we also state the resulting theorem which establishes consistency of our procedure with respect to graph structure recovery. In Section 4, we present some simulations studies to corroborate its theoretical efficiency.

Model and estimation procedure 2.1. Model

Let X be a Gaussian p-vector X ∼ N p (µ, Σ * ) where µ = (µ j ) j=1,...,p ∈ R p is the mean vector and Σ * = (Σ * jk ) 1≤j,k≤p ∈ M p (R) the covariance matrix. Let us consider the p-vector Y defined as: Y j = 1 aj ≤Xj ≤bj X j for all j ∈ {1, . . . , p}, where the points of truncation a j , b j ∈ R, a j < b j are known and depend on j. The Gaussian vector X is not directly observed, but it is observed through the zero-inflated truncated vector Y . Unlike what exists in the (multivariate) truncated Gaussian literature (for example, in [START_REF] Gupta | Recurrence relations for the moments of truncated multinormal distribution[END_REF], Lee, L. (1983) and Manjunath, BG and Wilhelm, S. (2009)), when one of the initial Gaussian variables falls outside its points of truncation, we observe a zero instead and we observe the rest of the vector according the same rule. In other words, our truncation does not consist in restricting to the observations of X which falls into [a

1 , b 1 ] × • • • × [a p , b p ] ⊂ R p .
Assume that µ j = 0 and Σ * jj = 1 for all j ∈ {1, . . . , p}. In practice, we can estimate them with existing techniques for doubly truncated univariate Gaussian vector (for example, (Cohen, Jr., A. C., 1957)) if variables are not centered and scaled.

Gaussian graphical model is particularly appropriate for conditional dependency graph inference. Indeed, the precision matrix Θ * := (Σ * ) -1 specifies the conditional dependency imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019 structure (see Hastie T., Tibshirani R., and Friedman J. ( 2001)). More precisely, the graph contains an edge between the variables X j and X k iff:

X j ←→ X k ⇐⇒ X j |= X k | (X l ) l =j,k ⇐⇒ cor(X j , X k | (X l ) l =j,k ) = 0 ⇐⇒ Θ * jk = 0.      (2.1)
Given a symmetric positive definite matrix M , let us denote:

E(M ) = {(j, k) ∈ {1, . . . , p} 2 , j = k, M jk = 0}. (2.2)
In particular, E(Θ * ) denotes the set of the edges of the theoretical graph.

The goal of this paper is to recover the latent graph structure of the variables of the Gaussian vector X from observations of the zero-inflated truncated vector Y .

Some theoretical tools

To explicit the model and exhibit some complexities, we give some theoretical tools and will focus here on bivariate marginal likelihood from the truncated vector Y .

Let (j, k) ∈ {1, . . . , p} 2 , j < k and let

f jk (x, y) = f (x, y, Σ * jk ) denotes the bivariate marginal log-likelihood function of the Gaussian couple (X j , X k ) ∼ N 2 0 0 , 1 Σ * jk Σ * jk 1 .
With these notations, the likelihood of (Y j , Y k ) is then L jk (Σ * jk , y) where y is an observation of the vector Y and:

L jk (σ, y) = 1 a,b=0
φ ab,jk (σ, y j , y k )n ab (y j , y k ),

(2.3) with :

• n ab (y j , y k ) = 1 ζj =a,ζ k =b where ζ l = 1 if y l ∈ [a l , b l ] \ {0}, 0 if y l = 0. • 1 a,b=0 n ab (y j , y k ) = 1 • φ 11,jk (σ, y j , y k ) = f (y j , y k , σ) • φ 01,jk (σ, y j , y k ) = φ 01,jk (σ, y k ) = [aj ,bj ] c f (x, y k , σ)dx • φ 10,jk (σ, y j , y k ) = φ 10,jk (σ, y j ) = [a k ,b k ] c f (y j , y, σ)dy • φ 00,jk (σ, y j , y k ) = φ 00,jk (σ) = [aj ,bj ] c ×[a k ,b k ] c f (x, y, σ)dxdy.
imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019

The likelihood (and log-likelihood) of a couple of variables of Y involves four terms according to the nullity of each of the components of the couple. In the same way, the likelihood of Y would involve 2 p terms by distinguishing all possible cases: the density of the Gaussian vector X (no component of Y is null), p simple integrals (only one null component), p 2 double integrals (two null components), p 3 triple integrals, ..., one p-multiple integral (all the components are null). Writing the likelihood of the vector Y becomes than quite complicated and tedious. This is why we choose to restrict to the study of the likelihoods of couples for the estimation.

In practice, we have a n-sample Y := (Y (1) , . . . , Y (n) ) of the vector Y . The likelihood of the n-sample (Y

(i) j , Y (i) k ) i=1,...,n is then L (n) jk (Σ *
jk , y) defined by:

L (n) jk (σ, y) = n i=1 L jk (σ, y (i) ), = n i=1 1 a,b=0
φ ab,jk (σ, y

(i) j , y (i) k )n ab (y (i) j , y (i) k ),
where y := (y (1) , . . . , y (n) ) is the realisation (value) of the n-sample Y. The log-likelihood is then

L (n)
jk (Σ * jk , y) where:

L (n) jk (σ, y) = n i=1 1 a,b=0
n ab (y

(i) j , y (i) k ) log φ ab,jk (σ, y (i) j , y (i) k ) = n i=1 i:y (i) j =y (i) k =0 log φ 00,jk (σ) + n i=1 i:y (i) j =0,y (i) k =0 log φ 01,jk (σ, y (i) k ) + n i=1 i:y (i) j =0,y (i) k =0 log φ 10,jk (σ, y (i) j ) + n i=1 i:y (i) j =0,y (i) k =0 log φ 11,jk (σ, y (i) j , y (i) k ) .

Estimation procedure

Our goal is to recover the latent graph structure of the Gaussian vector X, encoded in the precision matrix Θ * , from observations of the truncated vector Y . Our estimation procedure is a two step procedure. In the first instance, we estimate the covariance matrix of the Gaussian vector X. Then, we estimate the precision matrix by using the graphical lasso procedure (Friedman, J., Hastie, T. and Tibshirani, R., 2008) to recover the underlying graph structure.

Step 1: covariance matrix estimation

Estimating the covariance matrix Σ * of X as the empirical covariance matrix of the nsample Y would lead to poor results because of the zero-inflation.

Another idea could be to maximise the likelihood of the vector Y . But we have seen that this likelihood involves 2 p terms and is too tedious.

Because of these difficulties, we estimate the covariance matrix by estimating each of its entries separately using the likelihood of the couples (Y j , Y k ), j < k defined in (2.3). More precisely, we estimate Σ * by Σ (n) by estimating each of its entries Σ * jk by maximisation of the log-likelihood of the n-sample (Y

(i) j , Y (i) k ) i=1,...,n of the couple (Y j , Y k ), which is not convex. Definition 2.1 (Estimator Σ (n) of Σ * ). The estimator Σ (n) = ( Σ (n) jk ) 1≤j,k≤p of the covariance matrix Σ * is defined by: Σ (n) jk = argmax |σ|≤1 L (n) jk (σ, y) (2.4) = argmax |σ|≤1 1 n L (n) jk (σ, y),
for all j < k, where y := (y (1) , . . . , y (n) ) is the realisation of the n-sample Y.

Step 2: precision matrix estimation

As our goal is to recover the conditional dependency graph, it is natural to use the estimator of the precision matrix Θ * given by the graphical lasso (Friedman, J., Hastie, T. and Tibshirani, R., 2008). The graphical lasso is a procedure used in the Gaussian graphical model. It consists in estimating the precision matrix by maximising the penalised log-likelihood of the Gaussian model over the set p × p non-negative definite symmetric matrices: argmax

Θ 0 log det(Θ) -trace(ΘS) -λ n ||Θ|| 1,off ,
where ||Θ|| 1,off = p j,k=1 j =k |Θ jk |, S is the empirical covariance matrix of X and λ n > 0 the regularisation parameter. This optimisation problem is convex and has an unique solution (Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B., 2011).

In our case, the empirical covariance matrix of X is not directly available. Instead of obtaining this matrix as the empirical covariance matrix of Y , we replace the empirical covariance matrix S by the estimator Σ (n) of Σ * obtained at the step 1 (2.4).

Definition 2.2 (Estimator Θ (n) of Θ * ). The estimator Θ (n) of the precision matrix Θ * is defined as the unique solution of the following convex optimisation problem:

Θ (n) = argmax Θ 0 log det(Θ) -trace(Θ Σ (n) ) -λ n ||Θ|| 1,off .
(2.5)

Theoretical results of Section 3 relate the estimators Σ (n) and Θ (n) respectively defined in (2.4) and (2.5) when the points of truncation (a j ) 1≤j≤p and (b j ) 1≤j≤p are known.

Convergence results

The goal of this Section is to show that the estimation procedure proposed in Subsection 2.3 has strong theoretical guarantees. For that, we study theoretical properties of the estimator Θ (n) with regard to the recovery of the graph structure. Assume that the points of truncation (a j ) 1≤j≤p and (b j ) 1≤j≤p are known.

Covariance matrix estimator

Convergence rates in elementwise infinite norm

In a first place, we give a result about the estimator Σ (n) of the covariance matrix Σ * given by (2.4). Let us first set out two assumptions:

(H1) For all j < k, |Σ * jk | = 1. Thus, there exists δ > 0 such that for all j < k,

|Σ * jk | < 1 -δ.
(H2) Let j < k and consider the application g :

σ ∈ [-1 + δ, 1 -δ] → E L (n)
jk (σ, y) . Then, we assume that:

• -1 + δ and 1 -δ are not critical points of g,

• g has a finite number of critical points,

• every critical points of g, different from Σ * jk , are non-degenerate, i.e.:

for all σ = Σ * jk , g (σ) = 0 ⇒ g (σ) = 0.
Note that Σ * jk is a non-degenerate critical point of g. This will be proved in the proof of Proposition 3.1 (see equations (3.7)).

Proposition 3.1 states rate convergence results about the estimator Σ (n) of the covariance matrix Σ * by bounding Σ (n) -Σ * ∞ with high probability.

Proposition 3.1. Assume (H1) and (H2) and let 0 < ρ < 1. There exist some known constants B, C and D such that letting

n log n ≥ C log B ρ
, then the estimator of the covariance matrix Σ (n) defined by (2.4) satisfies:

P Σ (n) -Σ * ∞ ≥ D log n n log B ρ ≤ p(p -1) 2 ρ,
where

||A|| ∞ = max j,k∈{1,...,p}
|A jk | is the elementwise infinite norm of the matrix A.

Proof of Proposition 3.1

Proof relies on Theorem 2 of Mei, S., Bai, Y. and Montanari, A. (2017), who study the properties of the stationary points of non-convex empirical risk minimisation problems.

We begin with three auxiliary Lemmas, proved in Appendix, which all state properties about the bivariate marginal likelihood defined in (2.3) or components of it:

Lemma 3.1. There exists γ > 0 such that, for all j < k, if (y j , y k ) ∈ [a j , b j ] × [a k , b k ] and σ ∈ [-1 + δ, 1 -δ], then for all a, b ∈ {0, 1}, φ ab,jk (σ, y j , y k ) ≥ 1 γ . Lemma 3.2. There exist L 1 , L 2 and L 3 > 0 such that for all j < k, if (y j , y k ) ∈ [a j , b j ] × [a k , b k ] and σ ∈ [-1 + δ, 1 -δ], then for all a, b ∈ {0, 1}, ∂ m σ φ ab,jk (σ, y j , y k ) ≤ L m , for m ∈ {1, 2, 3}. Lemma 3.3. Let j < k. 1. For all σ ∈ [-1 + δ, 1 -δ] and for all l ∈ N * , R 2 ∂ l σ L jk (σ, y)dµ(y) = ∂ l σ R 2 L jk (σ, y)dµ(y) = 0,
where µ is the measure on R 2 defined by:

µ := δ 0 ⊗ δ 0 + δ 0 ⊗ λ + λ ⊗ δ 0 + λ ⊗ λ, (3.1)
where δ a denotes the Dirac measure in a ∈ R and λ the Lebesgue measure on R.

For all

σ ∈ [-1 + δ, 1 -δ] and for all l ∈ N * , ∂ l σ E Σ * jk log L jk (σ, Y ) = E Σ * jk ∂ l σ log L jk (σ, Y ) , i.e., ∂ l σ R 2 log L jk (σ, y)L jk (Σ * jk , y)dµ(y) = R 2 ∂ l σ log L jk (σ, y) L jk (Σ * jk , y)dµ(y).
imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019

Fix j < k. With notations of Mei, S., Bai, Y. and Montanari, A. (2017), let us set:

jk (σ, y) = log L jk (σ, y) = 1 a=0 1 b=0 n ab (y j , y k ) log φ ab,jk (σ, y j , y k ) (3.2) Rn (σ, Y) = 1 n L (n) jk (σ, Y) = 1 n n i=1 (σ, Y (i) ) (3.3) R(σ) = E Σ * jk Rn (σ, Y) = E Σ * jk (σ, Y ) . (3.4) Remarks 3.1.
• To lighten notations, we drop the underscripts jk and simply write , Rn and R instead of jk , Rn,jk and R jk . • Point 2 of Lemma 3.3 can be rewritten as:

For all σ ∈ [-1 + δ, 1 -δ] and for l ∈ N * :

R (l) (σ) = ∂ l σ E Σ * jk (σ, Y )) = E Σ * jk ∂ l σ (σ, Y )) . (3.5)
Theorem 2 of Mei, S., Bai, Y. and Montanari, A. ( 2017) requires four assumptions. Let us check these assumptions.

(i) Gradient statistical noise. The gradient of w.r.t. σ is τ 2 -sub-Gaussian for some τ > 0, i.e.:

∀σ ∈ [-1 + δ, 1 -δ], ∀λ ∈ R, E exp λ ∂ σ (σ, Y ) -E ∂ σ (σ, Y ) ≤ exp τ 2 λ 2 2 .
Indeed, for all y

∈ p j=1 [a j , b j ] and σ ∈ [-1 + δ, 1 -δ], ∂ σ (σ, y) = 1 a=0 1 b=0 n ab (y j , y k ) ∂ σ φ ab,jk (σ, y j , y k ) φ ab,jk (σ, y j , y k ) Thus: ∂ σ (σ, y) ≤ 1 a=0 1 b=0 n ab (y j , y k ) ∂ σ φ ab,jk (σ, y j , y k ) φ ab,jk (σ, y j , y k ) ≤ 1 a=0 1 b=0 n ab (y j , y k )γL 1 = γL 1 by Lemmas 3.1 and 3.2. This way, ∂ σ (σ, Y ) -E(∂ σ (σ, Y )
) is zero-mean and bounded by 2γL 1 . By Theorem 9.9 of Stromberg K. (1994)

, ∂ σ (σ, Y ) -E(∂ σ (σ, Y )) is then τ 2 -sub-Gaussian for τ = 2γL 1 . Assumption "Gradient statistical noise" is satisfied.
(ii) Hessian statistical noise. The second derivative of w.r.t. σ is τ 2 -sub-exponential (τ = 2γL 1 ), that is:

||∂ 2 σ (σ, Y )|| ψ1 ≤ τ 2 ,
where ||.|| ψ1 is the Orlicz ψ 1 -norm defined by

||X|| ψ1 := sup k≥1 1 k E X -E(X) k 1 k . For all y ∈ p j=1 [a j , b j ] and σ ∈ [-1 + δ, 1 -δ], ∂ 2 σ (σ, y) = 1 a=0 1 b=0 n ab (y j , y k ) ∂ 2 σ φ ab,jk (σ, y j , y k ) φ ab,jk (σ, y j , y k ) - ∂ σ φ ab,jk (σ, y j , y k ) φ ab,jk (σ, y j , y k ) 2 Thus: ∂ 2 σ (σ, y) ≤ 1 a=0 1 b=0 n ab (y j , y k ) ∂ 2 σ φ ab,jk (σ, y j , y k ) φ ab,jk (σ, y j , y k ) + ∂ σ φ ab,jk (σ, y j , y k ) φ ab,jk (σ, y j , y k ) 2 ≤ 1 a=0 1 b=0 n ab (y j , y k )(γL 2 + γ 2 L 2 1
) by Lemmas 3.1 and 3.2

= γL 2 + γ 2 L 2 1 . (3.6) Therefore, ∂ 2 σ (σ, Y ) -E(∂ 2 σ (σ, Y )) is bounded by 2(γL 2 + γ 2 L 2 1
) and for all k ≥ 1,

1 k E ∂ 2 σ (σ, Y ) -E(∂ 2 σ (σ, Y )) k 1 k ≤ 2 k (γL 2 + γ 2 L 2 1 ). Hence, ||∂ 2 σ (σ, Y )|| ψ1 ≤ 2(γL 2 + γ 2 L 2 1 ) ≤ τ 2 = 4γ 2 L 2 1 (we can possibly choose L 1 and γ larger). So, ∂ 2 σ (σ, Y ) is τ 2 -sub-exponential.
Assumption "Hessian statistical noise" is satisfied.

(iii) Hessian regularity.

1. The second derivative of R (defined in (3.4)) is bounded at one point:

there exists |σ * | ≤ 1 -δ and H > 0 such that R (σ * ) ≤ H.
2. The second derivative of w.r.t. σ is Lipschitz continuous with integrable Lipschitz constant (w.r.t. y), i.e.:

there exists

J * > 0 such that E[J(Y )] ≤ J * ,
where

J(y) = sup |σ1|,|σ2|≤1-δ σ1 =σ2 |∂ 2 σ (σ 1 , y) -∂ 2 σ (σ 2 , y)| |σ 1 -σ 2 | .
3. Constants H and J * satisfy: H ≤ τ 2 and J * ≤ τ 3 .

imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019

First, R (σ) = E Σ * jk ∂ 2 σ (σ, Y
) by the point 2 of Lemma 3.3 and (3.5). By (3.6),

∂ 2 σ (σ, Y ) ≤ γL 2 + γ 2 L 2 1 for all σ ∈ [-1 + δ, 1 -δ], thus any |σ * | ≤ 1 -δ and H = γL 2 + γ 2 L 2
1 are appropriate. Moreover, we have H ≤ τ 2 = 4γ 2 L 2 1 (with L 1 and γ possibly chosen larger).

For all y

∈ p j=1 [a j , b j ] and σ ∈ [-1 + δ, 1 -δ],
we have (with a slight lightening of notations):

∂ 3 σ (σ, y) = 1 a=0 1 b=0 n ab (y j , y k ) ∂ 3 σ φ ab,jk φ ab,jk -3 ∂ σ φ ab,jk ∂ 2 σ φ ab,jk φ 2 ab,jk + 2 ∂ σ φ ab,jk φ ab,jk 3
Thus: ∂ 3 σ (σ, y) ≤ 1 a=0 1 b=0 n ab (y j , y k )(γL 3 + 3γ 2 L 1 L 2 + 2γ 3 L 3 1 ) (Lemmas 3.1 and 3.2) = γL 3 + 3γ 2 L 1 L 2 2 + 2γ 3 L 3 1 .
Therefore, for all y

∈ p j=1 [a j , b j ], J(y) ≤ γL 3 + 3γ 2 L 1 L 2 2 + 2γ 3 L 3 1 by the mean value theorem. Taking J * = γL 3 + 3γ 2 L 1 L 2 2 + 2γ 3 L 3 1
with L 1 and γ possibly chosen larger, we have J * ≤ τ 3 = 8γ 3 L 3 1 . Assumption "Hessian regularity" is satisfied.

(iv) Morse. There exists > 0 and η > 0 such that R is ( , η) strongly Morse, i.e.:

1. |R (σ)| > for all σ such that |σ| = 1 -δ and, 2. for all σ such that |σ| < 1 -δ:

|R (σ)| ≤ ⇒ |R (σ)| ≥ η.
In other words, R satisfies this assumption if -1 + δ and 1 -δ are not critical points of R and if R has a finite number of critical points, which are moreover non-degenerate:

R (σ) = 0 ⇒ R (σ) = 0.
Assumption (H2) implies point 1. and point 2. for σ = Σ * jk . Let us prove that Σ * jk is a non-degenerate critical point by showing that Σ * jk is a global maximum of R. Indeed, for all σ such that |σ| < 1:

R(σ) ≤ R(Σ * jk ) ⇐⇒ E Σ * jk (σ, Y ) ≤ E Σ * jk (Σ * jk , Y ) ⇐⇒ E Σ * jk log L jk (σ, Y ) ≤ E Σ * jk log L jk (Σ * jk , Y ) since (σ, y) = log L jk (σ, y) (defined in (3.2)), ⇐⇒ E Σ * jk log L jk (σ, Y ) L jk (Σ * jk , Y ) ≤ 0.                    (3.7)
imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019

By Jensen inequality,

E Σ * jk log L jk (σ, Y ) L jk (Σ * jk , Y ) ≤ log E Σ * jk L jk (σ, Y ) L jk (Σ * jk , Y ) = log R 2 L jk (σ, y)dµ(y) = 0, ( 3
R (Σ * jk ) = E Σ * jk ∂ 2 σ (Σ * jk , Y ) by point 2 of Lemma 3.3 = E Σ * jk ∂ 2 σ L jk (Σ * jk , Y ) L jk (Σ * jk , Y ) - ∂ σ L jk (Σ * jk , Y ) L jk (Σ * jk , Y ) 2 = -E Σ * jk ∂ σ L jk (Σ * jk , Y ) L jk (Σ * jk , Y ) 2 , since E Σ * jk ∂ 2 σ L jk (Σ * jk , Y ) L jk (Σ * jk , Y ) = R 2 ∂ 2 σ L jk (Σ * jk , y)dµ(y) = 0 by point 1 of Lemma 3.3. If we assume that R (Σ * jk ) = 0, then ∂ σ L jk (Σ * jk , Y ) = 0 a.s.
, which contradicts the definition of L jk (Σ * jk , Y ) given in (2.3). Accordingly, there exists > 0 and η > 0 such that R is ( , η) strongly Morse. Assumption "Morse" is satisfied.

For each couple (j, k) such that j < k, Theorem 2 of Mei, S., Bai, Y. and Montanari, A. (2017) applied to the estimator Σ (n) jk yields: Let 0 < ρ < 1. There exists an universal constant C 0 such that letting

n log n ≥ 4C 0 log τ (1 -δ) ρ ∨ 1 τ 2 2 ∨ τ 4 η 2 ∨ τ 2 L 2 η 4 with τ = 2γL 1 and L = sup σ:|σ|≤1-δ R (3) (σ) , P Σ (n) jk -Σ * jk ≤ 2τ η C 0 log n n log τ (1 -δ) ρ ∨ 1 ≥ 1 -ρ.
It follows that, for 0 < ρ < 1 and n such that

n log n ≥ 4C 0 log τ (1 -δ) ρ ∨ 1 τ 2 2 ∨ τ 4 η 2 ∨ τ 2 L 2
η 4 , then:

P Σ (n) -Σ * ∞ ≤ 2τ η C 0 log n n log τ (1 -δ) ρ ∨ 1 ≥ 1 -ρ p(p -1) 2 ,
imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019

where ||A|| ∞ = max j,k∈{1,...,p}

|A jk | denotes the elementwise infinite norm of the matrix A.

This finishes the proof of Proposition 3.1.

Precision matrix estimator and graph recovery

Before giving a result about the estimator of the precision matrix Θ * , let us state a third and last assumption:

(H3) There exists some α ∈]0, 1] such that:

max e∈S c Γ * eS (Γ * SS ) -1 1 = Γ * S c S (Γ * SS ) -1 ∞ ≤ 1 -α,
where:

• if M ∈ M r,m (R), A ⊂ 1, r and B ⊂ 1, m , M AB denotes the matrix (m ij ) i∈A,j∈B , • S = S(Θ * ) := E(Θ * ) ∪ {(1, 1), . . . , (p, p)} where Θ * = (Σ * ) -1 and E(Θ * ) is the set of the edges of the theoretical graph (see (2.2)). Let s := |E(Θ * )|, hence |S(Θ * )| = |E(Θ * )| + p = s + p, • S c = S c (Θ * ) = 1, p 2 \ S(Θ * ), • Γ * = Σ * ⊗ Σ *
where ⊗ denotes the Kronecker matrix product. We have:

Γ * (j,k),(l,m) = cov(X j X k , X l X m ) and thus Γ * SS ∈ M s+p,s+p (R), • ||u|| 1 = d j=1 |u j | for all u ∈ R d is the 1 -norm, • |||U ||| ∞ = max j=1,...,d m k=1 |U jk | for all U ∈ M d,m (R).
The underlying intuition is that this assumption (H3) limits the influence that the non-edge terms, indexed by S c , can have on the edge-based terms, indexed by S (Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B., 2011). Proposition 3.2, set out below, provides an upper-bound for the elementwise maximum norm of the precision matrix estimator Θ (n) obtained by the graphical lasso procedure (2.5). It also shows its sparsistency with respect to graphical model structure recovery. Here are some preliminary notations:

• d is the maximum degree:

d = max j=1,...,p {k ∈ 1, p : Θ * jk = 0} .
(3.9)

• κ Σ * and κ Γ * are defined by: Proposition 3.2. Assume (H3) and assume that there exist some strictly positive constants B, C, D and c > 2 such that letting n log n ≥ C log Bp c , we have:

κ Σ * := |||Σ * ||| ∞ = max
P Σ (n) -Σ * ∞ ≥ D log n n log Bp c ≤ p(p -1) 2p c , (3.12)
where ||A|| ∞ = max j,k∈{1,...,p}

|A jk | denotes the elementwise infinite norm of the matrix A.

Assume that the sample size n is lower bounded as

n log n > D 2 log Bp c max √ C D , 6(1+ 8α -1 )d max{κ Σ * κ Γ * , κ 3 Σ * κ 2 Γ * } 2
, and denote Θ (n) the unique solution of (2.5) and λ n = 8D α log n n log Bp c the regularisation parameter involved in (2.5). Then, with probability greater than 1 -1 p c-2 , we have:

(a) The estimator Θ (n) of Θ * satisfies:

|| Θ (n) -Θ * || ∞ ≤ 2D(1 + 8α -1 )κ Γ * log n n log Bp c . (b)
The estimated edges set is a subset of the true edges set:

E( Θ (n) ) ⊂ E(Θ * ) and E( Θ (n)
) includes all edges (j, k) with:

Θ * jk > 2D(1 + 8α -1 )κ Γ * log n n log Bp c .
Proof relies on results of Theorem 1 of Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B. ( 2011), in which they study the precision matrix estimation problem in the multivariate Gaussian setting.

Proof. (Proposition 3.2) Let us check the two assumptions of Theorem 1 of Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B. (2011).

• Incoherence assumption. This assumption is exactly our assumption (H3).

• Control of sampling noise. A careful reading of Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B. ( 2011) reveals that the tail conditions of their Theorem 1 is not necessary. The required assumption, stated below, is in fact weaker, and is given in Lemma 8 of Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B.

(2011):

There exists v * > 0 such that for all c > 2 and n such that βf (n, p c ) ≤ 1 v * , we have: p 2 , Assumption (3.12) gives this result. Indeed:

P || Σ (n) -Σ * || ∞ ≥ βf (n, p c ) ≤ 1 p c-2 , imsart-
• βf (n, r) = argmax{β : f (n, β) ≤ r} = D log n n log Br • n log n ≥ C log Bp c ⇐⇒ βf (n, p c ) ≤ D √ C
Assumption "Control of sampling noise" is satisfied.

At last, let us set nf (β, r) := argmax{n : f (n, β) ≤ r}. Then, the condition n >

nf (β, r) is equivalent to n log n > log Br D 2 β 2 since f (n, β) ≤ r ⇐⇒ n log n ≤ log Br D 2 β 2 .
We complete the proof by applying Theorem 1 of Ravikumar, P., Wainwright, M. J., Raskutti, G. and Yu, B. ( 2011).

Finally, Propositions 3.1 and 3.2 provide the following theorem, which establishes consistency of the estimator Θ (n) in the elementwise maximum-norm: Theorem 3.1. Assume (H1), (H2) and (H3). Let c > 2, Θ (n) the unique solution of (2.5) and α, d, κ Σ * and κ Γ * respectively defined in (H3), in (3.9), in (3.

|| Θ (n) -Θ * || ∞ ≤ 2D(1 + 8α -1 )κ Γ * log n n log Bp c . (b) E( Θ (n) ) ⊂ E(Θ * ) and E( Θ (n)
) includes all edges (j, k) with:

Θ * jk > 2D(1 + 8α -1 )κ Γ * log n n log Bp c .
In other words, the graph structure of latent Gaussian encoded in Θ * is consistently recovered as long as:

Θ * jk > 2D(1 + 8α -1 )κ Γ * log n n log Bp c .
The parameter c of Theorem 3.1 is a user-defined parameter. The larger c is, the larger the probability for which results of Theorem 3.1 hold is. However, large values of this parameter lead to more stringent requirements on the sample size n.

Simulation studies

Simulation settings

In almost all of this Section (unless otherwise stated), we use the following simulation settings. We simulate n = 500 observations of a p = 100-Gaussian vector X centered and scaled. Graph structure is a chain, that is X 1 ←→ X 2 ←→ • • • ←→ X 100 . The graph contains then 99 edges. Data have been simulated with the R function huge.generator, option graph = "band" of the package huge.

Two different settings of the points of truncation are presented:

• identical points of truncation: a = -0.5 and b = 2, • decreasing points of truncation: a = -1, b = seq(2,0.5, length = p).

We then apply the estimation procedure described in Subsection 2.3. We assume that the points of truncation are known and we use the estimators Σ (n) and Θ (n) of the covariance and precision matrices, respectively defined at Step 1 (2.4) and at Step 2 (2.5).

For the estimation of the precision matrix with the graphical lasso procedure, we use the function huge of the same package, option method = "glasso". Unfortunately, theoretical results of Section 3 do not give an explicit choice of the penalty parameter. We thus choose the penalty parameter with "stars" and "ebic" methods, implemented in the package huge.

Simulations address several problems. Let first explicit the two different procedures used for these simulations, respectively called "our procedure" and "graphical lasso directly on truncated data" (shortened in "Glasso"). The first one is our procedure, which consists in replacing the empirical covariance matrix of X in the graphical lasso by our estimator Σ (n) . The second one is the graphical lasso directly applied to the truncated data, which consists in replacing the empirical covariance matrix of X in the graphical lasso by the empirical covariance matrix of the truncated vector Y . Here are the problems addressed in the following subsections:

• Efficiency of our procedure. Does using our estimator for the covariance matrix really improve graph estimation? • Impact of the points of truncation, that is how the values of the truncation points impact edges detection. • Is our procedure as efficient with other graph structures?

imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019 To study efficiency of these procedures on graph estimation, we make 50 i.i.d. repetitions of the procedure and we represent the detection rates of each of the 100 2 = 4950 potential edges.

Efficiency

In this subsection, we aim at illustrating the efficiency of our procedure. For that, we compare detection rates of each potential edges obtained with our procedure and with graphical lasso directly on truncated data.

Results and comments. Figures 1 and 2 illustrate detection rates for our procedure and for graphical lasso directly on truncated data Y . We represent detection rates of each potential edge by unfolding the matrix with upper.tri, which implies that the first potential edge whose detection rate is displayed is

X 2 ←→ X 1 , then X 3 ←→ X 1 , X 3 ←→ X 2 , X 4 ←→ X 1 , .
. . , X 4 ←→ X 3 and so on. The 99 true theoretical edges, that is edges X i ←→ X i+1 , are displayed with red triangles and in the following order X 1 ←→ X 2 , X 2 ←→ X 3 , ..., X 99 ←→ X 100 . Detection rates of edges of type X i ←→ X i+2 , which are not true edges, are displayed with blue crosses. We distinguish these edges because these interactions can be relatively strong because of the indirect link through X i+1 .

Figure 1 shows identical points of truncation setting. We can observe that our method gives better results: true edges are better detected and other edges are less detected than with Glasso. For example, true edges are all detected more than 96% whereas Glasso detects them between 80% and 100%: 66 of them are detected at most 90%. Edges of type X i ←→ X i+2 tend to be more detected with Glasso.

Figure 2 shows decreasing points of truncation setting. Efficiency of our procedure is even more convincing in this setting. Indeed, true edges are much better detected by our procedure (more thant 80% whereas only 60 true edges are detected at most 80% by Glasso). These differences are probably due to the zero rate in the truncated data Y which grows from 20% to 50% according to the variables. Consequently, edges between variables whose zeros rate is high (that is edges X i ←→ X i+1 for i close to 100) tend to be less detected. The phenomenon is even stronger with Glasso. Besides, the other edges (the false ones) tend to be slightly less detected with Glasso. For our procedure, the false edges are more detected when the truncation points induce a high rate of zero for the involved variables.

Impact of points of truncation

To expand this section, we illustrate the impact of the points of truncation values. For that, we briefly compare results obtained for both "identical" and "decreasing" settings. Observations of the underlying Gaussian vector are the same and we only change values of the points of truncation according to the chosen setting.

Results and comments. Figures 1(a) and 2(a) respectively show detection rates for "identical" and "decreasing" points of truncation settings. The zero inflation of truncated data Y of the first setting is around 33%. Zero inflation of truncated data of the "decreasing" setting decreases from 20% (for Y 1 ) to 50% (for Y 100 ). Detection rates of potential edges are represented by unfolding the matrix with upper.tri, which implies that the first potential edge whose detection rate is displayed is X 2 ←→ X 1 , then X 3 ←→ X 1 , X 3 ←→ X 2 , X 4 ←→ X 1 , . . . , X 4 ←→ X 3 etc. The 99 true theoretical edges, that is the edges X i ←→ X i+1 , are displayed with red triangles and in the following order X 1 ←→ X 2 , X 2 ←→ X 3 , . . . , X 99 ←→ X 100 . Thus, we can observe that edges involving variables whose zero inflation is close to 50% have worse detection rates: true edges are less detected whereas false edges have higher detection rates.

In short and as expected, zero inflation impacts detection rates: the more the zero inflation is, worse is the detection rate. An other phenomenon, not noticeable in Figures 1(a) and 2(a) occurs. This phenomenon is also linked to zero inflation and is noticeable in Figure 3. It points out that detection rates does not only depend on zero inflation but also on the observation window. Indeed, Figure 3 exhibits detection rates for "identical" points of truncation setting, that is a = -0.5 and b = 2, and for an other setting a = -1 and b = 1. For both settings, zero inflation is about 32%. However, results obtained for "identical" points of truncation are much better. The underlying idea is that the observation of the Gaussian variable gives more informations between -0.5 and 2 than between -1 and 1, especially for covariances estimation.

Other graph structures

Previously, we restrict to only one graph structure, the chain structure which tends to give satisfactory results in general. To fulfil these simulation studies, we present some results with other graph structures:

• The "random" structure. There exists an edge between two variables with probability 1/50. Data have been simulated with the R function huge.generator, options graph = 'random', prob = 1/50. Resulting graph has 103 edges. • The "hub" structure. Variables are split into 4 groups of 25. Inside each group, one of the variable is a "hub" and is connected to all the variables of its group. Data have been simulated with the R function huge.generator, options graph = 'hub', g = 4. Resulting graph has 96 edges.

These graphs are displayed in Figure 4. Points of truncation are set to identical (a = -0.5 and b = 2). We simulate n = 500 observations of p = 100 variables and we compare detection rates of each potential edge obtained with our procedure and with graphical lasso directly on truncated data like in Subsection 4.2.

Results and comments. Figure 5 displays these comparisons for "random" and "hub" graph structures.

Results obtained with our procedure are indeed slightly less satisfying than with "chain" structure. Detection rates with "hub" graph structure are a little better than with "random". For "random", all the false edges are detected less than 10% and true edges are all detected more than 46% except for the edge X 34 ←→ X 85 (22%). This edge is also slightly less detected (92% whereas at least 98% for the other true edges) when we directly apply graphical lasso on untruncated Gaussian data X (see Figure 6(a)). Thus, this low detection rate is not only due to our procedure but probably also to the data itself. "Hub" structure gives good results: false edges are detected at most 8% and true edges at most 84%, which is also the case when we apply graphical lasso on untruncated data X (see Figure 6(b)).

In comparison, Glasso (on truncated data Y ) always gives less good results. For "random", the edge X 34 ←→ X 85 is for example never detected and the other true edges are detected between 14% and 82%. "Hub" structure exhibits the most striking difference: true edges are detected at most 54% (but at least 22%).

Discussion

In this paper, we proposed a procedure for graph estimation in a zero-inflated Gaussian model. In this model, zero-inflation is obtained by double truncation (left and right) of Gaussian data. More precisely, the goal is to retrieve the underlying graph structure given by the precision matrix of the Gaussian data with the doubly truncated data. Our procedure includes two steps: the first one consists in estimating the covariance matrix terms to terms by maximising the corresponding bivariate marginal log-likelihood of the truncated vector. The second one relies on the graphical lasso procedure to obtain a lasso estimation of the precision matrix. We then proved some theoretical convergence guarantees with regard to graph estimation. The first result states rate convergence about the covariance matrix estimator. The second one provides sparsistency of the precision matrix estimator with respect to graph structure recovery. Practically, simulations studies also corroborate efficiency of our procedure. They also show that our procedure is more appropriate than using graphical lasso directly on truncated data (without a preliminary estimation of the covariance matrix).

However, results depend on graph structure and simulations exhibit well that some graph structures are more favorable to graph recovery.

This work only deals with a double truncation. Yet, our procedure seems to be practically efficient with a single truncation (right or left) but proof of our theoretical results requires the both right and left points of truncation and they do not hold in the unilateral setting. It would be interessant to address this case later, perhaps using different tools.

Hence, for all m ∈ {1, 2, 3}, ∂ m σ φ ab,jk is continuous on a compact of R 3 , which establishes the result. We naturally distinguish the four cases (see (A.2) for expressions of φ ab,jk ):

• a = b = 1: φ 11,jk is C 3 on ] -1, 1[×R 2 hence on [-1 + δ, 1 -δ] × [a j , b j ] × [a k , b k ]. • a = 0, b = 1: Since F is C ∞ hence C 3 on R, φ 01,jk is C 3 on [-1 + δ, 1 -δ] × [a k , b k ].
• a = 1, b = 0: Analogous to a = 0, b = 1.

• a = b = 0 : Let m ∈ {1, 2, 3}. We apply Lebesgue theorem for continuity and differentiability of integrals with parameters: ∂ l σ φ 11,jk (σ, x, y)dxdy.

Let us deal with each of these terms:

• For a = 0, b = 0: it is obvious.

For the following terms, we use Lebesgue theorem for differentiability of integrals with parameters.

imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019

• φ 11,jk (Σ * jk , x, y) log φ 11,jk (σ, x, y)dxdy.

On the other hand,

E Σ * jk ∂ l σ log L jk (σ, Y ) =

  imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019

  bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019where βf (n, r) := argmax{β : f (n, β) ≤ r} for some function f (n, β).Setting f (n, β)

  10) and in (3.11). There exists some known constants B, C and D such that letting n Bp c the penalisation parameter of the equation (2.5), we have, with probability greater than 1 -1 p c-2 : (a) The estimator Θ (n) of Θ * satisfies:

  Our procedure, identical points of truncation: a = -0.5 and b = 2. Glasso, identical points of truncation: a = -0.5 and b = 2.

Figure 1 .

 1 Figure1. Comparison of detection rates obtained by our procedure and by graphical lasso directly on truncated data. Identical points of truncation setting. Detection rates are obtained on 50 i.i.d. repetitions for n = 500 observations of p = 100 variables. True edges are represented with red triangles, (false) edges of type X i -X i+2 with blue crosses and other false edges with black circles.

  Our procedure, decreasing points of truncation: a = -1, b = seq(2,0.5, length = p). Glasso, decreasing points of truncation: a = -1, b = seq(2,0.5, length = p).

Figure 2 .

 2 Figure2. Comparison of detection rates obtained by our procedure and by graphical lasso directly on truncated data. Decreasing points of truncation setting. Detection rates are obtained on 50 i.i.d. repetitions for n = 500 observations of p = 100 variables. True edges are represented with red triangles, (false) edges of type X i -X i+2 with blue crosses and other false edges with black circles.

  Points of truncation setting: a = -1 and b = 1. Identical points of truncation: a = -0.5 and b = 2.

Figure 3 .

 3 Figure 3. Comparison of detection rates for two points of truncation settings for which the zero inflation is similar (around 32%). Detection rates are obtained on 50 i.i.d. repetitions for n = 500 observations of p = 100 variables. True edges are represented with red triangles.

  (a) "Random" graph structure. (b) "Hub" graph structure.

Figure 4 .

 4 Figure 4. Graphical representation of the two graphs used in this subsection: "random" and "hub".

Figure 5 .

 5 Figure 5. Comparison of detection rates obtained by our procedure (left) and by graphical lasso (right) directly on truncated data. Identical points of truncation setting. "random" and "hub" graph structures. Detection rates are obtained on 50 i.i.d. repetitions for n = 500 observations of p = 100 variables. True edges are represented with red triangles.

Figure 6 .

 6 Figure6. Detection rates obtained with graphical lasso on Gaussian data X. "random" and "hub" graph structures. Detection rates are obtained on 50 i.i.d. repetitions for n = 500 observations of p = 100 variables. True edges are represented with red triangles.

•L

  σ → φ 01,jk (σ, y) is C 3 on [-1 + δ, 1 -δ]. • Straigthforward calculations of derivatives of φ 01,jk (σ, y) w.r.t. σ show that, for σ ∈ [-1 + δ, 1 -δ]: ∂ m σ φ 01,jk (σ, y) ≤ C(a j , b j , a k , b k , δ, m) exp -y 2 2 ,whereC(a j , b j , a k , b k , δ, m) is a positive constant depending on a j , b j , a k , b k , δ and m and y → C(a j , b j , a k , b k , δ, m) exp -y 2 2 is integrable on [a k , b k ] c . It follows that φ 00,jk is C 3 on [-1 + δ, 1 -δ].Appendix C: Proof of Lemma 3.3 Proof. (Lemma 3.3) 1. Let l ∈ N * . First, for all σ ∈ [-1 + δ, 1 -δ], we have: R 2 L jk (σ, y)dµ(y) = jk (σ, y)dµ(y) = ∂ l σ φ 00,jk (σ) + bj ]×[a k ,b k ]

  For a = 0, b = 1 (and a = 1, b = 0): According to (A.2), it is clear that φ01,jk is C ∞ on the compact [1 -δ, 1 + δ] × [a k , b k ],which establishes the formula.• For a = 1, b = 1: Analogously, (A.1) shows that φ 11,jk is C ∞ on the compact [1 -δ, 1 + δ] × [a j , b j ] × [a k , b k ]. 2. Let l ∈ N * .Let us first clarify some notations:E Σ * jk log L jk (σ, Y ) = E Σ * jk (σ, Y ) = R(σ) = R 2 log L jk (σ, y)L jk (Σ * jk , y)dµ(y) = φ 00,jk (Σ * jk ) log φ 00,jk (σ) + b k a k φ 01,jk (Σ * jk , y) log φ 01,jk (σ, y)dy + bj aj φ 10,jk (Σ * jk , x) log φ 10,jk (σ, x)dx + [aj ,bj ]×[a k ,b k ]

  log L jk (σ, y) L jk (Σ * jk , y)dµ(y) = φ 00,jk (Σ * jk )∂ l σ log φ 00,jk (σjk (Σ * jk , y)∂ l σ log φ 01,jk (σ, y)dy + bj aj φ 10,jk (Σ * jk , x)∂ l σ log φ 10,jk (σ, x)dx+ [aj ,bj ]×[a k ,b k ]φ 11,jk (Σ * jk , x, y)∂ l σ log φ 11,jk (σ, x, y)dxdy.To show the equality∂ l σ E Σ * jk log L jk (σ, Y ) = E Σ * jk ∂ l σ log L jk (σ, Y) , we show the equality for each of the four terms• For a = 0, b = 0: it is obvious. For the three remaining terms, we use Lebesgue theorem for differentiability of integrals with parameters.• For a = 0, b = 1 (and a = 1, b = 0): log φ 01,jk (σ, y) = -The function log φ 01,jk is C ∞ on the compact [-1 + δ, 1 -δ] × [a k , b k ]. Therefore, for all σ ∈ [-1 + δ, 1 -δ] and y ∈ [a k , b k ], ∂ l σ log φ 01,jk (σ,y) is upper bounded by imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019 a constant which is integrable on the compact [a k , b k ] (with regard to the density y → φ 01,jk (Σ * jk, y)). • For a = 1, b = 1: log φ 11,jk (σ, x, y) = -log(2π) -1 2 log(1 -σ 2 ) -x 2 -2σxy + y 2 2(1 -σ 2 ) .Analogously, the function logφ 11,jk is C ∞ on the compact [-1 + δ, 1 -δ] × [a j , b j ] × [a k , b k ]. Thus, for all σ ∈ [-1+δ, 1-δ], x ∈ [a j , b j ] and y ∈ [a k , b k ], ∂ l σ log φ 11,jk (σ, x,y) is upper bounded by a constant which is integrable on the compact [a j , b j ] × [a k , b k ] (with regard to the density (x, y) → φ 11,jk (Σ * jk , x, y)).

imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019

Acknowledgements

We would like to thank Stéphane Robin and Stéphane Chrétien for their expert advice and comments.

Appendix A: Proof of Lemma 3.1

Proof. (Lemma 3.1) It is sufficient to show the existence of such a constant γ jk > 0 for

The proof naturally falls into four parts:

• a = b = 1:

is continuous and strictly positive on a compact of R 3 .

• a = 0, b = 1: An easy computation yields that

where F denotes the c.d.f. of N (0, 1). Since

,jk is hence continuous and strictly positive on a compact of R 2 .

• a = 1, b = 0: Analogous to a = 0, b = 1.

• a = b = 0:

Since δ 2 ≤ 1 -σ 2 ≤ 1, φ 00,jk is hence continuous and strictly positive on a compact of R.

Appendix B: Proof of Lemma 3.2

Proof. (Lemma 3.2) In the same manner as for the proof of Lemma 3.1, it is sufficient to show the result for j < k fixed. Fix j < k and show that for all a, b ∈ {0, 1}, the function φ ab,jk is

imsart-bj ver. 2014/10/16 file: "TruncGaus -arXiv".tex date: November 17, 2019


Rodolfo Dufo-López

Daniel Champier

Stéphane Gibout

Juan Lujano-Rojas

José Domínguez-Navarro

Juan Lujano

Thermoelectric generation as a part of the optimal solution in off-grid hybrid renewable system in places with low irradiation

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Thermoelectric generation as a part of the optimal solution in off-grid hybrid renewable system in places with low irradiation.

17th European Conference on Thermoelectrics

Sept 2019 / Limassol, Cyprus -Optimise the system to supply electricity in off-grid locations , minimizing the total cost during the system lifetime -Demonstrate the economic viability of including TEG (in a cook or heat stove) in the optimal system Conclusion TEG has been included into the simulation and optimisation software iHOGA, developed by researchers of the University of Zaragoza (Spain). TEG standard simplified model combined with the TEM manufacturer output power curves lead to good simulation results, compared to the experimental results obtainted by the SIAME of UPPA. TEG can be part of the optimal system to supply the electrical load in places with low irradiation. The optimal solution for a household in south Norway is a PV-TEG-FFG-BESS system. TEG contributes to enlarge the battery lifetime and reduce the fuel consumption of the FFG Contact : daniel.champier@univ-pau.fr

Results

 The model has been applied in the optimization of the electricity supply of a off-grid household in the south of Norway with heating pellet stove.  The optimal system is a PV-TEG-FFG-BESS system  TEG contributes to enlarge the battery lifetime and reduce the fuel consumption of the FFG.

Validation

Fig

  Fig.1 Possible components of the hybrid system Fig. 3. iHOGA software (https://ihoga.unizar.es/en/)

  López 1 , Daniel Champier 2 , Stéphane Gibout 3 , Juan M. Lujano-Rojas 1 , José A. Domínguez-Navarro 1 1. University of Zaragoza, Electrical Engineering Department, C/María de Luna, 3, Edificio Torres Quevedo, 50018 Zaragoza, Spain. 2. Univ Pau & Pays Adour/E2S UPPA, Laboratoire des Sciences de l'Ingénieur Appliquée à la Mécanique et au Génie Electrique, Fédération IPRA, EA4581, Technopôle Hélioparc Pau Pyrénées 2, Av du Président Angot, 64000 Pau, France. 3. Univ Pau & Pays Adour/E2S UPPA, Laboratoire de Thermique, Energétique et Procédés IPRA EA1932, Rue Jules Ferry, BP 7511, 64 000 Pau, France.

Objectives

Project supported by the Spanish Gobernment "Ministerio de Educación, Cultura y Deporte. Gobierno de España", Conseil Régional Nouvelle Aquitaine and Communauté d'Agglomération Pau-Pyrénées . Dufo-López R et Al "Optimisation of off-grid hybrid renewable systems with thermoelectric generator." Energy Conversion and Management, 2019 Rodolfo Dufo-