Benjamin Hivert
email: benjamin.hivert@u-bordeaux.fr

Denis Agniel

Rodolphe Thiébaut

Boris P Hejblum

Post-clustering difference testing: valid inference and practical considerations

Keywords: Clustering, hypothesis testing, double-dipping, circular analysis, selective inference, multimodality test, Dip Test

Clustering is part of unsupervised analysis methods that consist in grouping samples into homogeneous and separate subgroups of observations also called clusters. To interpret the clusters, statistical hypothesis testing is often used to infer the variables that significantly separate the estimated clusters from each other. However, data-driven hypotheses are considered for the inference process, since the hypotheses are derived from the clustering results. This double use of the data leads traditional hypothesis test to fail to control the Type I error rate particularly because of uncertainty in the clustering process and the potential artificial differences it could create. We propose three novel statistical hypothesis tests which account for the clustering process. Our tests efficiently control the Type I error rate by identifying only variables that contain a true signal separating groups of observations.

Introduction

Cluster analysis is ubiquitous in medical research (see [START_REF] Gj Mclachlan | Cluster analysis and related techniques in medical research[END_REF] for a comprehensive overview) to perform data classification, data exploration, and hypothesis generation [START_REF] Xu | Clustering[END_REF]. Clustering works by grouping homogeneous observations into disjoint subgroups or clusters. When multivariate data are clustered, it is common to seek to identify which variables distinguish two or more of the estimated clusters, in order to interpret the clustering structure and characterise observation groups and how they differ from each other.

Despite the widespread use of clustering, [START_REF] Hennig | Handbook of cluster analysis[END_REF] state there is no commonly accepted and formal definition of what clusters are. In fact, the definition of what a cluster should be varies depending on the context and the analysis specifics. Here we will use the definition from Everitt and Hothorn [2006], which includes only two criteria: i) homogeneity of observations within a cluster and ii) separability of observations between two different clusters. These two criteria are general enough to encompass the majority of the working definitions of clusters. Both can be quantified using various approaches such as distances or similarity metrics, shape of distribution [START_REF] Steinbach | The challenges of clustering high dimensional data[END_REF], multimodality [START_REF] Kalogeratos | Dip-means: an incremental clustering method for estimating the number of clusters[END_REF]Likas, 2012, Siffer et al., 2018], or distributional assumptions [START_REF] Liu | Statistical significance of clustering for high-dimension, low-sample size data[END_REF][START_REF] Patrick K Kimes | Statistical significance for hierarchical clustering[END_REF].

While clustering is a multivariate methodology that takes into account all variables, only a set of variables can be expected to differentiate two particular clusters (i.e. separate their observations, according to the second criterion of our definition above). This question of which variable separate clusters of individuals is particularly relevant for high-dimensional data such as omics data [Ntranos et al., 2019, Vandenbon and[START_REF] Vandenbon | A clustering-independent method for finding differentially expressed genes in single-cell transcriptome data[END_REF]. The current practice to identify such variables is often based on post-clustering hypothesis testing. It leads to a two-step pipeline (a first step of clustering and a second step of inference) that is actually testing data-driven hypotheses in a process sometimes referred to as "double dipping" [START_REF] Kriegeskorte | Circular analysis in systems neuroscience: the dangers of double dipping[END_REF]. This approach does not efficiently control the type I error rate when testing for differences between clusters. In fact, it is always possible to cluster the data using a clustering method, even if there is no real process separating groups of observations. In this case, the clustering method artificially enforces the differences between the observations by dividing them into different clusters. The significant differences between clusters identified during the inference process could just be an artifact of the previous clustering step. To illustrate this phenomenon, we consider data generated from a univariate Gaussian distribution with mean 0 and variance 1 (Figure 1 panel A). Two clusters can be built, e.g., using hierarchical clustering with Ward's method and Euclidean distance (Figure 1 panel B). These two estimated clusters are not separated clusters, since all observations come from the same Gaussian distribution. One way to infer their separation is to test for a mean shift between them, for example using the classical t-test. Since there is no real process separating these two clusters, the resulting p-values should be uniformly distributed. However, when we look at the p-values of the t-test for 2000 simulations of the data, the resulting p-values are too small, leading to false positives (Figure 1 panel B). This simple example illustrates how it is possible to infer a separation of two clusters, even if this separation is not explained by a real process in the data. Classical inference requires a priori hypothesis. In this toy example, the hypothesis, i.e the lack of separation of the two clusters, is based on clusters derived from the data. Moreover, here we force differences between groups of observations by clustering them, so the clustering results do not represent the true structure of the data. Thus, the discoveries are only the results of clustering algorithms and not those of a true biological signal due to this double use of the data and the bad structures forced by clustering. For example, in the context of RNA-seq data analysis, accounting for this clustering step during the inference step is one of the open problems in the eleven grand challenges in single-cell data science mentioned by [START_REF] Lähnemann | Eleven grand challenges in single-cell data science[END_REF].

Our goal is to propose new methods for post-clustering inference that take into account the clustering step and the potential artificial differences it may introduce. For any clustering method that can be applied to all features of the data to build clusters, we are interested in testing the null hypothesis that a particular feature does not truly separate two of the estimated clusters. In particular, this null hypothesis allows that the feature: i) is not involved in the separation of the two subgroups and is not affected by the clustering step, and ii) is only involved in this separation because the clustering method applied to the data forced differences.

