Martina Sundqvist
email: ma.sundqvis@gmail.com

Julien Chiquet
email: julien.chiquet@inrae.fr

Guillem Rigaill
email: guillem.rigaill@inra.fr

Adjusting the adjusted Rand Index -A multinomial story

Keywords: Clustering, Rand Index, Multinomial distribution, Statistical Inference

The Adjusted Rand Index (ARI) is arguably one of the most popular measures for cluster comparison. The adjustment of the ARI is based on a hypergeometric distribution assumption which is unsatisfying from a modeling perspective as (i) it is not appropriate when the two clusterings are dependent, (ii) it forces the size of the clusters, and (iii) it ignores randomness of the sampling. In this work, we present a new "modified" version of the Rand Index. First, we redefine the M RI by only counting the pairs consistent by similarity and ignoring the pairs consistent by difference, increasing the interpretability of the score. Second, we base the adjusted version, M ARI, on a multinomial distribution instead of a hypergeometric distribution. The multinomial model is advantageous as it does not force the size of the clusters, properly models randomness, and is easily extended to the dependant case. We show that the ARI is biased under the multinomial model and that the difference between the ARI and M ARI can be large for small n but essentially vanish for large n, where n is the number of individuals. Finally, we provide an efficient algorithm to compute all these quantities ((A)RI and M (A)RI) by relying on a sparse representation of the contingency table in our aricode package.

Introduction

