In the context of finite mixture models one considers the problem of classifying as many observations as possible in the classes of interest while controlling the classification error rate in these same classes. Similar to what is done in the framework of statistical test theory, different type I and type II-like classification error rates can be defined, along with their associated optimal rules, where optimality is defined as minimizing type II error rate while controlling type I error rate at some nominal level. It is first shown that finding an optimal classification rule boils down to searching an optimal region in the observation space where to apply the classical Maximum A Posteriori (MAP) rule. Depending on the misclassification rate to be controlled, the shape of the optimal region is provided, along with a heuristic to compute the optimal classification rule in practice. In particular, a multiclass FDR-like optimal rule is defined and compared to the thresholded MAP rules that is used in most applications. It is shown on both simulated and real datasets that the FDR-like optimal rule may be significantly less conservative than the thresholded MAP rule.

Introduction

Consider a sample X 1 , . . . , X n of independent observations in space X stemming from P populations. The distribution of such a sample can be modelled as a finite mixture of distributions

f (x) = P ∑ p=1 π p f p (x)
where each class is described by its own probability density distribution f p and a weight π p so that 0 ≤ π p ≤ 1 and ∑ p π p = 1. For each observation a label variable Z can be introduced, that equals p if the observation belongs to class p. In the following, the labels Z 1 , . . . , Z n are supposed to be independent, and the mixture model is assumed to be completely known, i.e. the number of classes P, the class distributions f 1 , . . . , f P and the weights π 1 , . . . , π P are known.

In most practical situations the labels are unobserved, and the goal is to find a suitable clustering of the sample. This is typically obtained by applying a classification rule, i.e. a function ψ that maps X into {1, . . . , P}. The key quantity for building such a function is the posterior probability for observation x to belong to class p:

τ p (x) = Pr (Z = p|X = x) = π p f p (x) ∑ P p =1 π p f p (x)
.

(1)

The Maximum A Posteriori (MAP) classification rule, defined as

ψ MAP (x) = argmax 1≤p≤P τ p (x)
is by far the most popular classification rule. In particular, it is known to minimize the classification error rate P(ψ(X) = Z), [START_REF] Mclachlan | Finite mixture models[END_REF].

While optimality in terms of classification error rate is a desirable property, two drawbacks of the MAP rule should be mentioned. First, optimality does not prevent against a high level of misclassification. When the classification task is difficult, as much as half of the observations can be misclassified by the MAP rule. Second, the MAP rule does not account for the asymmetry that may exist between classes: in some real applications only a small number among the P classes may be of major importance for the experimenter. One should then focus on the misclassification rate on these classes of interest. This situation arises in most cases involving comparisons between different conditions, where one is interested in identifying units having different behaviors in different conditions, but not units whose behavior is unchanged across conditions. For instance in [START_REF] Bérard | Unsupervised classification for tiling arrays: Chip-chip and transcriptome[END_REF], a methylation analysis was conducted to identify regions of the genome that are differentially methylated between different organs of Arabidopsis thaliana. This problem can be cast into a clustering problem where units are probes (locations) spread over the Arabidopsis genome, and classes are defined according to the methylation status of these probes in the different organs. Only classes corresponding to differentially methylated behaviors are then of interest.

Different solutions have been proposed to circumvent these two problems. One possible approach is to take into account the asymmetry between classes by selecting unequal misclassification cost functions that emphasize the cost of misclassification c k from a noninteresting class into a class k of interest. This has been investigated in [START_REF] Friedman | The elements of statistical learning: data mining, inference, and prediction[END_REF] for instance.

On the other hand, a high level of misclassification for the MAP rule may come from the fact that it classifies any observation x, regardless of the uncertainty of the classification of point x. This uncertainty may be quantified through the probability to not belong to the class x was classified into. A possible solution is to perform a partial classification, using e.g. the thresholded classification rule where 0 < α < 1 is a parameter to be chosen by the experimenter. The condition τ p (x) > 1α can be reformulated has 1τ p (x) < α which amounts to focusing on the classification of observations for which uncertainty is not too high.

ψ α (x) = p if p = argmax
The design of partial rules has been considered through the use of classification rules with a reject (or abstention) option. An often considered setting there is that the reject option has a fixed cost (less costly than a misclassification error); see [START_REF] Chow | On optimum recognition error and reject tradeoff[END_REF], Herbei and [START_REF] Herbei | Classification with reject option[END_REF], [START_REF] Pillai | Multi-label classification with a reject option[END_REF]. This setting has been also considered in conjunction with specific machine learning methods such as large margin classifiers [START_REF] Bartlett | Classification with a reject option using a hinge loss[END_REF][START_REF] Grandvalet | Support vector machines with a reject option[END_REF][START_REF] Wegkamp | Support vector machines with a reject option[END_REF], active learning [START_REF] Zhang | Beyond disagreement-based agnostic active learning[END_REF], and fair learning [START_REF] Schreuder | Classification with abstention but without disparities[END_REF].

Alternative partial classification rules include the tight clustering algorithm introduced in [START_REF] Tseng | Tight clustering: A resampling-based approach for identifying stable and tight patterns in data[END_REF] and extended to deal with large datasets in [START_REF] Karmakar | Tight clustering for large datasets with an application to gene expression data[END_REF]. In the latter work, when applying the extended tight clustering algorithm to a dataset consisting of more than 50,000 gene expression probes for individuals suffering from psoriasis, more than 30,000 probes were not classified in any of the six clusters that were identified.

All aforementioned approaches (partial rules, rules with different misclassification costs and rule with abstention costs) are attractive, however there exist very few guidelines about how to tune parameters such as α, c k , or the cost of abstaining, and what will be the impact of this tuning on the final misclassification rate. In this article, we propose a method that addresses both the problem of partial classification and the presence of classes of no interest. In particular, we propose to design the abstention region based on a directly interpretable constraint taking the form of an overall statistical confidence level on the effective classifications. Based on the control of type I error rate in statistical test theory, we introduce optimal classification rules that guarantee that as many observations as possible are classified, under the constraint that the rate of misclassification in the classes of interest is controlled. As such, our method can be understood as an extension of the Bayesian False Discovery Rate (BFDR) developed in [START_REF] Efron | Empirical bayes methods and false discovery rates for microarrays[END_REF] to general finite mixture models. It is shown that different classification rules should be considered to achieve optimality according to the number K of classes of interest (K = P or K < P). Sub-optimality of the thresholded classification rule ( 2) is also demonstrated.