Recently, some methodological work has been done on post-clustering inference. Since the data is used twice, many of them use selective inference [START_REF] Ryan J Tibshirani | Exact post-selection inference for sequential regression procedures[END_REF][START_REF] Lee | Exact post-selection inference, with application to the lasso[END_REF] to account for the clustering step. Selective inference aims to control the selective type I error. This is defined as the probability 2 [START_REF] Gao | Selective inference for hierarchical clustering[END_REF] for separating the two estimated clusters. The uniform distribution is also shown for comparison under the null of rejecting the null hypothesis, given that the model and the null hypothesis have been selected thanks to the data. When data splitting is not possible, [START_REF] Fithian | Optimal inference after model selection[END_REF] has proposed to condition on this selection event during statistical hypothesis testing. In applying this approach, we use two different types of data: the data, to construct the model and the hypothesis, and the data given the fact that it has been used, i.e. data not yet observed, to perform the test. This leads to statistical hypothesis tests that efficiently control the selective type I error. Selective inference was first proposed for linear regression, change point detection [START_REF] Jewell | Testing for a change in mean after changepoint detection[END_REF] and more recently for tree regression [START_REF] Neufeld | Tree-values: selective inference for regression trees[END_REF]. Clustering is also a framework in which selective inference has been applied recently. For post-clustering inference applied on RNA-seq data, [START_REF] Zhang | Valid post-clustering differential analysis for single-cell rna-seq[END_REF] have developed a truncated-normal statistic that use selective inference and leads to valid p-values under their null of no differential expression. However, in addition to selective inference, they use data splitting, which is only possible if the number of observations is large enough. They also use a supervised approach to predict the partition formed on half of the data on the remaining half. Instead of conditioning on the clustering event in their statistical hypothesis test, they condition on the fact that in the remaining half of the data, the labels of the observations are predicted thanks to a supervised approach. More recently, [START_REF] Gao | Selective inference for hierarchical clustering[END_REF] have developed a multivariate selective test to investigate whether two estimated clusters are truly separated or whether the observations they contain come from a single cluster. By using selective inference, they account for the clustering step. Their approach is suitable for cluster validation because their null hypothesis is the equality of two cluster centers. This method also leads to valid p-values under the null hypothesis (Figure 1 panel c). However, this method is not suitable for our purpose, since in this particular context the goal is to study the separation of two clusters at the feature level, i.e., in a univariate setting.

In this paper, we introduce three new methods for post-clustering inference. First, we adapt the method proposed by [START_REF] Gao | Selective inference for hierarchical clustering[END_REF] for univariate hypotheses to investigate whether individual features contain information about group (clustering) structure. In doing so, we use a data-driven and fixed clustering of the data to ensure interpretations. To deal with the multiple clusters case, we also present an extension of this first test based on an aggregation of its p-values. Second, we propose another approach using a test of multimodality that account for the clustering step by investigating the presence of a continuum in the distribution of the variable. The paper proceeds as follows. In the Methods section, we describe the methods we proposed for post-clustering inference. These approaches are then evaluated and compared in the Results section using extensive numerical simulations and a real ecological dataset. Some final comments can be found in the Discussion section.

Methods

In the following, let X be a n × p random variable of n observations of p features, with g th column X g . On X we apply a clustering method c() to create c(X), a partition of the n observations into K disjoint clusters C 1 , . . . , C K . We are interested in the ability of a given variable X g to separate two clusters C k and C l estimated using all the information contained in X with the clustering method c().

Selective test

To develop our statistical hypothesis testing, we first specify a generative model to the observations along X g . We assume that each of the n observations of X g comes from independent Gaussian distributions with unknown mean µ gi and known variance σ 2 g . Then, for all i ∈ {1, . . . , n}, X gi ∼ N (µ gi , σ 2 g ). Because of the independence between each X gi , the multivariate distribution of X g is a multivariate Gaussian distribution N n µ g , σ 2 g I n with mean µ g = (µ g1 , . . . , µ gn ) t and covariance matrix Σ = σ 2 g I n . Let x g be the realisation of X g observed in X. Now, for a cluster C k , let

µ C k g = 1 |C k | i∈C k µ gi and X C k g = 1 |C k | i∈C k X gi
be the true mean and empirical mean, respectively, of the variable X g in cluster C k . Testing for a mean shift between the two clusters is a straightfoward way to evaluate the separation of two clusters along X g . Thus, we define the two following hypotheses:

H 0 : µ C k g = µ C l g vs H 1 : µ C k g = µ C l g (1)
By introducing a contrast vector η ∈ R n defined by: η i = [START_REF] Gao | Selective inference for hierarchical clustering[END_REF], we can rewrite (1) above as:

1 i∈C k |C k | - 1 i∈C l |C l | ∀i = 1, . . . , n following Jewell et al. [2019],
H 0 : µ t g η = 0 vs H 1 : µ t g η = 0 (2)
H 0 in (2) is actually generated by a function of the data c(X), which clearly sets us in the context of selective inference. Conditioning on this clustering event within statistical inference procedures is thus required. In particular, we derive an adaptation of the p-value proposed by [START_REF] Jewell | Testing for a change in mean after changepoint detection[END_REF] (originally intended for change point detection) for our purposes of clustering:

p C k ,C l g ≡ P H0 |X t g η| > |x t g η| | C k , C l ∈ c(X) (3)
Here we condition on the estimation of C k and C l by c(X), which leads to the definition of H 0 , and the resulting p-values (3) account for the clustering as well as the uncertainty associated with the estimation of these two clusters. p C k ,C l g quantifies the probability that the mean difference between C k and C l is as large as the observed difference under H 0 given the observed clustering structure. Its calculation relies on all possible realisations of X g resulting in the same estimation of C k and C l when we apply c() to X. Yet, enumerating all such data sets X is hard. To get more tractable p-values, we follow [START_REF] Jewell | Testing for a change in mean after changepoint detection[END_REF] and [START_REF] Gao | Selective inference for hierarchical clustering[END_REF] in constraining the randomness in the random variable X g and we define our p-value as follows:

pC k ,C l g ≡ P H0 |X t g η| > |x t g η| |C k , C l ∈ c(X), π η ⊥ X g = π η ⊥ x g (4)
where

π ⊥ η = I n -ηη t η 2 2
restricts the random variable X g to a space defined by the scalar π η ⊥ x g without losing control of type I error [START_REF] Gao | Selective inference for hierarchical clustering[END_REF]. The p-value (4) can be rewritten as (see Supplementary Materials for the proof):

pC k ,C l g = P H0 |φ g | > |x t g η| | φ g ∈ S g (5)
where

S g = {φ g : C k , C l ∈ c (x (φ g ))}
is the set of perturbations of the g th variable from X where both C k and C l are conserved by c(), and φ g = X t g η H0 ∼ N 0, σ 2 g η 2 2 . X(φ g ) thus represents a perturbed version of the data X, where only the g th variable is perturbed:

x g - ηη t x g η 2 2 + ηφ g η 2 2
This perturbation has a clear interpretation: if |φ g | > |x t g η| data from the two clusters are split further apart along X g than is observed in the data; whereas if |φ g | < |x t g η| instead, they are brought closer together along X g (and if φ g = x t g η the data are actually not perturbed because in this case X(φ g ) = X). Note that (5) can be rewritten as P H0 |φ g | > |x t g η|, φ g ∈ S g /P H0 (φ g ∈ S g ). So if C k and C l can only be preserved when the observation are perturbed further apart, then (5) will be large since P H0 |φ g | > |x t g η|, φ g ∈ S g P H0 (φ g ∈ S g ). In conclusion, this selective test can be interpreted in terms of separability of the two clusters considered (even though it is based on a difference in means) as it boils down to quantifying the possibility to bring closer together the observations from the two clusters while preserving their separation.

In order to explicitly describe the set S g while retaining as much generality as possible about c(), we follow [START_REF] Gao | Selective inference for hierarchical clustering[END_REF] and use Monte-Carlo simulations to approximate pC k ,C l g . This strategy relies on (5) being rewritten as:

pC k ,C l g = E 1 |φ g | > |x t g η|, φ g ∈ S g E [1 {φ g ∈ S g }] (6)
Namely, we sample φ 1 g , . . . , φ

N g i.i.d ∼ N 0, σ 2 g η 2 2
for some large value N , and replace the expectations in (6) with the sums over all samples. This Monte-Carlo procedure avoids the need to formally describe S g . In order to enhance numerical efficiency, [START_REF] Gao | Selective inference for hierarchical clustering[END_REF] use an importance sampling approach originally proposed by [START_REF] Yang | Selective inference for group-sparse linear models[END_REF] to improve the likelihood of preserving the clustering in the perturbed data.

Our proposed estimation of pC k ,C l g is thus:

pC k ,C l g ≈ N i=1 π i 1 |ω i g | ≥ |x t g η|, C k , C l ∈ c(X(ω i g )) N i=1 π i 1 C k , C l ∈ c(X(ω i g )) (7)
where ω 1 g , . . . , ω N g ∼ N x t g η, σ 2 g η 2 2 , and π i = f1(ω i g ) f2(ω i g ) the importance sampling probabilities with f 1 the density of a N 0, σ 2 g η 2 2 distribution and f 2 the distribution of a N x t g η, σ 2 g η 2 2 distribution. Of note, we adapt the method from [START_REF] Phipson | Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn[END_REF] to obtain unbiased Monte-Carlo p-values estimations (see Supplementary Materials for details).

At the core of the above test is the scaling variance parameter σ 2 g , which represents the variance of each column of X g . While σ 2 g is assumed to be known in the test, it is not the case in practice and we propose to use the following plug-in estimate instead:

σ2 g = 1 |C k | + |C l | -1 i∈C k ,C l X gi -X C k ,C l g 2 with X C k ,C l g = 1 |C k | + |C l | i∈C k ,C l X gi
This variance estimate only takes into account observations from the two clusters of interest in the test, in line our null hypothesis of no separation of the two clusters (the variance itself can informs on the separation of the data [START_REF] Liu | Understanding of internal clustering validation measures[END_REF]). In some cases, this σ2 g could underestimates the variance of X g (particularly if the clustering induces strong artificial differences). Meanwhile, [START_REF] Gao | Selective inference for hierarchical clustering[END_REF] rely on a different estimate, that instead overestimate the variance in certain cases. Still they have showed type I error control is guaranteed, even with an overestimated variance, at the cost of being overly conservative (see Supplementary Figure S1 for additional details).