With the increasing amount of data available, development of clustering methods have become crucial in unsupervised learning to explore and find patterns in data sets. Despite the wealth of theoretical research on this subject, in practice selecting and validating a clustering is difficult. To answer these questions, one often resorts to a measure of clustering comparison: when the data is labeled, the quality of the clustering is evaluated by measuring the overlap with the original labeling; in the absence of labels, the reliability of the clustering can be assessed by evaluating its stability (see, e.g. [START_REF] Von | Clustering stability: an overview[END_REF]. This can be done by comparing several clusterings obtained by perturbing the initial data set (i.e. with resampling), or by running different clustering methods on the same data set. The idea of clustering stability is dug deeper in cluster ensembles [START_REF] Strehl | Cluster ensembles-a knowledge reuse framework for combining multiple partitions[END_REF] and its variants, which involve measures of clustering comparison in the construction of the clustering itself.

Among the many measures proposed for pairwise clustering comparisons (see [START_REF] Nguyen | Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance[END_REF], for an overview) one of the most popular is the Rand index (RI) [START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF] and its adjusted variant [START_REF] Hubert | Comparing partitions[END_REF][START_REF] Leslie | The measurement of classification agreement: An adjustment to the rand statistic for chance agreement[END_REF]. The RI is designed to estimate the probability of having a coherent pair, which is a pair for which its two observations are either in the same group in the two compared clusterings or in different groups. It is computed from the contingency table of the two classifications. However, the RI depends on the number of groups [START_REF] Leslie | The measurement of classification agreement: An adjustment to the rand statistic for chance agreement[END_REF] and is therefore difficult to interpret. To overcome this issue, the Adjusted Rand Index (in short ARI) is obtained by subtracting to the RI an estimator of its expected value obtained under the assumption of two independent clusterings.

To obtain such an estimator, a population distribution has to be assumed upon the two compared clusterings, or more specifically upon the marginals of the contingency table of the two clusterings. Considering either the clusters sizes fixed or not, the two natural hypotheses that arise are either the hypergeometric distribution or the multinomial distribution. In the literature, there is discordance as to which of these hypotheses to use.

The RI and ARI as defined by [START_REF] Robert | Measuring agreement when two observers classify people into categories not defined in advance[END_REF] and then adapted by [START_REF] Hubert | Comparing partitions[END_REF] are based on the hypergeometric distribution hypothesis. In fact, considering fixed cluster sizes makes calculations easier and the expected value of the RI deterministic. However, this is a strong assumption that is violated in all cluster studies since no clustering algorithm fixes cluster sizes (see [START_REF] Wagner | Comparing clusterings: an overview[END_REF], for a detailed discussion). Moreover, from a modeling perspective, it implicitly ignores any randomness of the sampling procedure and considers that the set of individuals that we observed is fixed. Hence under this model the (A)RI are post-hoc quantities for which no inference to a parental population can be done, which limits the interpretation exclusively to the observed data points. Assuming the marginal to be fixed certainly simplifies the calculations under the hypothesis of independence between clustering. However, modeling dependency between clusterings under this assumption is not straightforward and rather unnatural compared to the multinomial model. Yet one certainly hopes to compare clusterings that are alike or dependant.

In comparison, the multinomial model does not assume the size of the clusters to be fixed, by considering a sample observed from an infinite population. Modeling dependent clusterings and adjusting accordingly is then greatly simplified. For all these reasons we argue that the multinomial model is more natural from a statistical perspective. Note that [START_REF] Leslie | The measurement of classification agreement: An adjustment to the rand statistic for chance agreement[END_REF] already studied this model to propose an adjusted version of the RI. Nonetheless, as pointed out in [START_REF] Hubert | Comparing partitions[END_REF]; [START_REF] Steinley | Properties of the hubert-arable adjusted rand index[END_REF]; [START_REF] Steinley | A note on the expected value of the rand index[END_REF], Morey and Agresti made an error in their calculation of the expected value of the RI, assuming that the expected value of a squared variable is the square of the expected value, which is wrong in general. We are convinced that this error is the reason for the problem described in [START_REF] Steinley | A note on the expected value of the rand index[END_REF], advocating unfairly for the hypergeometric version of the (A)RI. §

In this work, we essentially make a rigorous statistical analysis of the RI under the hypothesis of a multinomial distribution. In details, our contributions are the following:

1. Define new versions of the RI and the ARI, denoted by M RI and M ARI (for "modified" (A)RI), only counting consistent pairs by similarity. Indeed, we show that counting consistent pairs by dissimilarity is unnecessary and blurs the interpretation. In terms of our newly defined M ARI, considering those pairs would simply result in a multiplication by 2. 2. Finalise the work of [START_REF] Leslie | The measurement of classification agreement: An adjustment to the rand statistic for chance agreement[END_REF] and derive an unbiased estimator of the expected value of the M RI under a multinomial distribution valid for data under H 1 (dependent clusterings) and H 0 (independent clusterings). 3. Provide an efficient algorithm to compute all these quantities ((A)RI and M (A)RI) by relying on a sparse representation of the contingency table.

The complexity is in O(n) time and space where n is the number of individuals. This is better than the usual O(n + KL) complexity, where K and L are the sizes of the two clusterings one which to compare, typically obtained when using the non-sparse contingency table. Our code is available in versions ≥ 1.0.0 of the R package aricode [START_REF] Chiquet | aricode: Efficient Computations of Standard Clustering Comparison Measures[END_REF]. 4. Investigate the difference with the hypergeometric Hubert and Arabie's ARI and show that it is biased under the multinomial distribution, even if the difference between the two estimators remains small. This is in contradiction with the results of [START_REF] Steinley | A note on the expected value of the rand index[END_REF] We introduce the two quantities c 1 ij and c 2 ij indicating whether i and j are in the same group for respectively classification C 1 and C 2 :

c 1 ij = 1 if c 1 i = c 1 j = k, 0 otherwise, and c 2 ij = 1 if c 2 i = c 2 j = , 0 otherwise.
Note that c 1 ij and c 2 ij are the realisations of Bernoulli random variables denoted by C 1 ij and C 2 ij that will prove useful later in our statistical analysis, while studying the RI and other similar quantities as random variables.

Using these two quantities we see that a pair is consistent by similarity if

c 1 ij c 2 ij = 1 and consistent by difference if (1 -c 1 ij )(1 -c 2 ij ) = 1
. Now considering all pairs, we get the following formula for the RI as defined by Rand:

RI(C 1 , C 2 ) = 1 n 2 i<j c 1 ij c 2 ij + i<j (1 -c 1 ij )(1 -c 2 ij ) = 1 + 1 n 2 2 i<j c 1 ij c 2 ij - i<j c 1 ij - i<j c 2 ij .
(1)

In Equation (1), we remark that only the product c 1 ij c 2 ij depends on the joint distribution of C 1 and C 2 : all other terms, coming exclusively from coherent pairs by difference, depend on the marginal distributions of C 1 and C 2 . These terms will thus be cancelled out in any adjusted version of the RI, correcting for what would happen if C 1 and C 2 were drawn independently. Hence, we argue that considering the consistent pairs by difference unnecessarily complicates the reasoning and the probabilistic analysis of the RI. For simplicity we thus redefine the index and refer to it as the M RI (for "modified" RI):

M RI(C 1 , C 2 ) = 1 n 2 i<j c 1 ij c 2 ij . ( 2
)
Remark 1 For the derivation of the expected value of M RI, RI and their adjusted version M ARI and ARI, using the definition involving c 1 ij and c 2 ij (or more exactly C 1 ij and C 2 ij in a probabilistic perspective) considerably simplify the calculations compared to their classical combinatorial formulations. These combinatorial formulations are recalled in the next section as they are classically used to compute the RI and its variants.

Computing the

2 . C 1 C 2 c 2 1 • • • c 2 • • • c 2 L Sums c 1 1 n 11 • • • n 1 • • • n 1L n 1. . . . . . . . . . . . . . . . . . . . . . c 1 k n k1 • • • n k • • • n 2L n k. . . . . . . . . . . . . . . . . . . . . . c 1 K n K1 • • • n K • • • n KL n K. Sums n .1 • • • n . • • • n .L k n k = n
Using basics combinatorics we get the following relations between n k , n k. , n . and c 1 ij , c 2 ij :

i<j c 1 ij = k n k. 2 , i<j c 2 ij = n . 2 and i<j c 1 ij c 2 ij = n k 2 . (3)
Expressions ( 1) and ( 2) of RI and M RI turn to

M RI(C 1 , C 2 ) = 1 n 2 k, n k 2 = 1 2 n 2 k, (n 2 k -n) (4) RI(C 1 , C 2 ) = 1 + 2 n 2 k,l n k 2 - 1 n 2 k n k. 2 + l n .l 2 . (5)
Using these formula, one can see that the minimum of the M RI is obtained when all n k are equal, which has a simple and straightforward interpretation (as two perfectly independent and balanced clusterings). On the other hand the minimum of the RI is obtained for an extremely unbalanced table, i.e.

when one of the two clustering consists of a single cluster and the other only of clusters containing single points. This makes the interpretation of the RI rather difficult (i.e. the lowest value is not obtained for two perfectly independent and balanced clusterings) and give more credibility to the definition of M RI that does not consider consistent pairs by difference.

Probabilistic model and properties of the Rand Index

So far, the (M )RI have been computed from the observed quantities c 1 ij , c 2 ij , or equivalently from the observed contingency table n k . From now, we aim to study the statistical properties of the M RI and consider its status of random variable 1 :

M RI(C 1 , C 2 ) = 1 n 2 i<j C 1 ij C 2 ij , (6)
where we recall that C 1 ij and C 2 ij are Bernoulli random variables indicating whether individual i and j are in the same groups in classification C 1 respectively C 2 .

To derive the probability of success associated to C 1 ij and C 1 ij , we need a probabilistic model for the classification of a given individual in C 1 and C 2 , that is, a counterpart for generating the two observed clusterings c 1 i and c 2 i for the n data points. We denote by C 1 i and C 2 i the corresponding random variables. A natural model is the multinomial model, which give the joint distribution of (C 1 i , C 2 i ) as follows: for all (k, ) ∈ {1, . . . , K} × {1, . . . , L},

P(C 1 i = k, C 2 i = ) = π k , s.c. K,L k, π k = 1.
The marginal probabilities of a given group is defined for 2 for a global picture. Compared to the hypergeometric model, the multinomial model easily deals with dependent classifications and does not force the size of the clusters.

k in C 1 by L π k = π k. and for in C 2 by K k π k = π . . See Table
Based on this multinomial model for C 1 i and C 2 i , it is then relatively straightforward to derive the joint distribution and marginals of C 1 ij and C 2 ij . In particular we have:

P(C 1 ij = 1) = k π 2 k. , P(C 2 ij = 1) = π 2 . and P(C 1 ij C 2 ij = 1) = k, π 2 k .
(7) However, in order to derive the expectation, variance and unbiased adjustment of the M RI under the multinomial model, one not only needs to characterize events on the classification C 1 and C 2 on (unordered) pairs of individual {i, j}, 1 By a slight abuse of notation, we use M RI for both its observed value and its definition as a random variable. We think that the context suffices for the reader to remove any ambiguity.

= P(C 1 i = k, C 2 i = ) and marginal distributions π k. = P(C 1 i = k) and π . = P(C 2 i = ) C 1 C 2 c 2 i = 1 • • • c 2 i = • • • c 2 i = L Sums c 1 i = 1 π 11 • • • π 1 • • • π 1L π 1. . . . . . . . . . . . . . . . . . . . . . c 1 i = k π k1 • • • π k • • • π 2L π k. . . . . . . . . . . . . . . . . . . . . . c 1 i = K π K1 • • • π K • • • π KL π K. Sums π .1 • • • π . • • • π .L k π k = 1
but also on pairs of pairs of individual {i, j} and {i , j }, with terms like the expectation of

C 1 ij × C 2 i j .
The following section derives a couple of technical -yet simple -lemmas, on events implying such random variables so that the final calculation of the moments of M RI under the multinomial model are straightforward.

Remark 2 To our knowledge most derivations of the expectation and variance of the RI found in the literature are based on the combinatorial formulation given in Equation ( 5): these derivations rely on general results on the moments of either the multinomial or the generalized hypergeometric distribution and involve tedious calculations. In contrast, our proofs, found in the next sections, are short, self-contained and easily accessible to any reader with some basic knowledge in probability and statistics. For this reason we argue that our proofs are interesting in their own rights.

Subsets of Pairs of Pairs -preparing the derivations of the moments of the M RI

Consider {i, j} and {i , j } the P × P set of unordered pairs of {1, . . . , n} 2 such that i < j and i < j . This set is composed by pairs of pairs, and can equivalently be seen as the set of all quadruplets of {1, . . . , n} 4 such that i < j and i < j . We partition this set into the three following subsets:

1. the unordered pairs P, 2. the ordered-triplets T , 3. the ordered quadruplets Q.

These three subsets P, T and Q makes a partition of P × P and in particular,

|P| 2 = |P| + |T | + |Q|.
We now study respectively P, T and Q in the three following lemmas: we derive their cardinality and compute some expectations involving these subsets and the C 1 ij , C 2 ij variables under the multinomial model. These three lemmas will be the building blocks for the characterization of the M RI.

Lemma 1 (Subset of unordered pairs P) With a slight abuse of notation, we consider P as a subset of P × P:

P = {{i, j, i , j } : |{i, j} ∪ {i , j }| = 2.} The cardinality of P is |P| = n 2 and E   i,j∈P C 1 ij C 2 ij   = n 2 k π 2 k . (8)
Proof. For any i, j ∈ P, we have from ( 7) that

E(C 1 ij C 2 ij ) = k π 2 k .
We just need to sum over all possible pairs to get the desired result.

Lemma 2 (Subset of ordered triplets T ) Consider the subset

T of P × P T = {{i, j, i , j } : |{i, j} ∪ {i , j }| = 3}. The cardinality of T is |T | = n(n -1)(n -2) and E T C 1 ij C 2 ij = n(n -1)(n -2) k π k π k. π . . (9)
E T C 1 ij C 2 ij C 1 ij C 2 ij = n(n -1)(n -2) k π 3 k . ( 10
)
Proof. For the cardinality of T , one can map to the set of arrangements of {1, . . . , n} 3 . For (9), remark that C 1 ij C 2 ij is a Bernoulli variable equal to 1 only when i and j are in the same cluster k in C 1 and i and j are in the same cluster in C 2 . Hence, j can be in any cluster in C 2 and j can be in any cluster k in C 1 . From here one easily get its expectation,

E(C 1 ij C 2 ij ) = k k π k π k π k = k k π k k π k π k = k π k π . π k.
and we get the desired result by summing over all triplets. For (10), remark that

C 1 ij C 2 ij C 1 ij C 2 ij is a
Bernoulli variable equal to 1 if and only if i, j and j are in the same clusters for both classifications. Summing over all T we get (10).

Lemma 3 (Subset of ordered quadruplets Q) Consider the following subset Q of P × P:

Q = {{i, j, i , j } : |{i, j} ∪ {i , j }|} = 4}.
The cardinality is |Q| = 6 n 4 and