The principle of controlling a specific error rate related to classical statistical confidence criteria has also been considered in, e.g. [START_REF] Scott | A Neyman-Pearson approach to statistical learning[END_REF] and [START_REF] Scott | The false discovery rate for statistical pattern recognition[END_REF] in the context of binary classification; see also [START_REF] Tong | A survey on neyman-pearson classification and suggestions for future research[END_REF] for a more recent survey. In that setting, an asymmetry is introduced between the classes: the goal is to maximize correct classification rate in class 1 subject to a fixed control at a prescribed level of the classification error in class 0, either in absolute-value (Neyman-Pearson classification) or in the sense of FDR. The setting we introduce here is multi-class and includes the possibility of abstention. Furthermore all the classes of interest have the same status (the classes of no interest have a different status, however even if all classes are of interest, our setting makes sense and abstention remains a possibility). The goal we aim at is to guarantee a certain confidence level over classified examples while minimizing the probability of abstention of classes of interest. The statistical criteria we introduce generalize the confidence/coverage tradeoff considered in El-Yaniv and Wiener (2010), [START_REF] Wiener | Agnostic pointwise-competitive selective classification[END_REF], [START_REF] Denis | Consistency of plug-in confidence sets for classification in semi-supervised learning[END_REF], who only considered the binary classification case where all classes are of interest.

The problem of classification with confidence has also been considered under a related but different angle by [START_REF] Lei | Classification with confidence[END_REF]. In Lei's approach, only the binary classification setting is considered and two regions C 0 and C 1 , where the respective classes are predicted, are constructed such that C i has a prescribed coverage 1α i for class i and the "ambiguous classification region" C 0 ∩ C 1 has minimal overall probability. We observe that the latter overlap region could be interpreted as similar to an abstention region in our setting. However, the criterion used by [START_REF] Lei | Classification with confidence[END_REF] is different from ours, and finding the optimal classification regions amounts to determine appropriate thresholds for two different level sets of the posterior class probability function. Furthermore, it is unclear how to extend the considered criterion to more than two classes. In our approach, we consider a global constraint which in our view is more readily interpretable in terms of classification confidence, is naturally adapted to more than two classes, and only necessitates the determination of a single level set of an appropriate function.

The paper is organized as follows. In Section 2, notions corresponding to type I and type II error rates are defined in the multiclass classification framework along with optimality for a classification rule. Optimal rules are then exhibited in Section 3, and heuristic procedures to estimate these rules are presented. The empirical behavior of the procedures is briefly investigated and an application of the proposed strategy to the analysis of differential methylation is presented in Section 4. Lastly, some discussion is developed in Section 5.

Definitions

Restricted classification rules

Because the classification task may be arbitrarily difficult for some points of X , we consider restricted rules ψ R (x), that classify points in a subset R ⊆ X only:

ψ R : X → {0, 1, . . . , P},
where 0 is the status "not classified" given to observations in R, the complementary set of R in X : ∀x ∈ R, ψ R (x) = 0.

As an illustration, consider the thresholded classification rule ψ α provided in equation ( 2), and define

τ * (x) = max 1≤p≤P τ p (x),
which corresponds to the maximal posterior probability for point x. Then classification rule ψ α corresponds to the restricted classification rule ψ R α where

R α = {x|τ * (x) ≥ 1 -α}.
In some cases, the number K of classes of interest may be lower than the total number of classes P in the mixture. In the following, and without loss of generality, we will assume that the classes of interest are the first K classes of the mixture. In this context where K < P, it will be convenient to consider classification rules that are also restricted to predict labels in 1, . . . , K:

ψ K,R : X → {0, 1, . . . , K}.
In this situation, the maximal posterior probability of interest is

τ * K (x) = max 1≤k≤K τ k (x).
Note that in the case K = P one has τ * P (x) = τ * (x). One can extend the definition of the thresholded classification rule to the case K < P by defining

ψ α K (x) = k if k = argmax 1≤k ≤K τ k (x) and τ * K (x) > 1 -α, 0 (not classified) otherwise, . (3)

Error rates

A "good" classification rule should classify as many observations as possible while controlling the classification error rate. As in the statistical testing setting, one can introduce the type I and type II error rates associated to any given classification rule. We consider two definitions of the type I error rate and one for the type II error rate in the multiclass classification context:

Definition 2.1 Let {1, . . . , K} be the classes of interest (possibly with K = P) and ψ K,R be a restricted classification rule. Quantities

MNPR(ψ K,R ) = P (ψ K,R (X) = Z, ψ K,R (X) > 0) = P (ψ K,R (X) = Z, X ∈ R) and MFDR(ψ K,R ) =      P (ψ K,R (X) = Z| ψ K,R (X) > 0) = P (ψ K,R (X) = Z| X ∈ R) , if P (X ∈ R) > 0 0, if P (X ∈ R) = 0,
are called the multiclass Neyman-Pearson error rate and the multiclass false discovery rate, respectively. The quantity

MFNR(ψ K,R ) = P (Z ∈ {1, . . . , K}, ψ K,R (X) = 0) = P (Z ∈ {1, . . . , K}, X ∈ R)
is called the multiclass false negative rate.

The MFNR quantifies the proportion of observations that belong to a class of interest and that were not classified by rule ψ K,R . This quantity corresponds to the false negative proportion in the testing setting, and should be as small as possible. The MNPR and MFDR correspond to possible classification error rates one may want to control at a given level α. Note that a criterion equivalent to the MFDR has been considered by [START_REF] El-Yaniv | On the foundations of noise-free selective classification[END_REF], [START_REF] Wiener | Agnostic pointwise-competitive selective classification[END_REF], [START_REF] Denis | Consistency of plug-in confidence sets for classification in semi-supervised learning[END_REF], albeit only in the situation K = P = 2.

To exemplify the previous definitions, here again one can consider the thresholded classification rule ψ α (= ψ α P ) defined above in the by-default context where all classes are of interest. One can show that the MFDR of ψ α is always lower than α, whatever the number of classes of interest. Indeed one has

MFDR(ψ α ) = P (ψ α (X) = Z| X ∈ R α ) = P (ψ α (X) = Z ∩ X ∈ R α ) P (X ∈ R α )
.

The numerator in the last expression can be reformulated as follows:

P (ψ α (X) = Z ∩ X ∈ R α ) = X E Z|X=x I {ψ α (x) =Z} I {x∈R α } f (x)dx = X E Z|X=x [Z = p * ] I {x∈R α } f (x)dx = X (1 -τ * P (x))I {x∈R α } f (x)dx
where p * = arg max p τ p (x). Plugging this expression into the previous equation leads to

MFDR(ψ α ) = 1 P (X ∈ R α ) R α (1 -τ * P (x)) f (x)dx ≤ 1 P (X ∈ R α ) R α α f (x)dx ≤ α.

Optimal classification rules

We now introduce a formal definition for the optimality of a restricted classification rule.

Definition 2.2 Given a level α (0 ≤ α ≤ 1) and a set of classes of interest {1, . . . , K} (possibly with K = P), a classification rule

ψ * K,R is MNPR-optimal at level α if ψ * K,R = arg min ψ K,R MFNR(ψ K,R ) u.c. MNPR(ψ K,R ) ≤ α. (4) Alternatively, ψ * K,R is MFDR-optimal at level α if ψ * K,R = arg min ψ K,R MFNR(ψ K,R ) u.c. MFDR(ψ K,R ) ≤ α. (5)
Note that the definition implies both a region R as large as possible and an optimal classification at each point x in region R.