A more powerful test in the presence of intervening clusters

The above selective test has been designed for comparing a pair of clusters. Yet, in practice there are often more than 2 clusters. In such case, the test could have very limited statistical power even for well separated clusters: if there is one or more additional cluster in-between the two clusters of interest, it quickly becomes impossible to perturbed them closer together without changing the clustering (see Supplementary Figure S2). To overcome this limitation, we extend the above selective test assuming that two estimated clusters C k and C l are separated on X g if and only if at least one of the adjacent cluster pairs in-between them are separated. This means that on the contrary, if there is a continuum on X g to go from C k to C l , then C k to C l are not separated. By testing only the separation of pairs of adjacent clusters we retain the statistical power of the selective test. We propose to combine all the in-between adjacent pair selective test p-values into a so called combined selective test to finally assess the separation of C k and C l on X g .

To identify the clusters in-between C k and C l on X g , we define the set:

C k:l g := C i , i = 1, . . . , K / X Ci g ∈ min X C k g , X C l g , max X C k g , X C l g
where

X Ci g = 1 |Ci| j∈Ci
X gj , implicitly sorting clusters according to their empirical mean on X g . We define two clusters C m1 and C m2 as adjacent on

X g if C m1:m2 g = {C m1 , C m2 }. So if |C k:l g | = M
, there are M -1 pairs of adjacent clusters in C k:l g that draw a path from C k to C l . We can then define:

p C k :C l g := f p 1 g , . . . , p M -1
g where f must be a merging function as described in [START_REF] Vovk | Combining p-values via averaging[END_REF]. Based on numerical simulations, we favor the use of the harmonic mean merging function, recommended by [START_REF] Vovk | Combining p-values via averaging[END_REF] for potential dependencies between the p-values -which is our case here since each cluster data contributes to two p-values -and features a good trade-off between type I error and statistical power (see Supplementary Figure S3).

Thus, we use:

p C k :C l g = min     e log (M -1) M -1 M -1 i=1 1 p i g , 1    
Of note, in order for all p 1 g , . . . , p M -1 g to be computed using the same variance estimate, we propose this time to plug-in an estimate of σ 2 g that accounts for all observations belonging to either one the adjacent clusters in C g:l g :

σ2 g = 1 |C k:l g | -1 C∈C k:l g i∈C X gi -X C k:l g g 2 with X C k:l g g = 1 |C k:l g | C∈C k:l g i∈C X gi

Multimodality test

The separation of two clusters according to a given variable is equivalent to this variable's distribution being multimodal. Following [START_REF] Kim | Marcopolo: a clustering-free approach to the exploration of differentially expressed genes along with group information in single-cell rna-seq data[END_REF], multimodality thus becomes a marker for the separation of clusters: each mode corresponds to a group of homogeneous observations (i.e., a cluster), separated by less dense regions of the distribution. But as with artificial mean differences arising from clustering, multimodality may also be an artefact caused by the clustering method c(). We propose to leverage this notion of continuum between two clusters: if C k and C l are separated, then there must be dip in the distribution of this variable at some point between the two (i.e. multimodality). On the other hand, if there is a continuum between these two clusters, then they cannot truly be separated (i.e. unimodality). Fortunately, such a continuum cannot be caused by the clustering method. This second proposal can be seen as a simplified version of our first selective test. Indeed, by perturbing the data in the selective test to see if we can bring the two clusters closer without changing the clustering, we assess how likely it would be to observe a continuum between C k and C l . If there is a continuum between C k and C l on X g , then its distribution must be unimodal. Thus, to investigate separability of those two clusters on X g , it suffices to apply a unimodality test to its distribution restricted only to the individuals from clusters of the set C k:l g . Indeed, if the X g separates C k and C l , then there are at least two clusters in C k:l g that are separated from each other, and in particular, since these clusters are between C k and C l , there is also a separation between them on X g .

A unimodality test compares the null hypothesis "distribution of X g is unimodal" to the alternative "distribution of X g is multimodal". In the context of unsupervised clustering, [START_REF] Kalogeratos | Dip-means: an incremental clustering method for estimating the number of clusters[END_REF] developed a clustering algorithm based on incremental unimodality testing, and [START_REF] Siffer | Are your data gathered[END_REF] developed unimodality a test to assess data clusterability based on their multivariate distribution. Ameijeiras-Alonso et al. [2021] give a recent overview on unimodality testing, but three tests are the most frequent: i) the Silverman test [START_REF] Bernard W Silverman | Using kernel density estimates to investigate multimodality[END_REF] based on the kernel estimate of the density f of the data, ii) the Dip Test [START_REF] John A Hartigan | The dip test of unimodality[END_REF] based on the cumulative distribution function F , and iii) the excess mass test [START_REF] Werner | Excess mass estimates and tests for multimodality[END_REF]. The Dip Test avoids the need for estimating of additional parameters or making any distributional assumption and has already been applied to clustering [START_REF] Kalogeratos | Dip-means: an incremental clustering method for estimating the number of clusters[END_REF][START_REF] Wasserman | Feature selection for high-dimensional clustering[END_REF][START_REF] Schelling | Dataset-transformation: improving clustering by enhancing the structure with dipscaling and diptransformation[END_REF]. Furthermore, compared to several multimodality tests available in the R package multimode Ameijeiras-Alonso et al. [2021], the Dip Test outperforms its competition both in terms of computation times and performances (see Supplementary Figure S4).