E Q C 1 ij C 2 i ,j = 6 n 4 k, π 2 k. π 2 . . (11)
Proof. There are n 4 ways to pick 4 distinct elements of {1, ..., n} 4 . We can then arrange those in 4 2 to get an element of Q. Hence, all together there are 6 n 4 quadruplets. We get E(C 1 ij C 2 ij ) using the fact that i, j, i , j are all different and that their classes are drawn independently. We then sum over Q.

Expectation and Variance of the Rand Index

With Lemmas 1, 2 and 3, we are now equipped to easily derive the moments of the M RI. We use E for stating the expectation understood under the multinomial model in general. With the additional assumption of independence between the classification, what we refer to as the null hypothesis, we use E H0 . This terms is classically used for adjusting the Rand index.

Proposition 21 (Expectations of the M RI). Let θ denote the expectation of the M RI and θ 0 the expectation under H 0 . Then,

θ = E(M RI) = k π 2 k , θ 0 = E H0 (M RI) = k π 2 k. π 2 .
Proof. By Definition 6 and Lemma 1 we obtain θ. For θ 0 , it suffices to replace π kl by π k. π .l in the previous formula.

Similarly, we derive the expectation of the "usual" RI.

Proposition 22. Let θ RI denotes the expectation of the RI and θ RI 0 the expectation under H 0 . Then,