Optimal classification rules

The goal of the present section is to exhibit the shape of the optimal classification rules for problems (4) and (5). To this end, we will first prove that for a fixed region R it is optimal to apply the restricted MAP rule. As a consequence, looking for an optimal rule ψ * R actually boils down to finding the optimal region R * where to apply the (restricted) MAP rule. The theoretical form of the optimal region will then be derived for problems (4) and ( 5). Each optimal region requires the tuning of a unknown threshold λ , for which we will provide a heuristic estimation method.

Optimal rule when R is fixed

Let us first define the restricted MAP classification rule ψ MAP K,R as

ψ MAP K,R (x) = arg max 1≤k≤K τ k (x) if x ∈ R 0 otherwise.
Let first consider a classification ψ K,R restricted to the region R. At any point x ∈ R one has

P (ψ K,R (x) = Z) = K ∑ k=1 (1 -τ k (x))I {ψ K,R (x)=k} ≥ (1 -τ * K (x)) ≥ P ψ MAP K,R (x) = Z
and therefore

MNPR(ψ K,R ) ≥ MNPR(ψ MAP K,R ) and MFDR(ψ K,R ) ≥ MFDR(ψ MAP K,R ).
Alternatively it is straightforward to observe that

MFNR(ψ K,R ) = MFNR(ψ MAP K,R ).
As a consequence, one can conclude that for a fixed region R the optimal restricted classification rule is ψ MAP K,R . Therefore finding the solutions of problems ( 4) and ( 5) boils down to finding the region R * to which the MAP classification rule should be applied.

In the following the upperscript MAP and the lowerscript K in ψ MAP K,R will be skipped for simplicity, and ψ R should be understood as ψ MAP K,R , with K = P or K < P depending on the context.

Optimal classification region for MNPR control

If all classes are of interest (K = P), one looks for a region R * such that

R * ∈ arg min R P (X ∈ R) u.c. P (ψ R (X) = Z, X ∈ R) ≤ α. If K < P, the region R * should satisfy R * ∈ arg min R P (Z ∈ {1, . . . , K}, X ∈ R) u.c. P (ψ R (X) = Z, X ∈ R) ≤ α. (6) Proposition 1 If MNPR(ψ X ) ≤ α, then R * = X is a solution of (6). Otherwise, any R * satisfying MNPR(ψ R * ) = α and such that x ∑ K k=1 τ k (x) 1 -τ * K (x) > λ ⊆ R * ⊆ x ∑ K k=1 τ k (x) 1 -τ * K (x) ≥ λ , (7)
for some λ ≥ 0, is a solution of (6).

If the marginal distribution of X is without atoms, there exists R * satisfying the above conditions.

In the particular case where K = P, condition (7) can be reduced to:

{ x | τ * (x) > λ } ⊆ R * ⊆ { x | τ * (x) ≥ λ } . ( 8
)
for some λ ≥ 0.

The proof is based on the original proof of the Neyman-Pearson theorem, [START_REF] Neyman | On the problem of the most efficient tests of statistical hypotheses[END_REF], and can be found in Appendix A.

Optimal classification region for MFDR control

If all classes are of interest, one looks for a region R * such that

R * ∈ arg min R P (X ∈ R) u.c. P (ψ R (X) = Z| X ∈ R) ≤ α. If K < P, R * should satisfy R * ∈ arg min R P (Z ∈ {1, . . . , K}, X ∈ R) u.c. P (ψ R (X) = Z| X ∈ R) ≤ α. (9) Proposition 2 If MFDR(ψ X ) ≤ α, then R * = X is a solution of (9). If τ * K (X) < 1 -α almost surely, then R * = / 0 is a solution of (9). Otherwise, any R * satisfying MFDR(ψ R * ) = α and such that x 1 -α -τ * K (x) ∑ K k=1 τ k (x) < λ ⊆ R * ⊆ x 1 -α -τ * K (x) ∑ K k=1 τ k (x) ≤ λ , (10)
for some λ ≥ 0, is a solution of (9); if the marginal distribution of X is without atoms, such an R * exists.

In the case K = P, condition (10) can be reduced to:

{ x | τ * (x) > 1 -α -λ } ⊆ R * ⊆ { x | τ * (x) ≥ 1 -α -λ } , (11)
for some λ ≥ 0.

The proof of Proposition 2 can be found in Appendix B. In the case K = P, the optimal region takes the form of a thresholded rule (2) for an appropriate threshold. This was observed by [START_REF] Denis | Consistency of plug-in confidence sets for classification in semi-supervised learning[END_REF] who studied optimal rules for the MFDR/coverage tradeoff in the case K = P = 2. Interestingly, when K < P, the optimal rule is different from a simple threshold rule.

The form of the optimal region R * in ( 8), ( 10) is C <λ ⊆ R * ⊆ C ≤λ , where C <λ ,C ≤λ are sublevel sets (in the sense of strict, resp. non-strict inequality) of a criterion depending on the posterior probabilities at point x. This "sandwiching" relation is theoretically relevant in situations where P (C <λ ) < P C ≤λ , for example if the criterion used is piece-wise constant, and the level set at value λ has nonzero probability. In this situation, strictly speaking to achieve exactly the target MNPR or MFDR rate, R * should only include part of the level set at value λ (and this part can be chosen arbitrarily provided the constraint is satisfied). The assumption of atom-free X-marginal could also be lifted to grant the existence of optimal rules in general, provided randomized rules are allowed, as in classical Neyman-Peason theory. Since such situations are obviously of little relevance for most applications, in the remainder of this work, to simplify exposition we assume that P (C <λ ) = P C ≤λ and that R * exactly coincides with a sublevel set.

Comparison between thresholded rule and the MFDR-optimal rule

The restriction region R α of the thresholded rule (3) can now be compared to the MFDRoptimal regions R * found in the previous section. One has

R α = {x |τ * K (x) ≥ 1 -α } with K ≤ P, R * = x|τ * K (x) ≥ 1 -α -λ K ∑ k=1 τ k (x) with λ ≥ 0, if K ≤ P.
According to the shape of R * , the threshold for classifying an observation x based on the maximal posterior probability τ * K (x) depends not only on α, but also on ∑ K k=1 τ k (x), that is the probability of belonging to the overall group of populations of interest. When K < P this can be interpreted as follows: when τ * K (x) is mild, the optimal rule may still classify x if the probability of x to belong to the overall group of the populations of interest is high. This interesting feature is not accounted for in the thresholded rule that classifies an observation only based on τ * K (x). In this sense, the thresholded rule is not based on the appropriate classification criterion when one is interested only in a subset of all possible populations. More generally, one can observe that whatever the number of populations of interest, both ψ R * and ψ α classify observations based on τ * (x), but using a threshold that is always lower for ψ R * . Consequently i) the MFNR will always be higher for ψ α , and ii) the MFDR of ψ α will be lower than the requested nominal level, i.e. the thresholded rule is a conservative classification strategy to control the MFDR at level α. In Section 4 we illustrate on simulated and real data the fact that the gap in term of MFNR between ψ R * and ψ α may be high.