The Dip Test from [START_REF] John A Hartigan | The dip test of unimodality[END_REF] relies on the dip statistic dip(F ) = min [1985] showed that the Uniform distribution is the unimodal distribution with the asymptotically largest dip statistic among the unimodal distributions (intuitively the least favourable candidate for unimodality): a distribution with a dip statistic larger than that of the uniform distribution cannot be unimodal. Thus p Dn is interpreted as the probability under the null case of unimodality that the uniform distribution has a dip statistic greater than the observed dip statistic of F n .

For our purposes, we apply the Dip Test to the distribution of the variable X g restricted to the individuals that are in the clusters of the set C k:l to test for a continuum between C k and C l .

Results

Numerical simulations study

We present here results evaluating the behaviour of our proposed tests in the Methods section both in terms of type-I error control and statistical power.

Behaviour in a two-dimensional setting

We generated two-dimensional data (p = 2) under two scenarios: (i) first with no separated clusters from a common standard Gaussian distribution N (0, 1); and (ii) second with three clusters from Gaussian distributions N (µ Cj , 1) and built-in mean differences µ C1 = (-5, 0), µ C2 = (5, 0), and µ C3 = (0, 10) (thus X 1 separated all three clusters while X 2 only separated C 3 from the rest, meaning X 2 was under the null when comparing C 1 and C 2 ). In both cases, we applied hierarchical clustering with Ward method and Euclidean distance to build three clusters. Figure 2A shows an example realisation for each scenario. In the first scenario, clusters were estimated by forcing differences between groups of observations, creating artificial differences between clusters, while in the second scenario, the estimated clusters represented the true structure of the data (each cluster weas a homogeneous and separate group of observations).

Figure 2B shows the results of the three proposed approaches compared to the p-values from the usual t-test for 2,000 repeated simulations each with a sample size of n = 200. For the no cluster scenario, the t-test yielded extremely small p-values which translates into a direct inflation of the type-I error. In fact, the t-test identified the artificial differences created during the clustering process. Taking into account this clustering step, the p-values of the selective test p C k ,C l g and the p-values resulting from its merging extension p C k :C l g were fairly uniformly distributed over [0, 1], ensuring a good calibration of the p-values and a control of type-I error. As for the multimodality test, its p-values were overly conservatives but consistent with no real separation of clusters. This was due to its reference being the Uniform distribution (the limit case for a unimodal distribution) while the data were generated from a Gaussian distribution (which has a lower dip statistic than the uniform distribution). Those good results were confirmed under the 3 clusters scenario when comparing C 1 and C 2 along X 2 . For all other comparisons under this scenario, all 4 tests correctly detected the separated clusters that are under H 1 .

Of note, if clustering does not artificially force differences between groups of observations, e.g by discovering the actual group structure in the data, the t-test also control the type-I error. This illustrates the connections between artificial differences and the estimation of the number of clusters. But, this process is still testing data driven hypothesis, which does not respect the classical inference setting where hypothesis must be specified without using the data.

Statistical power

We now generated data from a univariate mixture of two Gaussian distributions with equal proportions and variance: 0.5N (0, 1) + 0.5N (δ, 1), where the two components were separated by a mean difference of δalso called a contamination model [Laurent et al., 2018]. The two components distinguished two different clusters, and the magnitude of δ tuned their separability. We applied hierarchical clustering with Ward method and Euclidean distance to build either 2 or 4 clusters. Figure 3A displays an example realization of this simulation. In the 2-cluster case, the true structure of the data was uncovered, while in the 4-cluster case spurious clusters were introduced. We evaluated the statistical power of the three proposed tests to detect the separation between clusters 1 and 2, the two most extreme clusters in the distribution of the data, according to δ at significance levels α = 5% using N = 2000 Monte-Carlo replicates. Figure 3B displays the results. Intuitively, statistical power increases with δ. The multimodality test appeared the most powerful in this setting , especially when δ ≥ 3. [START_REF] Siffer | Are your data gathered[END_REF] has shown that δ = 3.04 is a threshold value above which bimodality begins to appear in such two-part Gaussian mixture. Moreover, since all clusters in the 4-cluster case were between cluster 1 and cluster 2, statistical power achieved by the multimodality test is exactly the same as in the 2-clusters case. In the two-clusters case, the direct selective test and the merging selective test had the same statistical power (since only two clusters were estimated, they were necessarily adjacent and therefore the direct selective test was exactly the same as the merging one). Meanwhile, the direct selective test failed in the 4-cluster case, regardless of the value of δ. Indeed, it was impossible to bring clusters 1 and 2 closer without mixing them with clusters 3 and 4. Fortunately, the merging selective test avoided this pitfall, because the direct selective inference test performed favourably on adjacent clusters, and carried over the separation between clusters 3 and 4. Of note, the merging selective test remains valid when the numbers of adjacent p-values increases.

Application to real ecological data