θ RI = E(RI) = 1 + 2 k π 2 k - k π 2 k. - π 2 . θ RI 0 = E H0 (RI) = 1 + 2 k π 2 k. π 2 . - k π 2 k. - π 2 .
Proof. Compared to the M RI, the only additional terms are 1 + i,j C 1 ij + i,j C 2 ij . Using (7) and summing over all pairs P we get the desired results.

We now continue with the variance of the M RI.

Proposition 23. Let σ 2 = V(M RI) be the variance of the M RI. Then,

σ 2 = 1 n 2    k, π 2 k -   k, π 2 k   2   + n(n -1)(n -2) n 2 2    k, π 3 k -   k, π 2 k   2    Proof.
To obtain the variance of the M RI, first rewrite the variance in terms of covariance:

σ 2 = 1 n 2 2 V P×P C 1 ij C 2 ij = 1 n 2 2 Cov P×P C 1 ij C 2 ij , P×P C 1 ij C 2 ij = 1 n 2 2 P×P Cov C 1 ij C 2 ij , C 1 i j C 2 i j
We then split this final sum using our partition of P × P. Also noticing that for all {i, j}, {i , j

} ∈ Q we have Cov(C 1 ij C 2 ij , C 1 i j C 2 i j ) = 0, we get, σ 2 = 1 n 2 2 P Cov C 1 ij C 2 ij , C 1 ij C 2 ij + T Cov C 1 ij C 2 ij , C 1 ij C 2 ij = 1 n 2 V C 1 ij C 2 i,j + n(n -1)(n -2) n 2 2 Cov C 1 ij C 2 ij , C 1 ij C 2 ij = 1 n 2    k, π 2 k -   k, π 2 k   2    + n(n -1)(n -2) n 2 2    k, π 3 k -   k, π 2 k   2   
We get the second line by enumerating the elements of P and T . We get the third line using the definition of the covariance (for any two variable X and Y : Cov(X, Y ) = E(XY ) -E(X)E(Y )) and Lemmas 1 and 2.