Estimation of parameter λ

We present here an heuristic strategy to choose λ from the data at hand. First note that all optimal regions defined in Sections 3.2 and 3.3 are of the form of level sets

R * = { x | Crit(x) > λ } ,
where Crit(x) is a criterion based on the posterior probabilities at point x, that depends on the risk one wants to control and on the number of classes of interest as follows:

if one aims at controlling the MNPR then

Crit(x) = ∑ K k=1 τ k (x) 1 -τ * K (x) (case K < P) or Crit(x) = τ * (x) (case K = P),
if one aims at controlling the MFDR then

Crit(x) = τ * K (x) + α -1 ∑ K k=1 τ k (x) (case K < P) or Crit(x) = τ * (x) + α -1 (case K = P)
where λ ≤ 0.

Assuming an n-sample D n = (X 1 , . . . , X n ) is available, the MNPR and MFDR can be estimated for any region R by

P (ψ(X) = Z, X ∈ R) = 1 n ∑ X i ∈R (1 -τ * K (X i )) and P (ψ(X) = Z| X ∈ R) = 1 n R ∑ X i ∈R (1 -τ * K (X i ))
respectively, where n R is the number of observations in D n belonging to R (and we recall that τ * K (x) = τ * (x) if K = P). This leads to the following general heuristic for the evaluation of threshold λ : Algorithm 1 Require: α, x 1 , . . . , x n 1/ Order observations x 1 , . . . , x n according to Crit(x i ):

Crit(x (1) ) ≥ . . . ≥ Crit(x (N) ), 2/ Find the largest index i max such that 1 M(i max ) i max ∑ i=1 (1 -τ * K (x (i) )) ≤ α, 3/ Set λ = Crit(x (i max ) ), 4/ Output λ .
If the goal is to control MFDR, then in the second step of the heuristic M(i max ) = i max , otherwise M(i max ) = n for a MNPR control.

Applications

The aim of this section is to illustrate the performance of the MFDR classification rule derived in the previous section. We first consider different scenarios based on simulated data. In a first scenario the true posterior probabilities are available, which corresponds to the theoretical setting of Sections 2 and 3. In this scenario the only parameter to be estimated is λ , and we empirically evaluate the ability of the previous algorithm to efficiently estimate λ and control the error rate at the required nominal level. The performance of the optimal rule is also compared to the performance of the thresholded rule. In a second scenario the true posterior probabilities are assumed to be unknown -the by-default setting of most application cases. This scenario allows us to evaluate the impact of estimating the posterior probabilities on the error rate control procedure. Lastly, we present an application on real data, where the MFDR control procedure is applied to the differential analysis of methylation profiles.

Simulation setting

Datasets are simulated from a mixture of 3 bidimensional Gaussian distributions whose mean vectors are (-1, 0), (0, D) and (1, 0), respectively. Here D is a parameter that tunes the distance between populations 1 and 3 and population 2. When D = 0 the three populations are highly overlapping, whereas a high value of D makes the populations more distinct. In each class the covariance matrix is diagonal, with an identical variance σ 2 . Weights π 1 , π 2 , π 3 are all fixed at 1/3. For each dataset, 200 observations per population are simulated. In this context the easiness of the classification task is ruled by parameters D and σ 2 : the higher D (respectively the lower σ 2 ), the easier the classification. Several values are considered for D (0, 1, 2, 3) and σ 2 (0.5, 1, 2). An illustration of different configurations ranging from an "easy" classification case (D = 3, σ 2 = 0.5) to a "hard" one (D = 0, σ 2 = 2) can be found in Appendix C. Finally, 100 datasets are generated for each configuration.

In a supplementary simulation study, datasets were generated from a mixture of three bivariate Student distributions; details can be found in Appendix F.

Case 1: posterior probabilities are known

In what follows the objective is to perform classification with a MFDR controlled at a nominal threshold fixed at α = 0.05.

In Section 2.2, a straightforward calculation has shown that the thresholded classification rule ψ α with a threshold fixed at 1α guarantees a MFDR control at a level α. However the results of Section 3.3 suggest that fixing the threshold at 1α may result in a high MFNR. The alternative procedure for the optimal rule described in Section 3.5 should yield better performance thanks to a less conservative choice for the threshold. The potential gain (in terms of MFNR) of the optimal rule compared to the 1α thresholded rule is investigated in this section.

The analysis is performed as follows: for each dataset the posterior probabilities are computed for each observation using the true parameters of the model. The procedure (10) is then applied to estimate the threshold λ in order to control the MFDR at nominal level α = 0.05. Observations are classified using either the MAP rule, the thresholded rule (with α = 0.05) or the estimated optimal rule (i.e. the optimal rule applied with estimated threshold λ ). For all rules the actual MFDR and MFNR are evaluated by comparing the predicted and true labels of the classified observations. Results are displayed in Figure 1 for the case where all classes are of interest. The case where only classes 1 and 3 are of interest is provided in Appendix D.

In the case where σ 2 = 0.5 one can observe that all classification rules achieve a low MFDR -that may still be much higher than 0.05 in the case of the MAP rule. Both the thresholded and the optimal rule efficiently achieve control of the MFDR at the nominal level, but with significant differences in terms of MFNR. As illustrated in Figure 2, the estimated threshold for the optimal rule may be much lower (below 0.5 in some configurations) than the 1α = 0.95 value used in the thresholded rule. As a consequence, the MFNR is twice higher for the thresholded rule than for the optimal rule, illustrating how conservative the thresholded rule can be.

In the case where σ 2 = 2 the classification task becomes quite difficult, and the MAP rule yields MFDR that are higher than 50%. In such configurations both the thresholded q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 1 2 q q q q q q q q q q q q q q q q q q q q q q q q q 0.4 and the optimal rules do not classify any observation in most cases, since the highest posterior probability observed (in a class of interest) is lower than 0.95.

MFDR

In conclusion, whatever the configuration, the optimal rule ψ R * with estimated threshold λ controls the misclassification rate as requested, whereas the thresholded classification rule is more conservative. When the classification task is easy or intermediate the MFNR of the optimal rule can be much smaller than the MFNR of the thresholded classification rule. Results are unchanged if the prior proportions of the mixture vary (not shown), and/or if only a subset of the three classes are of interest (see Appendix D).

Case 2: posterior probabilities are inferred

So far we assumed that the true parameters of the mixture are known, but in practice these parameters are usually estimated from the data at hand. To evaluate the impact of parameter estimation on misclassification rate, we performed the same analysis as in the previous section, except that parameters are now supposed to be unknown and are estimated using the mclust R package [START_REF] Scrucca | mclust 5: clustering, classification and density estimation using Gaussian finite mixture models[END_REF], fixing the number of classes to its true value (i.e. 3). Results are displayed in Figure 3.