To further assess our proposed approaches, we also analyzed real data available from the R package palmerpenguins [START_REF] Allison | palmerpenguins: Palmer Archipelago (Antarctica) penguin data[END_REF]. This benchmark dataset features p = 4 measured variables -bill length (mm), bill depth (mm), flipper length (mm), and body mass (g) -for 344 penguins. After removing observations containing missing values for at least one of the 4 variables, n = 333 observations were kept in our analysis. The penguins belonged three different species: Adelie, Chinstrap, and Gentoo (with 146, 68 and 119 observations respectively).

Negative control

We initially selected only female Gentoo penguins to create a negative control dataset. Since this dataset contained only observations of the same species and sex, there should be no real differences between observations. We applied hierarchical clustering using Ward method and Euclidean distance on scaled data to build 3 clusters (scaling avoids the variable with the largest variance to dominate the clustering). Since there was no information defining any group structure in this subset, the clustering artificially created differences. Table 1 presents the p-values of each of the 3 proposed test along with the ones from the t-test for all cluster pairs along each of the four measures. Once again, the t-test identified numerous spurious associations. Meanwhile, all 3 proposed tests behaved properly by not identifying any of the four measures to be significantly separating clusters.

Full data analysis

We now included all n = 333 penguins in our analysis and analyzed the data as if we didn't know the species of the penguins. Since all 3 species were present, we want to identify which features are actually separating them. Figure 4A displays the density distribution for the four (scaled) variables across all 3 species. Adelie and Chinstrap penguins appear harder to distinguish, as they only appear to differ in bill length: Chinstrap penguins have larger bills (comparable to those of Gentoo penguins) than Adelie penguins. The Gentoo penguin species is the easiest to identify, as it is clearly different in the other 3 measures. Once again, we applied hierarchical clustering to scaled data with Euclidean distance and Ward method. Figure 4B displays the results of this clustering where we cut the dendrogram to get three clusters. Those three estimated clusters recovered the true species: cluster 2 and cluster 3 each contained only penguins of the Gentoo and Chinstrap species (respectively) while cluster 1 contained a mixture of two species (100% of the Adelie penguins plus 11 Chinstrap penguins).