Remark 3 Importantly, for a fixed π k , σ 2 goes towards 0 when n grows to infinity: the larger n, the better the estimation of θ.

The Rand Index depends on the number of groups

In the multinomial model with uniform clusters (equal cluster size), [START_REF] Leslie | The measurement of classification agreement: An adjustment to the rand statistic for chance agreement[END_REF] showed that θ RI 0 depends on the number of groups in C 1 and C 2 . This is also true for M RI and easier to prove since it does not include the marginal terms of coherence by difference. We also prove the following lemma showing that if one splits a cluster of C 1 or C 2 into two, the M RI always decreases. Note that this latter lemma does not assume independence between classifications.

Lemma 4 Consider two classifications C 1 and C 2 in K + 1 respectively L clusters. Let C 1 be the classification obtained by fusing two clusters of

C 1 . Then, M RI(C 1 , C 2 ) ≤ M RI(C 1 , C 2 ).
Also, for any distribution on C 1 and C 2 we have

θ(C 1 , C 2 ) ≤ θ(C 1 , C 2 ))
Proof. Assuming without loss of generality that clusters 1 and 2 were merged, we get

M RI(C 1 , C 2 ) -M RI(C 1 , C 2 ) = 1 2 n 2 n 2 1 + n 2 2 -(n 1 + n 2 ) 2 = - 1 n 2 n 1 n 2 ≤ 0.
Since the expectation is linear, we can consider any particular model on C 1 and C 2 to get the final result.

The Adjusted version of the Rand Index