When the classification problem is tractable (low or moderate values of σ 2 and/or high values of D) both the thresholded and the optimal rules efficiently control the MFDR. As in the previous scenario, a significant gain in terms of MFDR is observed when using the optimal rule rather the thresholded rule. When the classification problem becomes too difficult (i.e. classes strongly overlap) the estimated posterior probabilities get inaccurate, and any classification rule based on these quantities becomes irrelevant. Although neither the thresholded nor the optimal rule control the MFDR at the nominal level, it is still worth to use one of these classification rules rather than the MAP rule. Considering the case where only a subset of classes is of interest leads to similar conclusions (see Figure 7 in Appendix D).

The performance of the optimal rule depends on the accuracy of the posterior probability estimates. In order to explore the robustness of our method under model misspecification that could possibly lead to biased posterior probabilities estimates, we generated datasets according to a mixture of three Student bidimensional distributions with the same centers and covariance matrices as in the main Gaussian simulation scheme of "intermediate" difficulty (σ = 1). We then estimated the posterior probabilities using the EM algorithm for Gaussian mixtures implemented in mclust. Results are shown in Appendix F. As expected, for higher degrees of freedom, the error rates closely match the ones in the central column (σ = 1) of Figure 3. For lower degrees of freedom, performance breaks down in the "difficult" configurations given by low values of D. For more tractable configurations, MFDR control is generally achieved and we do not observe notably differences q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 1 2 in terms of MFNR with respect to the Gaussian simulations.

MFDR

Application to transcriptomic data

We consider the unsupervised classification problem described in [START_REF] Bérard | Unsupervised classification for tiling arrays: Chip-chip and transcriptome[END_REF]. A methylation experiment was performed to compare the methylation profiles of two organs (leaf and seed) of Arabidopsis thaliana. The methylation profiles were measured on a same plant using a tiling array technology. The genomic sequence of Arabidopsis thaliana is represented on the array by approximately 5 × 10 5 probes, covering both genic and intergenic regions. For each probe, the methylation signal is measured in the two organs.

From a statistical point of view the sampled population is the population of probes, each of them being described by a bivariate signal (hybridization in leaf and seed), and the goal of the analysis is to identify differentially methylated probes, i.e. probes whose methylation signals in leaf and seed differ. In the following only probes corresponding to chromosome 4 (107,199 probes) are considered. A more thorough description of the data can be found in [START_REF] Bérard | Unsupervised classification for tiling arrays: Chip-chip and transcriptome[END_REF].

In [START_REF] Bérard | Unsupervised classification for tiling arrays: Chip-chip and transcriptome[END_REF] a constrained Gaussian bivariate mixture model was fitted to the data. Comparing two samples requires distinguishing four different classes of probes that can be biologically interpreted as follows : a class of probes with low methylation signals ( 1), a class of probes exhibiting a similar intensity level in both organs (2), a class of probes with lower methylation intensities in leaf compared to seed (3), and a symmetric class of probes with higher intensities in leaf compared to seed (4) (see Appendix E for a graphical representation of the model, and the original article for technical details). Note that in the initial article the number of classes of the mixture was directly deduced from the biological comparison to be performed. Based on this four component mixture model, posterior probabilities to belong to each class were computed for each probe. In the original article the MAP rule was applied to infer class memberships, and probes classified into one of the two classes of interest (under and over-methylated classes) were identified and further investigated. Starting from the same posterior probabilities as in the original publication, we performed probe classification into the classes of interest using the optimal classification rule defined in section 3. Here P = 4 and K = 2, and the MFDR is controlled at nominal level α = 0.1. Among the 15,801 probes initially classified as under or overmethylated, 13,065 are classified by the resulting MFDR classification rule, yielding an MFNR estimated at 0.0624.

The different classification rules are illustrated in Figure 4. The left panel represents the MAP rule used in the initial analysis (black = non-methylated, red = identically methylated, blue = under-methylated in leaf, green = over-methylated in leaf). The center panel corresponds to the optimal rule. Here the same colouring of the points is applied, with grey points corresponding to probes unclassified by the optimal rule and purple points corresponding to probes initially classified as over or under-methylated by the MAP rule but unclassified by the optimal rule.

As expected, all purple points are positioned on the boundaries between the two classes of interest and the two other classes. Although the posterior probability to be either overmethylated or under-methylated is higher than any other one for these probes, the actual value of the maximal posterior probability may be quite low, ranging between 0.25 and 0.6 (see Figure 4, right). In contrast, all probes classified by the optimal rule have posterior probabilities higher than 0.57. One can also notice that the optimal classification rule discards some obvious spurious classifications, such as the ones observed for very low values of the methylation signal, that were probably due to the constrained shapes of the covariance matrices of the adjusted model. Such points are discarded due to their low posterior probabilities.

In order to gain some additional insight regarding the boundary probes, we investigated to which extent the classification of a given probe is consistent with the ones of the adjacent (left and right) probes1 . For a given probe, the classification consistency score counts the number of adjacent probes sharing the same classification status. This score takes value 1: (a) Distribution of the probe counts in terms of MAP classification (rows) and classification consistency with neighbours (columns). The four classes correspond to "non-methylated" (1), "identically methylated" (2), "under-methylated in leaf" (3) and "over-methylated in leaf" (4). The classification consistency score counts the number of adjacent probes sharing the same classification status as the current probe. Numbers in brackets correspond to per row fractions. (b) Same distribution for the optimal classification. 0, 1 or 2, with 2 indicating a high classification consistency between the probe and its two neighbours.

Table 1 displays the distribution of probes in terms of class assignment and classification consistency, for the MAP rule (left) and the optimal rule (right). Focusing on class "over-methylated" one notices that the proportion of probes with score 1 or 2 is higher when using the optimal classification rule rather than the MAP classification rule. This illustrates the fact that most of the probes discarded by the optimal rule (corresponding to the purple points of Figure 4, center) are isolated probes exhibiting weak evidence for methylation. One can conclude that applying the optimal rule yields a more robust set of candidate probes by shaving inconsistent candidates.

Discussion

The methodology presented here builds on two central ideas: i) in many unsupervised settings only a subset of the classes is of interest for the practitioner, and ii) one would like to provide some guarantees (in terms of error rates) regarding the classification of observations into these specific classes of interest. To this end, classes of interest must be identified beforehand. In some contexts this identification is straightforward. In the methylation application of Section 4 both the a priori number of classes P and the number of classes of interest K are known, a situation that may arise whenever the application context corresponds to a differential analysis setting where two or more conditions are compared. Such examples correspond to the ideal application cases of MFDR control procedures. However, the procedure may also be directly applied to other contexts where P and/or K is unknown. Two examples of such applications are

• cases where P is unknown but there exists a clear H 0 class with known distribution, and the goal is to identify observations that do not belong to the H 0 class; here both P and K are initially unknown, and K can be set to P -1 once the fitting and selection of the model is performed. Note that in this context and for a fixed nominal level of MFDR, the MFNR could reach different values for different choices of P (hence K). In particular, one could expect the MFNR to increase with higher values of P.