Table 2 presents the p-values of each of the 3 proposed test along with the ones from the t-test for all cluster pairs along each of the four measures. Since the identified clusters corresponded to the three real penguin species, the clustering step was not expected to induce any artificial differences, and thus the t-test results can be used as reference. Only one comparison was not significant at the 5% level according to the t-test: bill depth did not separate Adelie (cluster 1) from Chinstrap penguins (cluster 3), which was visually coherent with Figure 4A. The multimodality test seemed to lack statistical power here, but inspection of the measures distribution depicted in Figure 4A showed that only a few comparison exhibited multimodality (namely cluster 1/Adelie compared to either two other clusters along flipper length. Both selective tests identified more significant differences (6/11 for the direct test and 7/11 for the merging test which is more robust when more than 2 clusters are identified). The missed separations can be explained by the lack of statistical power to detect small difference (see Supplementary Table 1 for additional details). By accounting the clustering step, our proposed test could have a reduced power compared to other tests like t-test (which is is the uniformly most powerful test [START_REF] Lehmann | Some history of optimality[END_REF]). But, this is because they are appropriately accounting for the variability and the uncertainties of the clustering step leading to results that are always valid. In this paper, we propose three new statistical tests for post-clustering inference that can be used to identify variables that separate two estimated clusters. We show that the double use of data (for clustering and for inference) and the failure to propagate the uncertainty associated with cluster estimation can lead to invalid p-values. This is particularly the case when too many clusters are estimated compared to the true underlying structure of the data. In this case, since there is no true process separating every estimated clusters, the clustering forces artificial differences between observations belonging to a common group of observations. Our three proposed tests, which take into account this clustering step and/or its possible impact on inference, give p-values that indicate a separation not induced by the clustering algorithm but emanating from the underlying data generating process, while controlling the Type-I error rate adequately. Our approaches can be used regardless of the chosen clustering algorithm and take into account many data analysis pipelines where clustering results are used post-hoc to describe and interpret clusters.

All three approaches test whether there is a separation between clusters along a given variable. The selective test is a rigorously defined test based on the concepts of selective inference, adapted from the seminal work of [START_REF] Gao | Selective inference for hierarchical clustering[END_REF] and makes a Gaussian assumption on the data. Although it tests a univariate mean difference between two clusters, it also exploits the multivariate structure of the data since the (perturbed) clustering uses all variables. The multimodality test, on the other hand, is based on the more intuitive concept of multimodality to characterise the separation of two clusters along a variable. It only relies on univariate considerations, as the separation of clusters is examined based on the distribution of each variable. Thus, unlike the selective test, which has longer computation times (dependent on the number of observations, the number of variables, and the number of Monte-Carlo simulations required to estimate p-values), it is very computationally efficient (see Supplementary Figure S5). However, this simplicity comes at the expense of a larger null hypothesis: the multimodality test requires a clear separation between clusters on the variable to work well, as it only uses the variable-level information and does not consider the entire structure of the data. Finally, since false-negative problems could occur with the selective test (particularly when the two clusters of interest are separated by other clusters), we also propose a merging method based on the aggregation of p-values. This method has the advantage of correcting these false positive problems while guaranteeing good statistical power. However, its computation cost is even greater than the selective test because this approach requires the computation of all the adjacent p-values between C k and C l .

The selective test rely on some distributional assumptions. In particular, because it uses the selective inference framework, it assumes Gaussian data to efficiently control the Type I error rate. We show that the selective test remains robust to other distributions (see Supplementary Figure S6). The multimodality test is based on the Dip Test which is a non-parametric test of unimodality. However, in practice, its p-value is computed using the Uniform distribution as the reference distribution under the null of unimodality. It could affect its statistical power but the control of the Type I error is still guaranteed.

The main limitation of our tests lies in the high dimensional setting. Due to the large number of variables and their correlation, perturbation-based approaches can fail. In our case, this result is amplified by the fact that the perturbation is only univariate. Thus, the selective inference test performs poorly in high dimension since it is exclusively based on perturbations. In addition, the calculation of the p-values is done using a Monte-Carlo approach requiring the clustering step to be repeated for each simulation and for each variable. So, if the number of variables is high, the computation time of the selective test can therefore be too long. Another problem of our approaches is related to the "signal vs. noise" ratio of the high dimension. In high dimension, a small signal repeated over a large number of variables is sufficient to create separated clusters in the high dimensional space [START_REF] Klawonn | What are clusters in high dimensions and are they difficult to find? In International workshop on clustering high-dimensional data[END_REF]. For example in a two-components gaussian mixture, for n = 100 observations, a mean difference δ = 1 repeated over p = 50 variables is enough to generate separated clusters on the first principal component of a PCA. However, the unimodality test is not powerful enough when the signal is too weak. The problem here is that the existence of clusters is only due to the repetition of the signal on a large number of variables, i.e. one has to take into account all the variables and the information they bring to explain the separation between clusters, but the unimodality test is purely univariate, being only interested in the information brought by the tested variable and this is why it lacks power in high dimension. Therefore, all the issues raised by the high-dimension constitute a natural path to apply and extend the results of our work presented here.

A Proof: Computation of the selective p-values

We want to compute the selective p-value given in (4):

p C k ,C l g ≡ P H0 |X g t η| > |x g t η| |C k , C l ∈ c(X), Π η ⊥ X g = Π η ⊥ x g
To compute (4), we have to write our data matrix X as a function of our statistic X g t η and the residual term Π η ⊥ X g where Π We also have :

η ⊥ = I n -ηη t η 2 2 Since, X g = Π η ⊥ X g + I n -Π η ⊥ X g , then : c(X) = c ([x 1 | . . . |X g | . . . |x p ]) = c ([
X t g η ⊥ Π ⊥ η X g
because Π ⊥ η is the orthogonal projection matrix onto the subspace orthogonal to span(η) [START_REF] Jewell | Testing for a change in mean after changepoint detection[END_REF][START_REF] Gao | Selective inference for hierarchical clustering[END_REF].

Finally, we have :

X g ∼ N n (µ g , σ 2 g I n ) ⇒ X t g η ∼ N (µ t g η, σ 2 g η 2 2 ) ⇒ φ g H0 ∼ N 0, σ 2 g η 2 2
Thus, the p-value (4) is equal to :

P H0 |φ g | > |X g t η| | φ g ∈ S g
with S g = {φ g : C k , C l ∈ c (X(φ g ))} and φ g H0 ∼ N 0, σ 2 g η 2

B Numerical computation of the selective p-value

Since the selective p-value given in ( 4) is intractable, in practice it is computed using a Monte-Carlo approach with important sampling resulting in the p-value described in (7). However, [START_REF] Phipson | Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn[END_REF] showed that the classical Monte-Carlo estimator of a p-value could be biased for near to zero p-values. In fact, very tiny p-values could be approximated by exactly 0 using the Monte-Carlo approach, leading to statistical hypothesis testing that does not efficiently control the type I error rate. To overcome this problem, they propose to correct the Monte-Carlo p-value by adding 1 in both the numerator and the denominator of the estimated p-value. With this correction, instead of having exactly 0 Monte-Carlo p-value, the near to zero p-values are approximated by 1 N +1 where N is the number of Monte-Carlo samples. Unfortunately, we showed using numerical simulation studies that this correction could not work with the selective p-value computed as in ( 7) for two reasons. The first one is because this p-value originally come from a conditional probability (4), so by definition, there are in fact two probabilities to compute and to correct. The second problem arrives because of the important sample approach. In fact, because under H 1 the π i in (7) are very small, adding 1 will drastically change the scale of the p-value. So, because of the important sampling approach, we need to correct our Monte-Carlo p-value by adding a constant in the same order of π.

We propose to add π = 1 N N i=1 π i 1 C k , C l ∈ c(X(ω i g )) in both the numerator and the denominator of ( 7). This correction is reasonable since for small p-value, that is under H 1 we have : i) |x t g η| is large because C k and C l are truly separated on X g ii)

N i=1 π i 1{|ω i g | > |x t g η|, C k , C l ∈ c(X(ω i g )} 0 since each π i = f1(ω i g ) f2(ω i g )
where f 1 is the density of a gaussian distribution with mean 0. Then, because ω i g ∼ N (x t g η, σ 2 η 2 ) where |x t g η| is large, f 1 is evaluated in a point that is far away of the mean, and that is why f 1 (ω i g ) 0.

So using i) and ii):