Since the (M )RI depends on the number of groups, it needs to be adjusted for chance. A way to do so, is to subtract its expectation under the null hypothesis H 0 (as motivated in [START_REF] Robert | Measuring agreement when two observers classify people into categories not defined in advance[END_REF][START_REF] Hubert | Comparing partitions[END_REF][START_REF] Leslie | The measurement of classification agreement: An adjustment to the rand statistic for chance agreement[END_REF]. Ideally one would like to get θ -θ 0 with their true values. Under our multinomial model this quantity is

θ -θ 0 = k π 2 k - k π 2 k. π 2 .
which is equal to zero under H 0 (independence of the classifications), that is, when π k. π . = π k for all k, . In practice, one can only estimate the quantities θ -θ 0 from observed classifications. Our goal is therefore to get an unbiased estimator of θ -θ 0 .

The M RI being by definition an unbiased estimator of θ, we only need an unbiased estimator of θ 0 , that is k π 2 k. π 2 . . However, under the alternative H 1 (i.e. when the compared classifications are not independent, the most natural case), deriving an unbiased estimator of θ 0 is trickier and depends on the model assumption. [START_REF] Leslie | The measurement of classification agreement: An adjustment to the rand statistic for chance agreement[END_REF] proposed a plug-in estimator for the multinomial model, but as pointed out by [START_REF] Hubert | Comparing partitions[END_REF]; [START_REF] Steinley | Properties of the hubert-arable adjusted rand index[END_REF]; [START_REF] Steinley | A note on the expected value of the rand index[END_REF], they made errors in their calculations. In the next section we continue their work by proposing an unbiased estimator for θ 0 . We also show that the hypergeometric estimator of [START_REF] Hubert | Comparing partitions[END_REF] for θ 0 , used as correction in the "traditional" ARI, is biased under our multinomial H 1 .

A new Adjusted Rand Index.We now define our own adjusted version of the M RI that we denote M ARI:

M ARI = θ -θ 0 . (12)
with

θ = P C 1 ij C 2 ij n 2 θ 0 = Q C 1 ij C 2 i j 6 n 4 .
and its observed value,

M ARI obs = P c 1 ij c 2 ij n 2 - Q c 1 ij c 2 i j 6 n 4 .
where we recall that c 1 ij and c 2 ij are the observed counterparts of C 1 ij , C 2 ij and P, Q are defined in Section 2.3.1.

Lemma 5 Under the multinomial model, the M ARI is unbiased, that is,

E(M ARI) = θ -θ 0 .
Proof. The proof is straightforward using Lemma 3.

Computing the M ARI from a contingency table.In practice, the comparison of two classifications is given as a contingency table as Table 1, and we thus need a formulation of the M ARI defined in ( 12) as a function of n k . We already gave in (4) an expression of θ as a function of n k . As we will see, θ 0 can as well be computed from the n k contingency Table 1 even if summing over all elements of Q rather than P is a bit less straightforward. To get Q c 1 ij c 2 i j , we will use the term k n 2 k. n 2 . from which we will, as a direct result of Definition (3), derive the ( P c 1 ij )( P c 2 i j ) terms. These latter can be decomposed as follows:

( P c 1 ij )( P c 2 i j ) = P c 1 ij c 2 ij + T c 1 ij c 2 ij + Q c 1 ij c 2 i j . (13)
It is then sufficient to subtract the terms of P and T from the left side of Equation ( 13) to get Q c 1 ij c 2 i j . All terms summing over P are easy to recover (see Definition 3). However, the terms involving elements of T are more tedious to obtain and are derived in Lemma 6. The terms of Q derived in Lemma 7.

Lemma 6 We have the following expression of

T c 1 ij c 2 ij in terms of n k : T c 1 ij c 2 i j = 2n + k, n k. n k n . - k, n 2 k - k n 2 k. - n 2 .
Proof. We need to consider all i in {1, ..., n}. Assuming for now that i is in classes (k, ), that is c 1 i = k and c 2 i = , let us consider all j, j such that

c 1 ij c 2 ij = 1. The term c 1 ij c 2 ij is equal to one if c 1 j = k and c 2 j = .
We then get different scenarios according to whether c 1 j = k or not and whether c 2 j = . Those scenarios are enumerated in Table 3.

Summing all terms of Table 3 we get n k. n . + 2 -n k -n k. -n . . To account for all i belonging to class (k, ) we then multiply by n k . Eventually we sum over all k, to recover

T c 1 ij c 2 ij = k, n k (2 + n k. n . -n k -n k. -n . ) = 2n + k, n k. n k n . - k, n 2 k - k n 2 k. - n 2 .
Table 3: Four scenarios to be considered for j and j in the calculation of the terms in T c 1 ij c 2 ij when i is in class (k, ).

j in j not in j in k (n k -1)(n k -2) (n k -1)(n . -n k ) j not in k (n k. -n k )(n k -1) (n k. -n k )(n . -n k ) Lemma 7 We have the following expression of Q c 1 ij c 2 i j in terms of n k : Q c 1 ij c 2 i j = k n 2 k. n 2 . -4 k n k 2 + 4(2n + k, n k. n k n . - k, n 2 k - k n 2 k. - n 2 . ) + 2n k n k. 2 + n k. 2 + n 2 4 Proof. From Equation (3) we can derive k n 2 k. n 2 . as a function of P×P c 1 ij c 2 i j and n, since, k n 2 k. = n + 2 P c 1 ij and n 2 . = n + 2 P c 2 ij with, k n 2 k. n 2 . = (2 i<j c 1 ij + n)(2 i <j c 2 i j + n) = 4 P×P c 1 ij c 2 i j + 2n P c 1 ij + P c 2 i j + n 2 (14)
Using equation ( 13), we decompose P×P c 1 ij c 2 i j into terms of P, T and Q and get,

Q c 1 ij c 2 i j = k n 2 k. n 2 . -4 P c 1 ij c 2 ij +4 T c 1 ij c 2 ij +2n P c 1 ij + P c 2 ij +n 2 4 = k n 2 k. n 2 . -4 k n k 2 +4(2n+ k, n k. n k n . - k, n 2 k - k n 2 k. - n 2 . ) + 2n k n k. 2 + n k. 2 + n 2 4
3 Implementation -package aricode

We implemented code for fast computation of the M RI and its adjusted version the M ARI, as well as a number of other clustering comparison measures in the R/C++ package aricode, which is available on CRAN. Computing these measures is straightforward by means of the whole K × L contingency table. However, the time and space complexity is in O(n + KL), which is somewhat inefficient when K and L are large. Our implementation in aricode is in O(n): the key idea is that, given n observations, at most n elements of the n k contingency matrix can be non zero. To recover these non zero elements one can proceed in two simple steps: first, all observations are sorted in lexicographical order in terms of their first and second cluster index. This can be done in O(n) using bucket sort [START_REF] Thomas | Introduction to Algorithms[END_REF] or radix sort (as implemented in R [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]). Note that once the observations are sorted, all i that are in clusters k and are one after the other in the data table. Thus, in a second step aricode counts all non zero n kl in a single path over the data table. Internally this is done using Rcpp [START_REF] Eddelbuettel | Rcpp: Seamless r and c++ integration[END_REF].

In Figure 1 we compare our implementation of the standard ARI with the implementation of mclust [START_REF] Scrucca | mclust 5: clustering, classification and density estimation using Gaussian finite mixture models[END_REF] (that uses the whole contingency table). As can be noted, the cost of the latter can be prohibitive for large vectors.

Hubert and Arabie's ARI

In this section we study the expectation of the 'standard' RI of [START_REF] Robert | Measuring agreement when two observers classify people into categories not defined in advance[END_REF] (by contrast with our M RI); the expression of which results from the hypergeometric model. This expression was used by [START_REF] Hubert | Comparing partitions[END_REF] for adjusting the RI and producing the usual ARI. We study this expected value when the expectation corresponds to the multinomial distribution. We show that this estimator is biased in general under the alternative hypothesis, that is, when the two compared clusterings are not independent.

Expectation of Hubert and Arabie's ARI

Consider the observed value of the ARI proposed by [START_REF] Robert | Measuring agreement when two observers classify people into categories not defined in advance[END_REF]; [START_REF] Hubert | Comparing partitions[END_REF]: in order to analyse this quantity in our multinomial setup, we first give its definition in terms of c 1 ij and c 2 ij , that is

ARI obs = 2 n 2 KL kl n kl 2 - 2 n 2 2 KL kl n k. 2 n .l 2 = 2 n 2 P c 1 ij c 2 ij - 2 n 2 2 P c 1 ij P c 2 ij ,
where we recall that c 1 ij ,c 2 ij are realisations of the Bernoulli variables C 1 ij , C 2 ij . In a probabilistic perspective, we consider the ARI as a random variable:

ARI = 2 n 2 P C 1 ij C 2 ij θ RI - 2 n 2 2 P C 1 ij P C 2 ij θ RI 0 , (15)
where, as for the M RI, we ignored the marginal terms in our definitions of θ RI and θ RI 0 that cancel in the ARI. We now claim the following proposition. Proposition 41. Under the multinomial model we have

E(ARI) = E( θRI ) -E( θRI 0 ), with E( θRI ) = 2 KL k π 2 k and E( θRI 0 ) = 2 n 2 2 n 2 KL k π 2 k + n(n -1)(n -2) KL k π k π k. π . + 6 n 4 KL k π 2 k. π 2 .
Assuming we are under the null this simplifies so that E H0 (ARI) = 0.

Proof. Using Lemma 1, we have

E( P C 1 ij C 2 ij ) = n 2 KL k π 2 k .
Using Definition (13) and Lemmas 1, 2, 3 we obtain

E( P C 1 ij P C 2 ij ) = n 2 KL k π 2 k +n(n-1)(n-2) KL k π k π k. π . +6 n 4 KL k π 2 k. π 2 . .
Under the null we have π 2 k = π 2 k. π 2 . and we get

E H0 ( P C 1 ij C 2 ij ) = n 2 KL k π 2 k. π 2 . E H0 ( P C 1 ij P C 2 ij ) = KL k π 2 k. π 2 . n 2 + n(n -1)(n -2) + 6 n 4 = n 2 2 KL k π 2 k. π 2 . .
The expectations E( θRI ) and E( θRI 0 ) are obtained by scaling respectively with 2/ n 2 and 2/ n 2 2 ; E(ARI) is their difference.

From these results we conclude that Hubert and Arabie's ARI is biased under the multinomial model in general, since the term used for the adjustment is biased as E( θ RI 0 ) = θ RI 0 . Note, however, that this estimator is not biased under the null H 0 .

Study of the bias Hubert and Arabie 's ARI

The quantity that we study in this section is bias

n (θ RI 0 ) = θ RI 0 -E( θ RI 0 ) = K,L k, π 2 k. π 2 . - n 2 K,L kl π 2 k + 6 n 3 K,L k π k π . π k. + 6 n 4 K,L k π 2 k. π 2 . n 2
Bias disappear when n goes to infinity. The bias can be rewritten as bias n (θ RI 0 ) = 4n -6 n(n -1)

K,L k, π 2 k. π 2 . - 2 n(n -1) K,L kl π 2 k - 4(n -2) n(n -1) K,L k π k π . π k.
From this expression we get Proof. As seen in Equation (4.2), the bias consist of three terms. The absolute value of the sum of these three terms is bounded by the sum of their absolute values. Then, using that k, π k = 1 and all π k ≥ 0, we bound k, π 2 k. π 2 . , kl π 2 k and k π k π . π k. by 1 and we get |bias n (θ RI 0 )| ≤ 4(2n-3) n(n-1) . We have, (2n -3) < 2(n -1) and the result follows.

Empirical bias. In the case of independence the bias is zero. In the case of dependence, using Lemma 8 we get that the bias is smaller than 0.04 for n larger than 200. Following the work of [START_REF] Steinley | A note on the expected value of the rand index[END_REF], we study the importance of the difference empirically for small value of n in the next paragraph. In summary for n larger than 64 we observe a small bias, typically smaller than 10 -2 . For smaller values of n the bias can be larger.

Simulation setting.We study the evolution of the bias by comparing two classifications with equal number of groups (K = L), with values varying in K ∈ {2, 4, 8, 16, 32, 64, 128} and a growing number of individuals. For drawing the two compared classifications under the multinomial model, see Table 2. We consider three scenarios described below where we tune the level of difficulty by controlling the balance between group sizes with the parameters . Scenario 1. In the first scenario we investigate a π kl distribution with a disproportionate diagonal. All other entries being null.

π k =      1 - 0 • • • 0 0 K-1 • • • 0 . . . . . . . . . . . . 0 0 • • • K-1     
Scenario 2. In the second scenario we investigate a π kl distribution with a proportional diagonal and extra diagonal dependency. All other entries being null.

π k =      (1 -)/K /K • • • 0 0 (1 -)/K • • • 0 . . . . . . . . . . . . /K 0 • • • (1 -)/K     
Scenario 3. In the third scenario we investigate a π kl distribution with one line and one column being disproportional and all other entries being null.

π k =      1 - K+L-2 • • • K+L-2 K+L-2 0 • • • 0 . . . . . . . . . . . . K+L-2 0 • • • 0     
Results. The results are shown in Figure 2 where the bias is shown in its absolute value with log 2 / log 10 scales. For the different scenarios, the parameter of imbalanceness ,is fixed to 0.3 and 0.8.

In the different scenarios, the bias remains moderates for most values of K and n. When the number of individuals is small however, the difference turns to be more important and using the (A)RI lead to misguiding conclusions.

Fig. 1 :

 1 Fig. 1: Timings comparing the cost of computing the ARI with aricode or with the commonly used function adjustedRandIndex of the mclust package.

  that used the faulty ARI of[START_REF] Leslie | The measurement of classification agreement: An adjustment to the rand statistic for chance agreement[END_REF]. The consistent pairs are all pairs where observations i and j are in the same group (consistent by similarity), or in different groups (consistent by difference) in C 1 and C 2 .

	2 Statistical Model
	2.1 A new Rand Index -counting only pairs consistent by similarity
	The Rand Index (RI) proposed by Rand (1971) counts all the consistent pairs
	in two given classifications. In details, let us consider two classifications C 1
	and C 2 in respectively K and L classes of the same n individuals. The labels
	of individual i are given by c 1 i ∈ [1, . . . , K] and c 2 i ∈ [1, . . . , L].

Table 1 :

 1 Rand Index from the n k contingency table The information from two observed classifications is usually summarized in a contingency table like Table 1, representing the number of observations n k in group k in C 1 and in group in C 2 . Contingency Table between clusterings C 1 and C 2 ; each entry n k corresponds to the number of observations in group k in C 1 and group in C

Table 2 :

 2 Multinomial model: probabilistic distributions π k

Acknowledgements This work is supported by (1) allocations doctorales sur domaines ciblés (ARDoc) de la Région île de France and (2) an ATIGE grant from Genopole. The IPS2 benefits from the support of the LabExSaclay Plant Sciences-SPS.

As a conclusion, we argue that one should always prefer our M (A)RI to the (A)RI. There are four main reasons for this.

-The adjustment of the RI is based on a hypergeometric distribution which is unsatisfying from a modeling perspective. In particular, it forces the size of the clusters to be the same and it ignores randomness of the sampling (see the introduction). A multinomial model of the M ARI does not force the size of the clusters and properly model randomness. Furthermore, the model easily extends to the dependant case. -The difference between the ARI and M ARI can be large for small n but essentially vanish for large n (see Section 4.2). -The M (A)RI can be computed just as fast as the (A)RI in only O(n) rather than O(n + KL) using our aricode package. -The M (A)RI does not take into account pairs coherent by difference which -as argued in Section 2.1 -unnecessarily complexify the analysis and interpretation of the (A)RI.

Conflict of interest

We declare that we have no conflict of interest.