Fixing a maximum level for the MFNR could then guide the model selection, by e.g. choosing P as the largest P satisfying the constraint -a larger P leading to a finer granularity of the classification.

• cases where both P and K are unknown but prior information about some observations belonging to a class of interest is available. For instance, in a genomics context one may identify classes of genes based on their expression profiles, then i) identify classes where some genes have a known biological function, and ii) apply the MFDR control procedure to classes corresponding to these functions.

In such application cases where either P or K are unknown, the MFDR procedure may be impacted by the model selection, something that was not investigated here. To what extend this impact will affect the procedure is difficult to quantify and will directly depend on the quality of the a priori knowledge one has at hand. For instance in the genomics application mentioned above, applying the MFDR procedure to a well-characterized class of genes may be highly efficient even if the number of classes P is poorly estimated or if some (other) classes are poorly fitted by the inferred model.

While the present paper focused on finding optimal classification rules in an unsupervised framework, a same motivation may exist in the context of (semi-) supervised classification. From a practical point of view, since the optimal classification rules derived in this article only depend on posterior probabilities, they can be applied to any statistical method that yields such probabilities, e.g. logistic regression or discriminant analysis, and could also be extended to methods for which pseudo posterior probabilities can be obtained (see [START_REF] Tao | Posterior probability support vector machines for unbalanced data[END_REF]Wang, 2005, Grandvalet, Mariéthoz, and[START_REF] Grandvalet | A probabilistic interpretation of svms with an application to unbalanced classification[END_REF]references therein).

Future work also includes the definition and derivation of optimal restricted classification rules in contexts where the labels are not assumed to be independent, such as hidden Markov models and/or latent variables models for network data (stochastic block model and latent block model, see [START_REF] Matias | Modeling heterogeneity in random graphs through latent space models: a selective review[END_REF] where existing results are restricted to the binary classification case [START_REF] Sun | Large-scale multiple testing under dependence[END_REF]. Additionally the consistency of the MFDR estimates presented here could be investigated using techniques similar to the ones developed in [START_REF] Denis | Consistency of plug-in confidence sets for classification in semi-supervised learning[END_REF], such that their theoretical results could be extended to the case where P > 2 and only some classes are of interest.

Appendix A: Proof of Proposition 1

We denote ψ the MAP rule amongst class 1 to K, and ψ R this rule restricted to a region R. If MNPR(ψ X ) ≤ α, since MFNR(ψ R ) ≥ MFNR(ψ X ) for any region R, obviously R * = X is a solution.

Otherwise, let R * be a region satisfying the conditions of the theorem, and let R be any region such that the MAP classification rule restricted to R satisfies MNPR(ψ R ) ≤ α = MNPR(ψ R * ). Let us prove that MFNR(ψ R ) ≥ MFNR(ψ R * ). Observing that R = R ∩ (R * ∪ R * ) and R * = R * ∩ (R ∪ R), one has :

P (ψ(X) = Z, X ∈ R) ≤ P (ψ(X) = Z, X ∈ R * ) ⇔ P (ψ(X) = Z, X ∈ R ∩ R * ) ≤ P (ψ(X) = Z, X ∈ R * ∩ R) ⇔ R∩ R * (1 -τ * K (x)) f (x)dx ≤ R * ∩ R(1 -τ * K (x)) f (x)dx,
where the last equation follows from the application of the MAP rule on R and R * respectively. Similarly,

P (Z ∈ {1, . . . , K}, X ∈ R) -P (Z ∈ {1, . . . , K}, X ∈ R * ) ≥ 0 ⇔ P (Z ∈ {1, . . . , K}, X ∈ R ∩ R * ) -P (Z ∈ {1, . . . , K}, X ∈ R * ∩ R) ≥ 0. Therefore P (Z ∈ {1, . . . , K}, X ∈ R ∩ R * ) = R∩R * K ∑ k=1 τ k (x) f (x)dx ≥ λ R∩R * (1 -τ * K (x)) f (x)dx ≥ λ R * ∩R (1 -τ * K (x)) f (x)dx ≥ R * ∩R K ∑ k=1 τ k (x) f (x)dx = P (Z ∈ {1, . . . , K}, X ∈ R * ∩ R) ,
where the first and third inequalities above follow from the property of R * (and thus R * ) given in equation ( 7).

Assuming that the marginal distribution of X is without atoms, we now establish that there exists a region R * satisfying (7) for some λ ≥ 0, and satisfying the MNPR constraint with equality, provided MNPR(ψ X ) > α. Define the event A = {ψ(X) = Z}, and Q the measure on the space X as

Q(R) = P({X ∈ R} ∩ A) = MNPR(ψ R ). Since the marginal distribution of X is without atoms, so is Q. Denote ξ (x) = (∑ K k=1 τ k (x))/(1 -τ * K (x))
, and define the level sets for any λ ≥ 0:

L(λ ) = {x ∈ X : ξ (x) ≥ λ }; L > (λ ) = {x ∈ X : ξ (x) > λ }.
Observe that it holds L > (λ ) ⊆ L(λ ) for any λ , and L(λ

) = λ <λ L > (λ ), L > (λ ) = λ >λ L > (λ ), so that lim sup λ λ Q(L > (λ )) = Q(L > (λ )) ≤ Q(L(λ )) = lim inf λ λ Q(L > (λ )). ( 12
) Define λ * = sup{λ ≥ 0 : Q(L > (λ )) ≥ α} ∪ {0},
Observe that it holds L < (λ ) ⊆ L(λ ) for any λ , and L(λ

) = λ >λ L < (λ ), L < (λ ) = λ <λ L < (λ ), so that lim sup λ λ Q(L < (λ )) = Q(L < (λ )) ≤ Q(L(λ )) = lim inf λ λ Q(L < (λ )). (13)
Define λ * = inf{λ ≥ 0 : Q(L < (λ )) ≥ 0}; we take λ * = ∞ if the above set is empty. By the definition of λ * , it holds

Q(L < (λ * )) ≤ 0 ≤ Q(L(λ * )); this is true by (13) if 0 < λ * < ∞; if λ * = ∞, it holds as well because Q(L(∞)) = Q(X ) > 0 since MFDR(ψ X ) > α; if λ * = 0, it still holds because Q(L < (0)) ≤ 0 in general. Therefore, since Q is without atoms, there exists a set R ⊆ X such that L > (λ * ) ⊆ R ⊆ L(λ * ) with Q(R) = 0.
This will imply MFDR(ψ R ) = α, as required, provided P (X ∈ R) > 0, which we still have to check to finish the proof.

Recall we assumed 0), this means P (X ∈ L(0)) > 0. If λ * > 0, it holds L(0) ⊆ R and we are finished. In the special case λ * = 0, it holds both Q(L(0)) ≤ 0 by definition of L(0), and Q(L(0)) ≥ 0 by definition of λ * . Hence Q(L(0)) = 0 in this case, and we can take R = L(0).