N i=1 π i 1 |ω i g | ≥ |x t g η|, C k , C l ∈ c(X(ω i g )) + π N i=1 π i 1 C k , C l ∈ c(X(ω i g )) + π = 0 + π N i=1 π i 1 C k , C l ∈ c(X(ω i g )) + π 1 N N i=1 π i 1 C k , C l ∈ c(X(ω i g )) N i=1 π i 1 C k , C l ∈ c(X(ω i g )) + 1 N N i=1 π i 1 C k , C l ∈ c(X(ω i g )) = 1 N N i=1 π i 1 C k , C l ∈ c(X(ω i g )) N i=1 π i 1 C k , C l ∈ c(X(ω i g )) 1 + 1 N = 1 N 1 + 1 N = 1 N + 1
So, by correcting our Monte-Carlo p-value by adding π we obtain the estimator proposed by [START_REF] Phipson | Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn[END_REF] for small p-value.

  G), where ρ(F, G) = sup x |F (x) -G(x)| and U is the class of unimodal distributions. Thus, the dip statistic is the distance of F to the set of unimodal functions and it measures the deviation of our distribution from unimodality. If F is unimodal, then dip(F ) = 0, and conversly if F is multimodal, then dip(F ) > 0. In practice p-values can be computed as:p Dn := P d Un ≥ D nwhere d Un is the dip statistic computed for a n-sample drawn from U[0, 1] (the standard uniform distribution), D n is the observed dip statistic, and n is the sample size. Hartigan et al.

Figure 2 :

 2 Figure 2: Validity of p-values returned by our proposed tests, comparison with t-test. panel A Data generation process. Two cases are studied: a case under a global null of no clusters in the data (No cluster) and a case with three real clusters (3 clusters). In both cases, hierarchical clustering with Ward method and Euclidean distance is used to build 3 clusters. panel B Resulting p-values for each possible cluster pair for each variable for 2000 simulations of the data.

Figure 3 :

 3 Figure 3: Statistical power of the proposed tests. panel A Data generation process: data are generated according to a univariate gaussian mixture with two components (with equal probability and variance) separated by a mean difference δ (contamination model). Two cases are studied: a case where the true numbers of clusters is estimated (2 clusters) and a case where more clusters are estimated (4 clusters). The orange dashed line represents the mean of the component. panel B Statistical power (5% level) of the proposed tests for the separation of Cluster 1 and Cluster 2 according the mean difference δ separating the two components of the gaussian mixture.

Figure 4 :

 4 Figure 4: Clustering of the 333 completely observed palmer penguins. A Distribution of each scaled variable according to the three true species of penguins. B 3 clusters are build thanks to hierarchical clustering. Colors in the dendogram represent the estimated clusters. C Clustering results reveals that estimated clusters correspond to true species

Table 1 :

 1 P-values for all cluster pair tests along each of the 4 variables from the negative control real data.

	Cluster pair tested	Selective test (direct) Merging selective test Multimodality test	t-test
	Variable tested
	Cluster 1 vs Cluster 2
	bill length	0.4082	0.4110	0.4899	0.0759
	bill depth	0.6478	0.6400	0.1478	0.4802
	flipper length	0.1160	0.1154	0.0992	0.0017*
	body mass	0.3321	0.3425	0.8320	0.0000*
	Cluster 1 vs Cluster 3
	bill length	0.1748	0.4995	0.6345	0.0001*
	bill depth	0.2914	0.3025	0.5242	0.0000*
	flipper length	0.3361	0.3206	0.6146	0.0005*
	body mass	0.3404	0.3868	0.2918	0.1190
	Cluster 2 vs Cluster 3
	bill length	0.2096	0.2120	0.9140	0.0041*
	bill depth	0.1867	0.6618	0.2376	0.0000*
	flipper length	0.2101	0.4322	0.1337	0.0000*
	body mass	0.1573	0.7967	0.6759	0.0000*

* highlights significant p-values at the α = 5% level

  x 1 | . . . |0 n | . . . |x p ] + [0 n | . . . |X g | . . . |0 n ]) = c [x 1 | . . . |0 n | . . . |x p ] + 0 n | . . . |Π ⊥ η X g + (I n -Π η ⊥ )X g | . . . |0 n = c [x 1 | . . . |0 n | . . . |x p ] + 0 n | . . . |Π ⊥ η X g + I n -I n + X g | . . . |0 n = c [x 1 | . . . |0 n | . . . |x p ] + 0 n | . . . |Π ⊥ η X g + ηη t η 2 X g | . . . |0 n = c [x 1 | . . . |0 n | . . . |x p ] + 0 n | . . . |Π ⊥ η X g + = c [x 1 | . . . |0 n | . . . |x p ] + 0 n | . . . |Π ⊥ η X g +

		ηη t
		η 2 2
	ηφ g 2 η 2	| . . . |0 n	with φ g = X g	t η
	ηφ g η 2 2 | . . . |0 ηη t x g η 2 2 + ηφ g 2 η 2 | . . . |x p

n = c x 1 | . . . |x g -

C Supplementary

 5,3,6,12,15,18,19,19.5}. panel B Statistical power at the α = 5% level of the selective test computed using 500 simulation of the data as described in panel A according to δ, the mean difference between Cluster 1 and Cluster 2. For each simulation, the selective test is applied to test the separation of Cluster 1 vs Cluster 2 and Cluster 1 vs Cluster 3 only on X 1 .

E Supplementary