P (τ * K (X) ≥ 1 -α)) > 0. Since {x ∈ X : τ * K (X) ≥ 1 -α)} = L(
Appendix C: Illustration of different configurations of the simulation study q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 2 4 -2 -1 0 1 2 X1 X2 Label q q q 1 2 3 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 0 2 4 -4 -2 0 2 X1 X2 Label q q q 1 2 3 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -4 0 4 -4 0 4 X1 X2 Label q q q 1 2 3 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 1 2 scenario of "intermediate" difficulty, i.e. ν ν-2 σ 2 = 1. We considered several degrees of freedom ν (5, 10, 20, 50), and the same values of D as in our Gaussian simulations (0, 1, 2, 3). We generated 100 simulations for each configuration.

MFDR

Figure 9 shows the MFDR and MFNR of the considered rules using exact posterior probabilities in the case when all the classes are of interest, while Figure 10 shows error rates when the posterior probabilities are estimated using the EM algorithm for Gaussian mixtures implemented in mclust. Note that these results should be compared to the results shown in the central column of Figures 1 and3 (σ = 1) in the main text.

  1≤p ≤P τ p (x) and τ p (x) > 1α,

Figure 1 :

 1 Figure 1: Performances of the MAP, optimal and 1α thresholded classification rules in terms of realized MFDR and MFNR. Columns correspond to different values of the variance σ 2 . The black line on the top graphs corresponds to the nominal level of 5%. This figure appears in color in the electronic version of this article, and color refers to that version.

Figure 2 :

 2 Figure 2: Boxplots of the estimated threshold for the optimal rule. Colors correspond to levels of the variance σ 2 . The horizontal line indicates the threshold 0.95 used in the thresholded rule. This figure appears in color in the electronic version of this article.

Figure 3 :

 3 Figure 3: Same figure as Figure 1, except that posterior probabilities are now estimated. This figure appears in color in the electronic version of this article.

Figure 4 :

 4 Figure 4: Left: Probes of chromosome 4 coloured according to the MAP classification rule. Center: Probes of chromosome 4 coloured according to the optimal classification rule. Right: Histogram of the maximum posterior probabilities for purple probes classified by the MAP rule but not by the optimal rule. This figure appears in color in the electronic version of this article.

Figure 5 :

 5 Figure 5: Three examples of simulated data, with different parameter values: an easy case (left) corresponding to D = 3 and σ 2 = 0.5, and intermediate case (center) corresponding to D = 2 and σ 2 = 1 and a difficult case (right) corresponding to D = 0 and σ 2 = 2. Colors correspond to class labels.

Figure 6 :

 6 Figure6: Performances of the MAP, optimal and 1α thresholded classification rules in terms of estimated and true MFDR and MFNR. All rules are based on the true posterior probabilities. Only classes 1 and 3 are assumed to be of interest.

Figure 9 :

 9 Figure 9: Performances of the MAP, optimal and 1α thresholded classification rules on mixture of three bivariate Student distributions. All rules are based on the true posterior probabilities. Columns correspond to different values of the degrees of freedom ν. All classes are considered of interest.

Although not accounted for in the initial publication, the information about the genomic position of the probes is available.

Acknowledgments

G. Blanchard acknowledges support from Agence Nationale de la Recherche (ANR) via the project ANR-19-CHIA-0021-01 (BiSCottE), and the project ANR-16-CE40-0019 (SansSouci); and from the Franco-German University through the binational Doktorandenkolleg CDFA 01-18. GQE and IPS2 benefit from the support of the LabEx Saclay Plant Sciences-SPS (ANR-17-EUR-0007).

By the definition of λ * , it holds

this is true by (12) if λ * > 0, and, if λ * = 0, it holds as well since we assumed Q(L(0)) = Q(X ) = MNPR(ψ X ) > α. (It can be also checked that it still holds if λ * = ∞, which can happen in principle when P (ξ (X) = ∞) > 0.) Therefore, since Q is without atoms, there exists a set R ⊆ X such that L > (λ * ) ⊆ R ⊆ L(λ * )) with Q(R) = MNPR(ψ R ) = α, as required.

Appendix B: Proof of Proposition 2

Thus, in this case only a region R with null probability (in particular R * = / 0) satisfies the MFDR constraint since by convention MFDR(Ψ / 0 ) = 0. Otherwise, let R * be a region satisfying (10) and the MFDR constraint with equality. Consider R any region satisfying

On the one hand, one has :

On the other hand,

Analogously to the proof of Proposition 1, we establish the existence of R * satisfying the required conditions, provided we exclude the edge cases, i.e. we assume that MFDR(ψ X ) > α, and also

Let Q be the signed measure on the space X defined as

, and define the level sets for any λ ≥ 0:

Appendix D: Analysis of simulated data when only some classes are of interest.

Here the same simulation setting is considered as in the main article, the only difference being that it is now assumed that only classes 1 and 3 are of interest. Figure 6 displays the performance of the different classification rules when the true parameters of the model are known. Figure 7 displays the same result when the parameters of the model are inferred from the data.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.5 1 2 Appendix F: Simulations based on a mixture of Student bivariate distributions

MFDR

For this supplementary simulation study, datasets were simulated from a mixture of three bidimensional Student distributions. We recall that given a Gaussian vector y ∼ N (0, Σ), a location vector µ and u ∼ χ 2 ν , the vector x = y u ν + µ follows a multivariate Student distribution with parameters Σ, µ, ν. For each of the three components in the mixture, we simulated 200 observations with the same locations µ 1 = (-1, 0), µ 2 = (0, D) and µ 3 = (1, 0) as in our principal Gaussian simulations. We took the scale matrix Σ to be the diagonal matrix σ 2 I where σ was chosen so that the covariance matrix ν ν-2 Σ is the same as the one in the Gaussian


Sébastien Légaré

Jean Krivine

Russ Harmer

Jérôme Feret

Modelling Systems Biology Wide and Deep ⋆

Keywords: Systems biology, Rule-based modelling, Causality

Biological systems involve a wide amount of different molecular interactions. Each interaction can in turn present a deep level of mechanistic details. Most modelling methods can encompass only one of those two dimensions, that is the width or depth of biological systems. Rule-based modelling is a powerful method in that regard as it can model large systems with mechanistically detailed interactions. The advantage of including mechanistic details is that they allow a clear tracking of causal chains of events between biomolecules. This translates into a precise pinpointing of upstream regulators and a better understanding of complex biological systems.

Introduction

Interaction networks [START_REF] Barabási | Network medicine: a network-based approach to human disease[END_REF] are one of the most popular objects of systems biology. They allow the representation of a very large amount of molecular interactions on a single graph. That is, interaction networks can encompass the width of biological systems. However, many functional relationships between biomolecules are more complicated than what can be represented in simple graphs [START_REF] Klamt | Hypergraphs and Cellular Networks[END_REF]. Such networks can hence not depict in detail the depth of molecular interactions.

A consequence of this lack of detail is that wrong interpretations can occur when analyzing interaction networks. For example, consider the analysis of upstream regulators of a chosen molecule. In intricate graphs like signalling networks, it is not uncommon to have every output downstream of every input [START_REF] Rowland | Crosstalk and the evolvability of intracellular communication[END_REF]. In these conditions, ignoring mechanistic and quantitative details usually leads to the perplexing result that the chosen molecule is regulated by all the other molecules in the network with indistinguishable intensities.

To include more details in the description of molecular interactions, reactionbased models like Petri nets [START_REF] Heiner | Petri Nets for Systems and Synthetic Biology[END_REF] can be used. Such methods properly address the depth of molecular interaction mechanisms. However, this comes at the cost of a combinatorial explosion as the system to model grows. Very large models of biological organisms can hence not be built in a reaction-based setting.

Most organisms of interest to systems biology are large in terms of number and complexity of their constituent interactions. Methods are hence needed to deal with both the width and depth of biological systems simultaneously. Rulebased modelling [START_REF] Danos | Rule-Based Modelling of Cellular Signalling[END_REF][START_REF] Blinov | Graph Theory for Rule-Based Modeling of Biochemical Networks[END_REF] can achieve such feat.

In this work, we present a large rule-based model example. Then, we show how causality analysis can be used to extract pathways to observables of interest. Unlike upstream analysis in interaction networks, those pathways leverage mechanistic details to provide a real causal explanation of molecular events.

Rule-Based Model

Rules distinguish themselves from reactions by focusing only on the parts of molecules that change during a transition, rather than fully defining the species involved. It allows the representation of mechanistic details while avoiding the vast majority of combinatorial issues. The Kappa [START_REF] Boutillier | The Kappa platform for rule-based modeling[END_REF] rule-based modelling language was used in this work.

Causality Analysis

Causality analysis can be used to extract the upstream events which cause any observable of interest from a rule-based model, like the activation of a chosen protein. Since this analysis takes mechanistic details into account, it provides a clearer and more reliable picture that upstream analysis in interaction networks.

The first step of causality analysis is to extract stories to the observable of interest. A story consists in a graph that represents all the necessary events that lead to a single instance of the observable of interest during an execution of the model. Stories are extracted using either operational research techniques [START_REF] Danos | Graphs, Rewriting and Pathway Reconstruction for Rule-Based Models[END_REF] or heuristic approaches [START_REF] Laurent | KaFlow[END_REF].

1 SRC (sh2[.]) FYN (sh2[.]) 2 SRC (act1{t}) FYN (act1{t}) 3 BMX (Y566[1]), SRC (tyr_kin[1]) BCR (Y177[1]), FYN (tyr_kin[1]) PTK2B (Y580[1]), FYN (tyr_kin[1]) 4 BMX (Y566{p}) BCR (Y177[.]{p}) PTK2B (Y580{p}) 5 BMX (act1{t}) BCR (Y177[1]), FES (sh2[1]) PTK2B (act{t}) 6 BMX (Y216{p}) FES (act{t}) SRC (Y419[1]), PTK2B (tyr_kin[1]) 7 BMX (act2{t}) SRC (Y419{p}) BCR (Y276[1]
), FES (tyr_kin [START_REF] Barabási | Network medicine: a network-based approach to human disease[END_REF])

BCR (Y177[1]
), FES (tyr_kin [START_REF] Barabási | Network medicine: a network-based approach to human disease[END_REF])

BCR (Y283[1]
), FES (tyr_kin [START_REF] Barabási | Network medicine: a network-based approach to human disease[END_REF])

BCR (Y246[1]
), FES (tyr_kin [START_REF] Barabási | Network medicine: a network-based approach to human disease[END_REF])

BCR (Y279[1]), FES (tyr_kin[1]) 8 PTK2 (Y576[1]), BMX (tyr_kin[1]) SRC (act2{t}) BCR (Y276[.]{p}) BCR (Y283[.]{p}) BCR (Y246[.]{p}) BCR (Y279[.]{p}) 9 PTK2 (Y576{p}) BCR (Y276[1]), ABL1 (sh2[1]
)

SRC (Y419[1]), SRC (tyr_kin[1]
)

BCR (Y283[1]), ABL1 (sh2[1]
)

BCR (Y246[1]), ABL1 (sh2[1]
)

BCR (Y279[1]), ABL1 (sh2[1]) 10
PTK2 (act{t}) ABL1 (act{t})
11

PXN (Y118[1]), PTK2 (tyr_kin[1]
)

SRC (Y419[1]), PTK2 (tyr_kin[1]) TP73 (Y99[1]), ABL1 (tyr_kin[1]) EGFR (Y1197[1]), ABL1 (tyr_kin[1]) PXN (Y118[1]), ABL1 (tyr_kin[1]) BCR (Y177[1]), ABL1 (tyr_kin[1]) 12 PXN (Y118{p}) TP73 (Y99[.]{p}) EGFR (Y1197[.]{p}) 13 PXN-y118p TP73 (Y99[1]), ABL1 (sh2[1]) EGFR (Y1197[1]), ABL1 (sh2[1]) Fig. 2.
Example of a pathway obtained from causality analysis of the pYnet model. Node labels contain protein names in bold. Sites modified at each step are written following the Kappa syntax [START_REF] Boutillier | The Kappa platform for rule-based modeling[END_REF] in normal font. The width of edges represents their quantitative contribution.

Then, all the stories to the observable of interest are folded into a single quotient graph. Proper folding requires finding the correct trade-off in the amount of context to consider from story events [START_REF] Légaré | Distinguishing Context Dependent Events in Quotients of Causal Stories[END_REF]. Too much contextualization leads to unnecessarily large pathways while too few looses information. The contextualization and folding steps are carried by the Python package KappaPathways [START_REF] Légaré | KappaPathways[END_REF]. The resulting graph shares a lot of similarities with biological pathways.

Fig. 2 presents an example of a pathway obtained from causality analysis on the pYnet model. It allows a clear understanding of the steps which lead to the phosphorylation of tyrosine 118 from protein Paxillin (PXN). In comparison, going upstream from the node representing PXN in an interaction network would not have allowed a distinction between the structural features involved at each backward step. It would hence have most likely resulted in a much larger graph.

Conclusion

Rule-based modelling allows the construction of wide and deep models of biological systems. Causality analysis of those models provides detailed explanations for the occurrence of molecular events. Combining rule-based modelling with causality analysis could be the best method to truly understand complex biological systems.

Fig. 1 .

 1 Fig. 1. Graph representation of the pYnet model built with Kami [7]. On the left is a contact map of the model, which bears resemblance to interaction networks. The inset on the right zooms in and highlights the structural details of selected proteins.

Fig. 1

 1 Fig.1shows a graph representation of a rule-based model which goes both wide and deep in its description of a biological system. Dubbed the pYnet, this model depicts tyrosine phosphorylation in human cell signalling. It contains 175 different proteins and about 1800 interactions taken from the literature and databases. For each interaction, the precise protein residues and domains involved are specified.

⋆ Supported by French ANR project DCore ANR-18-CE25-0007.